声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2354|回复: 3

[FFT] 学写程序之四(源代码)---比值校正法之时域平移法

[复制链接]
发表于 2008-6-10 00:25 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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 编辑 ]
回复
分享到:

使用道具 举报

发表于 2008-6-10 23:42 | 显示全部楼层
题目写错了,是相位差法不是比值法
 楼主| 发表于 2008-6-11 12:08 | 显示全部楼层
原帖由 yangzj 于 2008-6-10 23:42 发表
题目写错了,是相位差法不是比值法

太不小心了貌似改不了了
发表于 2012-5-20 20:32 | 显示全部楼层
看了很有收获,谢谢
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 17:27 , Processed in 0.094036 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表