路面不平度与功率谱密度

nocturne

这个是我生成路面不平度的程序：
" g- `( d( h7 m* \clear all
nn=0.0048
7 P6 B/ D6 [' P3 d) GGxn0=262144
! n7 [& o2 E# n6 D4 d+ |3 pN=5200
ll=0.04
for ii=1:N/2
  W% o* A' ^! K: b4 v8 }5 Y' ~  nk=ii*nn
Gx(ii)=Gxn0*(nk/0.1)^(-2)
Xk(ii)=sqrt(N*Gx(ii)/2/ll)
end
1 }- m+ A* ]6 j7 D" C- J%R=2*pi*normrnd(0,1,1,N/2)
9 n2 U  O! R" f1 RR=2*pi*rand(1,N/2)
%x2=sqrt(2)*exp(i*(3*pi/4))
for ii=1:N/2
! @  Q$ p9 F; Q5 j8 J( \1 \: m% {( J  Xkf(ii)=Xk(ii)*exp(i*R(ii))
   Xkf(N-ii)=Xk(ii)*exp(-i*R(ii))
end
( ^- B1 M0 b1 PXkf(N/2)=0
& G3 Q' Z! t5 F' t& \) w% [7 c
%tt=Xk(2)*sin(R(2))
* d/ F, O4 l; [% V$ Gfor jj=1:N
" g# k4 I+ K: Z, \7 [6 S' D    ii=(1:1:N-1)
    bb=exp(i*2*pi*ii*(jj-1)/N)
    cc=Xkf.*bb
" r+ x8 `' q$ b: P% a    aa=sum(cc)
%aa=0
%for ii=1:N-1
% aa=aa+Xkf(ii)*exp(i*2*pi*ii*(jj-1)/N)
%end
# l9 d' @/ C* ~9 Z* d7 l) L! Qxmm(jj)=aa/N
* l/ n) F& |) z- P+ t6 r! G; a4 uend
. Q1 e. ^* W3 M0 X4 N8 e%bb
tt=[ll:ll:N*ll]
plot(tt,xmm)
怎么用welch法生成功率谱密度，坐标为双对数坐标。
  K  L1 I2 Z  B0 \

nocturne

怎么生成这个图形。

S2_hard

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

hfiwioq

感觉你的原始信号的生成有点问题，注意复数计算和for循环中i，j的使用，i，j当作复数单位是可以写成1i，1j，
clear all
" A8 q) {* G( `# r4 O3 uclc
n=0.0048;
6 q" a, n( `7 |+ w% G7 ~9 I+ hGxn0=262144;
N=5200;
  `0 o: ?+ O& V) k; Z( ^ll=0.04;
for k=1:N/2
    nk=k*n;
    Gx(k)=Gxn0*(nk/0.1)^(-2);
3 ^* V" o# S: v# K$ u3 {$ T    Xk(k)=sqrt(N*Gx(k)/2/ll);
end
R=2*pi*rand(1,N/2);
3 m0 X( u) f0 bfor k=1:N/2
    Xkf(k)=Xk(k)*exp(1i*R(k));
( x# p; X9 }/ U  _& d8 b( I! j    Xkf(N-k)=Xk(k)*exp(-1i*R(k));
end
# P0 E% D% n; T9 a' _6 L5 gXkf(N/2)=0;
! {  r7 Y1 f/ s2 u: n# e/ \* H. R! m" hfor j=1:N
    k=(1:1:N-1);
    bb=exp(1i*2*pi*k*(j-1)/N);
    cc=Xkf.*bb;
2 K, M  ~7 O1 ?1 d$ s    aa=sum(cc);
    Xm(j)=abs(aa/N);                    %%此处修改了！%%%
1 ]" P5 B* @( L: O) L* ?end
  b: |7 A$ ]# D% W9 S2 q& jt=(1:N)*ll;
, f$ Q' S& o7 c( d7 |2 [figure(1)
( V7 V# e4 Z7 b  i, cplot(t,Xm)
/ q  d8 @& T% rL=length(Xm);
nfft=2^nextpow2(round(length(Xm)/4));
; [& C- k; z3 `+ r6 ?2 n: b9 Z( }window=chebwin(nfft);
: B# @- W3 b& K0 b  T1 y9 \" m$ Uoverlap=round(length(Xm)/8);
7 v5 G! k4 V( N& l8 c" S$ s2 ons=ceil(N/L);
[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
figure(2)
plot(log10(f),log10(pxx))