学写程序之二－－能量重心校正法

zhaoyixu

第二贴，不知道有没有错。欢迎大家拍砖头阿，多多指教
参考自丁康老师《离散频谱校正技术》一书和版主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 编辑 ]

zhaoyixu

哎呀，居然没有一个人看，伤心了。:@L :'(

zhwang554

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 编辑 ]

terrywang12

f0=((f0fz/f0fm)-1)*Fs/N  这里为什么还要减1

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

谢谢谁帮我解释下

robustcococole

ip学习了。。。

无限剑制

不错不错，别伤心啊:lol

blue1122

不错，学习了

随风飞阳

回复 4 # terrywang12 的帖子

我也是想问你的问题，不知你找到答案了没？