声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2814|回复: 7

[FFT] 学写程序之二--能量重心校正法

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

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

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

x
第二贴,不知道有没有错。欢迎大家拍砖头阿,多多指教
参考自丁康老师《离散频谱校正技术》一书和版主yangzj的相关贴子。
%能量重心校正法
clear all;clc  
n=input('请输入能量重心法校正的点数')
Fs=1024;N=1024;
t =(0:N-1)/Fs;
tt=0:N-1;
hanning=0.5-0.5*cos(2*pi*tt/N);
x=5.3*cos(2*pi*123.4*t+20*pi/180);
x_hanning=5.3*cos(2*pi*123.4*t+20*pi/180).*hanning;
y = fft(x);
yh = fft(x_hanning); %后加“h”为汉宁窗
Y=y(1:N/2)/N*2;
Yh=yh(1:N/2)/N*2;
f = Fs/2*linspace(0,1,N/2);
A=abs(Y);
Ah=abs(Yh);
subplot(221);stem(f,A(1:N/2));title('加矩形窗未校正');grid on;
subplot(222);stem(f,2*Ah(1:N/2));title('加汉宁窗未校正');grid on %作图函数中的2倍是幅值恢复系数
G=A.^2;
Gh=Ah.^2;
%加矩形窗
[Gmax,k]=max(A);
phmax=angle(Y(k));
f0fz=0;
f0fm=0;
for i=-n:n
    f0fz=f0fz+(k+i)*G(k+i);
    f0fm=f0fm+G(k+i);   
end
f0=((f0fz/f0fm)-1)*Fs/N
am=sqrt(1*(f0fm))   %矩形窗的幅值恢复系数为1
ph=(phmax+pi*(k-1-f0))*180/pi
A(k)=am;f(k)=f0;
subplot(223);stem(f,A(1:N/2));title('加矩形窗校正后');grid on
%加汉宁窗
[Gmaxh,kh]=max(Ah);
phmaxh=angle(Yh(kh));
f0hfz=0;
f0hfm=0;
for i=-n:n
    f0hfz=f0hfz+(kh+i)*Gh(kh+i);
    f0hfm=f0hfm+Gh(kh+i);
    f0h=f0hfz/f0hfm;
end
f0h=(f0h-1)*Fs/N
amh=sqrt(2.667*(f0hfm))   %2.66为汉宁窗的幅值恢复系数
phh=(phmaxh+pi*(kh-1-f0h))*180/pi
Ah(kh)=amh;f(kh)=f0h;
subplot(224);stem(f,Ah(1:N/2));title('加汉宁窗校正后');grid on
运行结果:
取n=5
矩形窗结果
f0 =123.2966   
am =5.2087
ph =38.4834
汉宁窗结果
f0h =123.4000
amh =5.3003
phh =20.0018

[ 本帖最后由 zhaoyixu 于 2008-5-25 14:43 编辑 ]

评分

1

查看全部评分

回复
分享到:

使用道具 举报

 楼主| 发表于 2008-6-7 19:56 | 显示全部楼层
哎呀,居然没有一个人看,伤心了。:@L :'(
发表于 2008-6-17 08:28 | 显示全部楼层

能量重心校正法无窗时相位校正值误差大的原因

1楼程序中,信号为
x=5.3*cos(2*pi*123.4*t+20*pi/180);

运行结果:
取n=5
矩形窗结果
f0 =123.2966   
am =5.2087
ph =38.4834

汉宁窗结果
f0h =123.4000
amh =5.3003
phh =20.0018

从以上结果看, 无窗时相位ph误差太大, 应该是20度, 校正是ph=38.4834
分析原因, 开始以为程序是否有问题, 仔细分析, 程序没有问题

原因是ph除从相位谱测得phmax外,还需校正, 校正中需频率校正值f0,  无窗时f0校正不准, 123.4-123.2966=0.1034, 引起相位校正不准

加窗phh误差很小,只要将ph和phh作比较可找出原因
从程序知 ph=(phmax+pi*(k-1-f0))*180/pi,
                phh=(phmaxh+pi*(kh-1-f0h))*180/pi

其中 phmax = 1.60333906443151
         phmaxh =1.60570291030911        这二项差不多

         f0 =      123.296562035017
         f0h =    123.399989767377           这二项差0.1034

从公式知,0.1034*180=18.612度,引起18度的误差

原因是无窗时泄漏大,频率139.4偏差0.4,泄漏很大,引起校正误差
当n取74时,才较好
n =    74
f0 =   123.399153841853
am =  5.29179779620624
ph =   20.0168699868383
f0h =  123.399999999996
amh = 5.30033123962369
phh =  19.9999999132367

[ 本帖最后由 zhwang554 于 2008-6-17 10:19 编辑 ]

评分

1

查看全部评分

发表于 2008-10-16 19:09 | 显示全部楼层
f0=((f0fz/f0fm)-1)*Fs/N  这里为什么还要减1

ph=(phmax+pi*(k-1-f0))*180/pi 这项中的pi*(k-1-f0)表示什么

谢谢谁帮我解释下
发表于 2009-12-22 21:24 | 显示全部楼层
ip学习了。。。
发表于 2010-2-11 14:13 | 显示全部楼层
不错不错,别伤心啊:lol
发表于 2011-3-9 13:54 | 显示全部楼层
不错,学习了
发表于 2012-5-4 09:50 | 显示全部楼层
回复 4 # terrywang12 的帖子

我也是想问你的问题,不知你找到答案了没?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-13 15:21 , Processed in 0.092238 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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