|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
第四贴,欢迎大家拍砖头啊,还请各位多多指教。参考自丁康老师《离散频谱校正技术》一书和版主yangzj的相关贴子。
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Fs:采样频率 N:做谱点数 L:平移点数
clear all;clc
Fs=1024;N=1024;L=100;
t =(0:N+L-1)/Fs;
tt=0:N-1;
hanning=0.5-0.5*cos(2*pi*tt/N);
windowtype=input('请选择加窗类型1.矩形窗2.汉宁窗');
x=10.343*cos(2*pi*298.30453*t+240*pi/180);%.*hanning;%L+N点时间序列
%取出两段时间序列
for i=1:N
x1(i)=x(i);
x2(i)=x(i+L);
end
if windowtype==1
y1=fft(x1);%对第一段信号加矩形窗做FFT变换
y2=fft(x2);%对第二段信号加矩形窗做FFT变换
elseif windowtype==2
y1=fft(x1.*hanning);%对第一段信号加汉宁窗做FFT变换
y2=fft(x2.*hanning);%对第二段信号加汉宁窗做FFT变换
else
error('选择有误,请重新选择');
end
Y1=abs(y1(1:N/2)/N*2);%第一段信号幅值
Y2=abs(y2(1:N/2)/N*2);%第二段信号幅值
f=(1:N/2)*Fs/N;
subplot(211);
if windowtype==1
plot(f,Y1);
xlabel('f');ylabel('A');title('加矩形窗校正前');grid on
elseif windowtype==2
plot(f,2*Y1);
xlabel('f');ylabel('A');title('加汉宁窗校正前');grid on
end
[Y1Amax,k1]=max(Y1);
[Y2Amax,k2]=max(Y2);
phase1=angle(y1(k1));%第一段峰值所对应的相位
phase2=angle(y2(k2));%第二段峰值所对应的相位
if windowtype==1
A_uncorrect=Y1Amax
elseif windowtype==2
A_uncorrect=Y1Amax*2
end
f_uncorrect=(k1-1)*Fs/N
phase_uncorrect=phase1*180/pi
delt=mod(phase2-phase1-2*pi*(k1-1)*L/N,2*pi);
%将delt调整到(-pi,pi)之间
if delt<-pi
delt1=delt+2*pi;
elseif delt>pi
delt1=delt-2*pi;
else delt1=delt;
end
deltf=delt1/(2*pi*L/N);
f_correct=(k1-1+deltf)*Fs/N
phase_correct=(phase1-deltf*pi)*180/pi
Y_correct=zeros(1,N/2);
if windowtype==1
A_correct=Y1Amax/sinc(deltf)
Y_correct(k1)=A_correct;
f(k1)=f_correct;
subplot(212);stem(f,Y_correct);grid on
xlabel('f');ylabel('A');title('加矩形窗校正后');
elseif windowtype==2
A_correct=2/sinc(deltf)*Y1Amax*(1-deltf^2)
Y_correct(k2)=A_correct;
f(k2)=f_correct;
subplot(212);stem(f,Y_correct);grid on
xlabel('f');ylabel('A');title('加汉宁窗校正后');
end
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
运行结果:
理论值:
幅值:10.343 频率:298.30453 相位:240
加矩形窗
校正前 校正后
A_uncorrect =8.830978191997250 A_correct =10.362781259858098
f_uncorrect =298 f_correct =2.983068229594135e+002
phase_uncorrect =-65.283647874590685 phase_correct =-1.205117805690294e+002
加汉宁窗
校正前 校正后
A_uncorrect =9.739025234720254 A_correct =10.342999982203294
f_uncorrect =298 f_correct =2.983045299993504e+002
phase_uncorrect =-65.184599901407665 phase_correct =-1.199999997844760e+002
[ 本帖最后由 zhaoyixu 于 2008-6-10 00:28 编辑 ] |
|