学写程序之二－－能量重心校正法
第二贴，不知道有没有错。欢迎大家拍砖头阿，多多指教参考自丁康老师《离散频谱校正技术》一书和版主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;
%加矩形窗
=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
%加汉宁窗
=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 编辑 ] 哎呀，居然没有一个人看，伤心了。:@L :'(
能量重心校正法无窗时相位校正值误差大的原因
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 编辑 ] f0=((f0fz/f0fm)-1)*Fs/N这里为什么还要减1
ph=(phmax+pi*(k-1-f0))*180/pi 这项中的pi*(k-1-f0)表示什么
谢谢谁帮我解释下 ip学习了。。。 不错不错，别伤心啊:lol 不错，学习了 回复 4 # terrywang12 的帖子
我也是想问你的问题，不知你找到答案了没？
页:
[1]