路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
clear all
nn=0.0048
; |2 p6 b' t7 o, \5 I) |Gxn0=262144
N=5200
ll=0.04
  T! y# |' T* S% B. Kfor ii=1:N/2
; ?: t) q6 X$ \( R+ x( C1 g  nk=ii*nn
7 b! [: H/ m) e. V+ s0 HGx(ii)=Gxn0*(nk/0.1)^(-2)
5 O% P9 U2 E) ^# x# nXk(ii)=sqrt(N*Gx(ii)/2/ll)
end
# |" p8 N4 V% k* q+ o& l; `%R=2*pi*normrnd(0,1,1,N/2)
R=2*pi*rand(1,N/2)
' n* }, {% L7 p3 o& \# n%x2=sqrt(2)*exp(i*(3*pi/4))
for ii=1:N/2
6 D( f# d% s/ |- B  Xkf(ii)=Xk(ii)*exp(i*R(ii))
2 o' S, M" I! S; K2 Z  Q  }   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
, w- V( u% G5 _6 h; P3 H) bend
Xkf(N/2)=0

%tt=Xk(2)*sin(R(2))
for jj=1:N
) Q3 m9 t' H& S  {9 G) Y    ii=(1:1:N-1)
# W/ F$ x* W) `! ~; u  c# m) ^; J. s    bb=exp(i*2*pi*ii*(jj-1)/N)
    cc=Xkf.*bb
. K% l+ c0 ^4 Q6 B    aa=sum(cc)
' G! j  U6 i, V( d- |6 j; c0 N1 }- D' w%aa=0
8 ]7 E/ B' y. [& S" a%for ii=1:N-1
- L# e8 `) T% U2 O, P3 x% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
%end
8 I" G& @# o& r1 Z* J. Nxmm(jj)=aa/N
end
%bb
tt=[ll:ll:N*ll]
plot(tt,xmm)
# s; I' v, G* C" A( W+ V怎么用welch法生成功率谱密度，坐标为双对数坐标。
5 u& i( K* s( ^9 O" W3 h! N. O4 c

nocturne

怎么生成这个图形。

: Q$ x0 ^; f- D0 R; @# T* h

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
1 c' b; s+ Q7 _5 V  t) Vclear all
, Y, ^% ~; t+ a' {( Qclc
( Z1 w+ n4 c5 fn=0.0048;
Gxn0=262144;
N=5200;
  T- H0 c4 r& tll=0.04;
6 r* ~+ @2 f( C) M: sfor k=1:N/2
    nk=k*n;
3 Z3 s" q0 v3 o5 R, ^$ n# J    Gx(k)=Gxn0*(nk/0.1)^(-2);
    Xk(k)=sqrt(N*Gx(k)/2/ll);
end
0 D0 b2 ?6 i- IR=2*pi*rand(1,N/2);
, S: Y/ k( O& o/ D5 O& F  e& N0 Sfor k=1:N/2
    Xkf(k)=Xk(k)*exp(1i*R(k));
+ z$ R- [; C9 T3 N% k    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
6 s( a& }! s( \" Mend
. y8 K& `8 ?% R- e1 BXkf(N/2)=0;
* l1 ^- R, p7 b1 r% _for j=1:N
    k=(1:1:N-1);
8 N# x4 {- q" W    bb=exp(1i*2*pi*k*(j-1)/N);
    cc=Xkf.*bb;
    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
, U2 C; m" P3 o: Q8 x: jend
3 s8 q; y( k" A+ Z  c6 ut=(1:N)*ll;
figure(1)
) O% J& u2 E% A' k5 aplot(t,Xm)
L=length(Xm);
nfft=2^nextpow2(round(length(Xm)/4));
window=chebwin(nfft);
- P7 z" @# z8 a# n2 |9 [& ?overlap=round(length(Xm)/8);
ns=ceil(N/L);
0 o' f6 G. d, R1 D8 x0 w$ c2 h[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
figure(2)
; ~( d" P  u/ i5 xplot(log10(f),log10(pxx))