路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
" _% f- P& S- tclear all
nn=0.0048
Gxn0=262144
- B3 x2 ^& B3 s  [% ]) BN=5200
ll=0.04
for ii=1:N/2
  nk=ii*nn
& A2 C4 A. ?6 V) k$ Q' g/ K. m$ x' VGx(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)
& V( _) R8 F, H, L7 M  ]; U+ F%x2=sqrt(2)*exp(i*(3*pi/4))
( D1 _/ S' m2 u! e8 S9 `for ii=1:N/2
! b0 W5 t8 D, ~7 p- j  Xkf(ii)=Xk(ii)*exp(i*R(ii))
   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
end
Xkf(N/2)=0

) k0 ?: [% o; I  O%tt=Xk(2)*sin(R(2))
for jj=1:N
    ii=(1:1:N-1)
    bb=exp(i*2*pi*ii*(jj-1)/N)
& J, N' _: ^% ]" W    cc=Xkf.*bb
( b) d0 F3 O" g  W) i* M0 R3 k- S    aa=sum(cc)
( y3 I- m( K/ `$ `9 J# L%aa=0
%for ii=1:N-1
1 S0 `5 C( _1 C$ w; X% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
) P. p) p* k, k$ E. ]+ m; V9 V%end
xmm(jj)=aa/N
end
%bb
tt=[ll:ll:N*ll]
plot(tt,xmm)
+ G$ E. V; s* h) R8 Q怎么用welch法生成功率谱密度，坐标为双对数坐标。

nocturne

怎么生成这个图形。

& v# W5 s! |( M$ l

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
clear all
/ ^: W# i4 g5 P# o: w0 Lclc
n=0.0048;
& s0 O/ J  G9 j2 G3 h) O' A, b/ pGxn0=262144;
N=5200;
% L0 O& ~6 A7 pll=0.04;
- V; Z9 D% b& @) T6 Dfor k=1:N/2
    nk=k*n;
" O! J- n  u% u4 q2 _$ I3 W5 _7 ]    Gx(k)=Gxn0*(nk/0.1)^(-2);
    Xk(k)=sqrt(N*Gx(k)/2/ll);
end
7 V' Y) Z1 P/ Q: Y7 uR=2*pi*rand(1,N/2);
for k=1:N/2
7 D# s$ O! p. [( g2 r8 X- K, b; d    Xkf(k)=Xk(k)*exp(1i*R(k));
9 y4 m8 n6 ^. Z) t% M# p8 r    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
end
Xkf(N/2)=0;
for j=1:N
  G9 l* @: w! }- y' f    k=(1:1:N-1);
. |# a& z' M5 F1 u  w    bb=exp(1i*2*pi*k*(j-1)/N);
0 @& j, l5 D5 o9 ^    cc=Xkf.*bb;
    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
end
( E( ~/ C2 _+ W- S, `6 j# c" c+ Dt=(1:N)*ll;
9 Z) t" B2 @' N; J" ]6 V3 Cfigure(1)
* C- z( j5 R+ q5 p, }plot(t,Xm)
/ G* p% |- }+ f, _' gL=length(Xm);
nfft=2^nextpow2(round(length(Xm)/4));
window=chebwin(nfft);
overlap=round(length(Xm)/8);
ns=ceil(N/L);
; A! j, L, ^' y4 a6 z2 E7 H+ Q[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
figure(2)
- {. @# O- O& l9 ?1 M" hplot(log10(f),log10(pxx))