TA的每日心情 | 开心 2022-1-24 15:10 |
|---|
签到天数: 1 天 [LV.1]初来乍到
|
感觉你的原始信号的生成有点问题,注意复数计算和for循环中i,j的使用,i,j当作复数单位是可以写成1i,1j,) h4 z0 X9 j Q2 d
clear all
/ ^: W# i4 g5 P# o: w0 Lclc* U8 J! y& F# Y
n=0.0048;
& s0 O/ J G9 j2 G3 h) O' A, b/ pGxn0=262144;5 {8 h. }9 V/ M5 b$ J
N=5200;
% L0 O& ~6 A7 pll=0.04;
- V; Z9 D% b& @) T6 Dfor k=1:N/2$ U- g. G8 v( b u( `
nk=k*n;
" O! J- n u% u4 q2 _$ I3 W5 _7 ] Gx(k)=Gxn0*(nk/0.1)^(-2);, ^0 U4 ? S6 R( I
Xk(k)=sqrt(N*Gx(k)/2/ll);$ A0 Z8 l, z) v2 q
end
7 V' Y) Z1 P/ Q: Y7 uR=2*pi*rand(1,N/2);- ^( z' J) |/ f R7 X5 X
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));7 B0 j" L$ d. y1 y- ]0 V* s( h
end" F5 c+ y$ n) b9 Q6 S- l
Xkf(N/2)=0;% g' D" i) ]' \
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;- n7 G& q! x5 d I$ Z9 W
aa=sum(cc);( R- g% U0 F1 o- _: J2 d
Xm(j)=abs(aa/N); %%此处修改了!%%%" S/ [% a1 {* J- g: g7 B
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);* N" e y6 R8 K0 v% W. S
nfft=2^nextpow2(round(length(Xm)/4));; X o, S/ u* m4 B8 w
window=chebwin(nfft);9 y1 K2 I+ }) ^4 q; O
overlap=round(length(Xm)/8);3 d4 y% S/ S: |2 u2 Z% p
ns=ceil(N/L);
; A! j, L, ^' y4 a6 z2 E7 H+ Q[pxx,f] = pwelch(Xm,window,overlap,nfft,ns,'oneside');) C+ n# `; m" c, S! E2 W5 o6 d
figure(2)
- {. @# O- O& l9 ?1 M" hplot(log10(f),log10(pxx)) |
|