路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
clear all
nn=0.0048
- d8 m8 w# b: e6 [+ f2 nGxn0=262144
N=5200
" o- ]" R  h9 J& J5 ell=0.04
for ii=1:N/2
! S, ^) b- \' d$ J  nk=ii*nn
Gx(ii)=Gxn0*(nk/0.1)^(-2)
Xk(ii)=sqrt(N*Gx(ii)/2/ll)
end
%R=2*pi*normrnd(0,1,1,N/2)
) P9 ?1 X0 h. dR=2*pi*rand(1,N/2)
# e/ G6 {  V5 D6 J2 t- x%x2=sqrt(2)*exp(i*(3*pi/4))
! g; G9 s! o" i3 p; C. xfor ii=1:N/2
  Xkf(ii)=Xk(ii)*exp(i*R(ii))
- I0 p! Y% d, D1 f9 l   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
end
Xkf(N/2)=0
# c  O8 |) B9 U
9 t: R+ [1 k3 V9 O, o; u+ S%tt=Xk(2)*sin(R(2))
! ~/ r7 e, E9 B0 m' y  xfor jj=1:N
- r. E0 g# r5 n8 R5 L/ D* Z    ii=(1:1:N-1)
' l# n9 ^+ L. |4 v+ `2 `) @2 ^; A$ e    bb=exp(i*2*pi*ii*(jj-1)/N)
    cc=Xkf.*bb
    aa=sum(cc)
5 \3 s5 @% E: H% G%aa=0
%for ii=1:N-1
* s/ \9 E1 E+ R$ C1 I  |9 W% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
%end
xmm(jj)=aa/N
end
%bb
tt=[ll:ll:N*ll]
plot(tt,xmm)
$ z& h9 u1 E+ _, ?6 y" O- f怎么用welch法生成功率谱密度，坐标为双对数坐标。

nocturne

怎么生成这个图形。

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
' {- i* K: R$ ^6 D$ l4 }clear all
clc
n=0.0048;
7 h9 B0 h4 U" R& D9 C* D, EGxn0=262144;
N=5200;
ll=0.04;
for k=1:N/2
3 ?3 J0 A8 p! a' r0 f; W/ `% k    nk=k*n;
* t$ K( \/ C! }) q    Gx(k)=Gxn0*(nk/0.1)^(-2);
5 ]: s8 A! o  q    Xk(k)=sqrt(N*Gx(k)/2/ll);
+ u& p0 Z' J# v% n+ Q, Aend
R=2*pi*rand(1,N/2);
for k=1:N/2
    Xkf(k)=Xk(k)*exp(1i*R(k));
' O, G2 K7 Y# A    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
end
Xkf(N/2)=0;
for j=1:N
/ j7 ?# W0 z. \) R, m* N) V    k=(1:1:N-1);
    bb=exp(1i*2*pi*k*(j-1)/N);
" f, s; d/ j& s$ P0 z. K, g    cc=Xkf.*bb;
  A1 D4 v& N4 K+ Y7 [5 y  Z1 V) E2 ?    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
end
t=(1:N)*ll;
figure(1)
! J. @  F( [/ [& G3 \8 ^plot(t,Xm)
! c8 v/ T. i: ]: N9 P5 Q1 @L=length(Xm);
3 c. x' l2 }. S" G# j+ s2 ^. anfft=2^nextpow2(round(length(Xm)/4));
2 `) ], |) ~; H! ^' xwindow=chebwin(nfft);
+ c- c6 [# f, @5 \: ?4 Z- soverlap=round(length(Xm)/8);
ns=ceil(N/L);
[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
3 R* a( H! c+ ]figure(2)
plot(log10(f),log10(pxx))