路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
clear all
# U! Q, C/ a9 ]) o# d. N, m6 Wnn=0.0048
Gxn0=262144
: p+ w8 D/ [1 O" |" p' l) AN=5200
ll=0.04
for ii=1:N/2
+ I( l* j0 U. `5 W: F9 Z4 A* b% t  l  nk=ii*nn
- e, I: B7 j' O3 O3 Z( \" b% c! X- qGx(ii)=Gxn0*(nk/0.1)^(-2)
Xk(ii)=sqrt(N*Gx(ii)/2/ll)
end
1 _) W0 k/ v6 N( f! f%R=2*pi*normrnd(0,1,1,N/2)
, Q" b7 a" }5 N- ^; b% }/ k$ {R=2*pi*rand(1,N/2)
# B. ^/ J0 |# {0 H* g$ A3 d%x2=sqrt(2)*exp(i*(3*pi/4))
  Y6 c9 Z9 m& G, k8 Nfor ii=1:N/2
2 s7 p4 ?. p& @6 W3 u9 w% L7 F8 e: B  Xkf(ii)=Xk(ii)*exp(i*R(ii))
   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
1 [) t# w; z' l* q# @end
5 I& s( u* \* GXkf(N/2)=0

0 c6 f% d# ?8 r%tt=Xk(2)*sin(R(2))
for jj=1:N
    ii=(1:1:N-1)
/ q/ B( |. X2 H+ m    bb=exp(i*2*pi*ii*(jj-1)/N)
    cc=Xkf.*bb
    aa=sum(cc)
%aa=0
%for ii=1:N-1
% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
  C3 Q( [& {4 H  I%end
xmm(jj)=aa/N
end
( r! R. U: z3 N1 w$ l9 ^2 ]%bb
4 z% W& K3 V! O/ K) X  [tt=[ll:ll:N*ll]
% [" Y6 m& t) b1 hplot(tt,xmm)
, r* `. G- o  S( i1 J7 F怎么用welch法生成功率谱密度，坐标为双对数坐标。

nocturne

怎么生成这个图形。

. `9 Q+ a. {4 M: i! b, s. z0 w

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
clear all
4 b5 I) z5 X- Pclc
n=0.0048;
Gxn0=262144;
N=5200;
! r5 W  ~1 [  z( I# P0 Wll=0.04;
for k=1:N/2
* M4 C: V3 s% U2 p4 k4 p; M5 m0 [    nk=k*n;
    Gx(k)=Gxn0*(nk/0.1)^(-2);
    Xk(k)=sqrt(N*Gx(k)/2/ll);
end
R=2*pi*rand(1,N/2);
for k=1:N/2
. f3 L0 P( R" ~1 S    Xkf(k)=Xk(k)*exp(1i*R(k));
    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
4 Q$ d7 A& g8 @  q# Fend
Xkf(N/2)=0;
for j=1:N
# p% Y* z# V! C( H    k=(1:1:N-1);
    bb=exp(1i*2*pi*k*(j-1)/N);
    cc=Xkf.*bb;
    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
' ~+ r! z. t" o/ c6 Iend
- |2 u6 I* j/ d, {; w2 bt=(1:N)*ll;
figure(1)
& p( `. i8 g8 J' A- S/ ]2 V. D7 dplot(t,Xm)
( `( M! E$ j$ b, x$ q9 `3 ZL=length(Xm);
nfft=2^nextpow2(round(length(Xm)/4));
( a! l6 W% c3 Q7 e& {window=chebwin(nfft);
overlap=round(length(Xm)/8);
ns=ceil(N/L);
4 m7 O1 o* E2 f  @* S9 p[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
figure(2)
plot(log10(f),log10(pxx))