路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
6 v3 h; o' o  s" `3 Jclear all
8 ^2 u/ ?2 ]) X- h2 Inn=0.0048
) f5 N3 l. X8 b6 {% }$ kGxn0=262144
' c; z6 d  a8 }' lN=5200
9 h0 |: v: y5 Zll=0.04
for ii=1:N/2
  nk=ii*nn
) u3 \+ u; ~4 p. `Gx(ii)=Gxn0*(nk/0.1)^(-2)
Xk(ii)=sqrt(N*Gx(ii)/2/ll)
end
9 H! Q& R1 C3 b9 @%R=2*pi*normrnd(0,1,1,N/2)
R=2*pi*rand(1,N/2)
2 n) v) M2 i% U1 Q5 n1 w  N%x2=sqrt(2)*exp(i*(3*pi/4))
! Z5 Y1 B) `4 I1 d0 [4 |$ _/ u7 jfor ii=1:N/2
: v% ^9 V" c( F$ `. A  Xkf(ii)=Xk(ii)*exp(i*R(ii))
- ?  N4 O! ~& ?6 \, q  Z   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
end
7 s  }8 |) J9 n4 ^Xkf(N/2)=0

%tt=Xk(2)*sin(R(2))
2 a5 ?9 m9 s3 ^' g( c: O0 yfor jj=1:N
    ii=(1:1:N-1)
2 }) R, w) \" v3 g. w    bb=exp(i*2*pi*ii*(jj-1)/N)
) R5 p) c$ R; e# E& @    cc=Xkf.*bb
/ t# Y+ V. ^( }- Z; R$ [    aa=sum(cc)
%aa=0
%for ii=1:N-1
& c5 @. F$ V: s9 B* {6 j% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
) L+ h8 B7 N  Y1 f%end
xmm(jj)=aa/N
7 c. b" ~% i! Q  c6 }end
%bb
! M7 S2 t: c8 n( Z% }tt=[ll:ll:N*ll]
plot(tt,xmm)
怎么用welch法生成功率谱密度，坐标为双对数坐标。

nocturne

怎么生成这个图形。

' r4 X* w7 |$ V  a5 G: i1 [0 o$ {

S2_hard

你画的图就默认把虚部去掉了，如果单纯想把虚部去掉，real（xmm）即可，但是不一定合理，是不是取模更合适，abs（xmm）

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
: i6 x/ Z1 ?# [clear all
clc
* v5 {8 A" `& x# C' O2 Zn=0.0048;
3 P! S0 B: i. s  O2 wGxn0=262144;
N=5200;
ll=0.04;
  r7 \2 {% _9 m" p! afor k=1:N/2
    nk=k*n;
    Gx(k)=Gxn0*(nk/0.1)^(-2);
5 y4 Z) l5 L' a3 q* f$ x! W4 h/ l    Xk(k)=sqrt(N*Gx(k)/2/ll);
0 [* x8 b5 a, \  n3 vend
R=2*pi*rand(1,N/2);
9 C; a8 f9 }  o; S' _for k=1:N/2
3 C* B7 z: ^  T4 x, v    Xkf(k)=Xk(k)*exp(1i*R(k));
    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
end
# M, D6 b! h! {% I" {Xkf(N/2)=0;
for j=1:N
: R' Z# }' Z- [9 s9 o5 R2 L    k=(1:1:N-1);
    bb=exp(1i*2*pi*k*(j-1)/N);
6 `* m! J: v7 f# L- w+ B, d" O6 a3 G. I    cc=Xkf.*bb;
    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
end
/ C8 D2 V' b3 q8 c+ Gt=(1:N)*ll;
' n0 }3 q2 O+ G& {figure(1)
) L; a" h3 F- v9 k/ {" L1 t% `plot(t,Xm)
  b1 J9 T5 D# q' q+ gL=length(Xm);
nfft=2^nextpow2(round(length(Xm)/4));
3 b  V; ?- v0 l" Mwindow=chebwin(nfft);
overlap=round(length(Xm)/8);
ns=ceil(N/L);
[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
figure(2)
plot(log10(f),log10(pxx))