|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
f1=6.3;
fs=400;
det=1/400;
t=0:det:1023*det;
N=length(t);
y1=sin(2*pi*f1*t);%+rand(1,N);%+200*sin(2*pi*6.3*t);
f=abs(fft(y1));
%plot(f);
[fz,Az]=jiaozheng(y1,fs)
以上是信号模型
function [fz Az]=jiaozheng(s,fs)
index=0;
N=800;
L=200;
x=s(1:800);
xf=fft(x);
A=abs(xf);
A(1)=0;
[Amax,index]=max(A);
A0=Amax*2/N;
phmax=angle(xf(index));
f=(index-1)*fs/N;
x1=s(201:1000);
[Rk,Ik,Wk]=ft(x1,index-1);
phmax1=atan2(Ik,Rk);
dph=phmax1-phmax-2*pi*L*index/N;
dph=mod(dph,2*pi);
if dph < (-1)*pi
dph=dph+2*pi;
end
if dph > pi
dph=dph-2*pi;
end
df=(dph)/(2*pi*L/N);
fz=(index+df)*fs/N;
[Rk1,Ik1,Wk1]=ft1(df);
Az=Amax*2/Wk1/N/N;
上面是校正函数
function [Rk1,Ik1,Wk1]=ft1(f)
Rk1=0.0;
Ik1=0.0;
Num=800;
win=hanning(Num);
win1=win/sum(win);
for i=1:Num
Rk1=Rk1+win1(i)*cos(2*pi*(i-1)*(f));
Ik1=Ik1+win1(i)*sin(2*pi*(i-1)*(f));
end
Ik1=-Ik1;
Wk1=+sqrt(Rk1*Rk1+Ik1*Ik1)/Num;
plot(win1);
上面是ft1函数
function [Rk,Ik,Wk]=ft(x,k)
Rk=0.0;
Ik=0.0;
Num=length(x);
for i=1:Num
Rk=Rk+x(i)*cos(2*pi*(i-1)*k/Num);
Ik=Ik+x(i)*sin(2*pi*(i-1)*k/Num);
end
Ik=-Ik;
Wk=sqrt(Rk*Rk+Ik*Ik)/Num;
上面是单点ft函数
为什么校正的幅值总是和实际相差很大呢?
我是参考《离散频谱分析校正理论与技术》这本书的4.5.2做的
如果频率是整数,则幅值为一 如果不是的话 就相差很大!! |
|