路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
clear all
nn=0.0048
0 ~5 I1 C$ r+ n2 MGxn0=262144
* m3 {" ~" Z+ Q' g; d* [; l# rN=5200
& \; Q1 o: E- e7 q7 I4 d' ~ll=0.04
# i1 N1 ]# u5 K0 \" ?( W7 s* mfor ii=1:N/2
6 h* W9 s3 a% N8 A  nk=ii*nn
  j0 z* f) v9 @! k# z: c  O" XGx(ii)=Gxn0*(nk/0.1)^(-2)
Xk(ii)=sqrt(N*Gx(ii)/2/ll)
end
' F: M- C5 K2 W; p" ]# O%R=2*pi*normrnd(0,1,1,N/2)
R=2*pi*rand(1,N/2)
: A3 b4 C* T% {+ h: L5 q( p%x2=sqrt(2)*exp(i*(3*pi/4))
for ii=1:N/2
0 Q9 o+ a/ W4 J3 u3 }9 l  Xkf(ii)=Xk(ii)*exp(i*R(ii))
   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
end
Xkf(N/2)=0

  R) j' w/ g# L, M$ t; x%tt=Xk(2)*sin(R(2))
# p, u, d7 N" W0 P- B, ~for jj=1:N
, n  [% k/ V3 ~7 a8 l$ B    ii=(1:1:N-1)
" v5 b; w! G& J- c    bb=exp(i*2*pi*ii*(jj-1)/N)
$ D5 }, t+ v3 y8 f9 E! T9 ~8 i4 A8 c    cc=Xkf.*bb
    aa=sum(cc)
% S0 k; T$ P  h, b%aa=0
8 S8 A9 Z& p, |%for ii=1:N-1
% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
%end
: t* x( w. c/ y5 T0 bxmm(jj)=aa/N
8 z7 |% e. U* ]7 t9 jend
%bb
tt=[ll:ll:N*ll]
plot(tt,xmm)
2 }9 x+ {9 m2 U7 n/ u* H, d1 I' H4 K怎么用welch法生成功率谱密度，坐标为双对数坐标。
& t, n# Y' n/ f- x+ ^0 b4 N

nocturne

怎么生成这个图形。

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
+ S3 h  y9 B0 a/ q) r9 C; D$ Qclear all
- E$ @* l" x* x+ K  U/ w" Jclc
5 L  n4 i, S  t* I  i' F+ T9 u- @n=0.0048;
Gxn0=262144;
3 h/ g. s0 i$ j% KN=5200;
ll=0.04;
% m, B! y3 I* m/ p( p+ Z( Q9 H) y5 afor k=1:N/2
3 V& b* t+ @' j( b9 {: a; f9 C    nk=k*n;
/ @* l; A+ i2 Z4 B- z    Gx(k)=Gxn0*(nk/0.1)^(-2);
    Xk(k)=sqrt(N*Gx(k)/2/ll);
! R) r3 k! G& n: d  Dend
5 t- J% a9 z" N" q' wR=2*pi*rand(1,N/2);
for k=1:N/2
    Xkf(k)=Xk(k)*exp(1i*R(k));
: Y" t' E$ y) ]4 n+ {& ^    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
end
6 A' J/ G6 k& |7 ^" y* s% u: ~Xkf(N/2)=0;
3 A7 X" Z0 e8 k6 m0 `for j=1:N
    k=(1:1:N-1);
, L# A4 T# V- v" |/ `    bb=exp(1i*2*pi*k*(j-1)/N);
    cc=Xkf.*bb;
    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
( O! E0 N  a1 i8 }( Kend
t=(1:N)*ll;
$ Q7 p. Y" Z2 v" Z7 ?3 _5 Vfigure(1)
& u' v" ], M7 G: w0 t2 W: Z9 uplot(t,Xm)
0 O( L: ~8 a5 ]* |7 d$ yL=length(Xm);
nfft=2^nextpow2(round(length(Xm)/4));
window=chebwin(nfft);
overlap=round(length(Xm)/8);
ns=ceil(N/L);
[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
! j6 s" Q, [; n! _" B# Ifigure(2)
plot(log10(f),log10(pxx))