路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
* I# o  n( }: T$ J/ r$ qclear all
nn=0.0048
Gxn0=262144
N=5200
ll=0.04
for ii=1:N/2
  nk=ii*nn
Gx(ii)=Gxn0*(nk/0.1)^(-2)
- c/ g3 Q2 W( X4 pXk(ii)=sqrt(N*Gx(ii)/2/ll)
end
! u5 l1 V2 `! _7 e- P- \%R=2*pi*normrnd(0,1,1,N/2)
R=2*pi*rand(1,N/2)
%x2=sqrt(2)*exp(i*(3*pi/4))
for ii=1:N/2
  X  X/ Q0 @& h" m* U7 ~. Z' G. W  Xkf(ii)=Xk(ii)*exp(i*R(ii))
   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
1 ^) e: R3 F2 y0 d5 kend
Xkf(N/2)=0
5 V4 V, Y* y0 ^
2 }+ a! C& J  z& j6 u  X%tt=Xk(2)*sin(R(2))
- P! q4 L  A* C! t; r1 {7 Hfor jj=1:N
- k" L9 Q) R& \8 D8 C# G3 m+ @# y1 o    ii=(1:1:N-1)
    bb=exp(i*2*pi*ii*(jj-1)/N)
    cc=Xkf.*bb
* b! r* N  f: [0 C! h( \" s4 v    aa=sum(cc)
3 K* B3 H. i# u( D9 p%aa=0
%for ii=1:N-1
4 O4 ?9 U  ^" n( L4 ~- u0 z1 J% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
%end
xmm(jj)=aa/N
, Y" a0 a3 O6 V0 D* g7 i2 q3 kend
%bb
& _) t- m; S! [# Y3 htt=[ll:ll:N*ll]
; N: Y# Z, k7 U4 P# }* ]' y; ^plot(tt,xmm)
# T! b  e/ ~2 G怎么用welch法生成功率谱密度，坐标为双对数坐标。

nocturne

怎么生成这个图形。

* t- e' ~- p4 C3 J: b

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
4 P7 D  G7 @( T; Y" J* A* zclear all
clc
9 H4 w/ N- x$ O9 M4 rn=0.0048;
Gxn0=262144;
N=5200;
" @" t3 M" O6 a" ]7 w  ]0 o0 @ll=0.04;
for k=1:N/2
    nk=k*n;
$ M* w7 u' {" U: v& I* H' J    Gx(k)=Gxn0*(nk/0.1)^(-2);
$ S& `  A& R+ s& W! }    Xk(k)=sqrt(N*Gx(k)/2/ll);
! z3 x5 B5 y. G# i) T# eend
R=2*pi*rand(1,N/2);
for k=1:N/2
. Z" u  a8 T4 e; m    Xkf(k)=Xk(k)*exp(1i*R(k));
    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
end
Xkf(N/2)=0;
+ m4 Q- g' W; y( @4 }+ `( m. u# lfor j=1:N
    k=(1:1:N-1);
4 T& T1 P8 ~5 ~5 U; t' K    bb=exp(1i*2*pi*k*(j-1)/N);
! ?8 ]% t/ t4 V    cc=Xkf.*bb;
, g  |" Z$ `0 f) C    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
: u2 W/ [: A& q) j3 Kend
t=(1:N)*ll;
figure(1)
6 E9 y3 L$ V6 t6 D$ jplot(t,Xm)
L=length(Xm);
& y; [0 F- A$ m, x2 Unfft=2^nextpow2(round(length(Xm)/4));
8 x' M- m# Q! r+ a/ L% ]6 wwindow=chebwin(nfft);
overlap=round(length(Xm)/8);
ns=ceil(N/L);
[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
" C, p8 C; I' l7 `" vfigure(2)
plot(log10(f),log10(pxx))