路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
; f: \: \5 l( u7 p* Vclear all
  G. o9 Q: s5 k7 o% snn=0.0048
Gxn0=262144
: N" ?3 m( h$ U; ?  P/ t# l3 }N=5200
( y/ v  d) F& v% b, G4 m8 q8 S1 V4 tll=0.04
for ii=1:N/2
) @  z5 o8 N% x0 x; C& r" _- c  nk=ii*nn
6 V; ~% [! I; B8 O- i2 \1 rGx(ii)=Gxn0*(nk/0.1)^(-2)
Xk(ii)=sqrt(N*Gx(ii)/2/ll)
end
%R=2*pi*normrnd(0,1,1,N/2)
R=2*pi*rand(1,N/2)
%x2=sqrt(2)*exp(i*(3*pi/4))
9 @3 a  u( {2 c+ T" A, p6 ofor ii=1:N/2
  Xkf(ii)=Xk(ii)*exp(i*R(ii))
& X! T; g. b' v   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
end
" t* O$ Y# W% ?% P, FXkf(N/2)=0
; N% |1 l1 |7 L5 A$ Y& E% G" F' T3 k
$ V# V  s! G$ w& m7 t, g: f/ [%tt=Xk(2)*sin(R(2))
6 Q) P( M3 p- V. H' }( S1 Hfor jj=1:N
/ }: c$ ^. ^/ g% b2 {9 G" a    ii=(1:1:N-1)
9 b1 _+ h5 y* A6 J' l, g    bb=exp(i*2*pi*ii*(jj-1)/N)
    cc=Xkf.*bb
    aa=sum(cc)
: w$ X9 j, c  A. h% V  S; I%aa=0
%for ii=1:N-1
% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
) ?, M3 D$ q% s9 O! v%end
/ h1 b& O# a& [0 C3 E9 g4 x( pxmm(jj)=aa/N
; G9 o4 u6 ~# x7 ?* Tend
! \" e! d4 v% w; {3 r4 r%bb
tt=[ll:ll:N*ll]
plot(tt,xmm)
怎么用welch法生成功率谱密度，坐标为双对数坐标。
% g. u4 a, S% p* H1 J

nocturne

怎么生成这个图形。

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
; G2 L* M$ N. g) J5 D2 zclear all
& {' B; Z" }. o1 a6 w. Hclc
5 H# q8 W* v5 R& F+ r( N  s" en=0.0048;
Gxn0=262144;
3 z6 x: f+ @% w3 L$ h( n4 aN=5200;
6 \) T  F% Q7 Y' a% s& h( |ll=0.04;
& e. S: c: M& `, ~  R* U5 @) b4 a# Y1 Ufor k=1:N/2
# |# _% ]1 i+ J8 Z+ b. ?5 x    nk=k*n;
    Gx(k)=Gxn0*(nk/0.1)^(-2);
    Xk(k)=sqrt(N*Gx(k)/2/ll);
5 ^8 A( t, j' s& [  U4 y7 nend
1 F4 r4 V3 R5 X) D! qR=2*pi*rand(1,N/2);
for k=1:N/2
    Xkf(k)=Xk(k)*exp(1i*R(k));
: R) Z5 _& V+ I% o; i* f    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
% r( G# h5 T. j& a6 eend
Xkf(N/2)=0;
$ k+ B. J( [3 X7 r% R' gfor j=1:N
    k=(1:1:N-1);
    bb=exp(1i*2*pi*k*(j-1)/N);
1 A% R$ n8 Z! N8 O( u5 j    cc=Xkf.*bb;
    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
end
t=(1:N)*ll;
figure(1)
plot(t,Xm)
. n% L  J! [- [6 b' wL=length(Xm);
! I& n; s8 A: R# G4 m$ hnfft=2^nextpow2(round(length(Xm)/4));
window=chebwin(nfft);
overlap=round(length(Xm)/8);
) e% K# q) V' W& Nns=ceil(N/L);
( b8 e- W% y" Z; Z! c0 M[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
4 _& {1 \2 v5 S6 A9 bfigure(2)
plot(log10(f),log10(pxx))