TA的每日心情 | 开心 2022-1-24 15:10 |
|---|
签到天数: 1 天 [LV.1]初来乍到
|
感觉你的原始信号的生成有点问题,注意复数计算和for循环中i,j的使用,i,j当作复数单位是可以写成1i,1j,
' {- i* K: R$ ^6 D$ l4 }clear all* f7 H' h8 l5 l8 V1 Q. }2 g" R
clc! s) E5 v0 K$ h V
n=0.0048;
7 h9 B0 h4 U" R& D9 C* D, EGxn0=262144;" z" d/ I2 g! G
N=5200;8 l- l, B) n3 k# q
ll=0.04;# \& O! d- z6 Y7 a5 L& w2 ~/ U0 K
for k=1:N/2
3 ?3 J0 A8 p! a' r0 f; W/ `% k nk=k*n;
* t$ K( \/ C! }) q Gx(k)=Gxn0*(nk/0.1)^(-2);
5 ]: s8 A! o q Xk(k)=sqrt(N*Gx(k)/2/ll);
+ u& p0 Z' J# v% n+ Q, Aend: d3 \! k) A, W* L$ N+ @, T
R=2*pi*rand(1,N/2);7 y, L9 j R" \1 W
for k=1:N/2' R E# [- V5 b6 |- E' N
Xkf(k)=Xk(k)*exp(1i*R(k));
' O, G2 K7 Y# A Xkf(N-k)=Xk(k)*exp(-1i*R(k));" R1 q( z$ d* V) s$ N0 [2 V! `
end/ j% l; r0 Y4 l$ W- G6 ]
Xkf(N/2)=0;4 Y# g; G5 c& h1 ~% A9 Q. w
for j=1:N
/ j7 ?# W0 z. \) R, m* N) V k=(1:1:N-1);( K7 G8 B" i1 l: U1 r" X, k& ]& T
bb=exp(1i*2*pi*k*(j-1)/N);
" f, s; d/ j& s$ P0 z. K, g cc=Xkf.*bb;
A1 D4 v& N4 K+ Y7 [5 y Z1 V) E2 ? aa=sum(cc);9 A( b& q; O$ v. m: Y2 j
Xm(j)=abs(aa/N); %%此处修改了!%%%% K1 q$ A8 u- T: p" y- \' K( r
end) V5 z1 _) y, J2 {4 [! T
t=(1:N)*ll;4 M I5 N' v$ l7 M/ J' ~2 ?* y
figure(1)
! J. @ F( [/ [& G3 \8 ^plot(t,Xm)
! c8 v/ T. i: ]: N9 P5 Q1 @L=length(Xm);
3 c. x' l2 }. S" G# j+ s2 ^. anfft=2^nextpow2(round(length(Xm)/4));
2 `) ], |) ~; H! ^' xwindow=chebwin(nfft);
+ c- c6 [# f, @5 \: ?4 Z- soverlap=round(length(Xm)/8);, X; b* r- b1 ~- d, Y8 J) R* L
ns=ceil(N/L);, x) i x, ]" c+ N" Q/ Q& M
[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');
3 R* a( H! c+ ]figure(2)3 g5 D2 ^8 [/ I& d% e
plot(log10(f),log10(pxx)) |
|