学写程序－比值校正法

zhaoyixu

第一次发贴，正学习中。
编写了一个单频率谐波函数的比值校正法校正程序。
也不知道写得对不对~~~~~~~~~~~~~~~~~请多多指教~~~~~~~~~~
参考自丁康老师的《离散频谱校正技术》一书和版主yangzj的贴子，深表感谢。

close all;clear all;clc  
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*253.5453*t+240*pi/180);%被分析函数
x_hanning=5.3*cos(2*pi*253.5453*t+240*pi/180).*hanning;%加汉宁窗
y = fft(x);
yh=fft(x_hanning);
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('加矩形窗未校正');
subplot(222);stem(f,2*Ah(1:N/2));title('加汉宁窗未校正');%作图函数中的2倍是幅值恢复系数
%加矩形窗
[Amax,k]=max(A);
phmax=angle(Y(k));
if A(k-1)>A(k+1);
    deltf=1/(1+A(k)/A(k-1));
    f0=(k-1-deltf)*Fs/N    %校正后频率
    am=Amax/sinc(deltf)     %校正后幅值
    ph=(phmax+pi*deltf)*180/pi  %校正后相位
else A(k-1)<A(k+1);
    deltf=1/(1+A(k)/A(k+1));
    f0=(k-1+deltf)*Fs/N
    am=Amax/sinc(deltf)
    ph=(phmax-pi*deltf)*180/pi
end
A(k)=am;f(k)=f0;
subplot(223);stem(f,A(1:N/2));title('加矩形窗校正后');
%加汉宁窗
[Amaxh,kh]=max(Ah);
phmaxh=angle(Yh(kh));
if Ah(kh-1)>Ah(kh+1);
    deltfh=(Ah(kh)/Ah(kh-1)-2)/(1+Ah(kh)/Ah(kh-1));
    f0h=(kh-1+deltfh)*Fs/N
    amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)
    phh=(phmaxh-pi*deltfh)*180/pi
else Ah(kh-1)<Ah(kh+1);
    deltfh=(Ah(kh)/Ah(kh+1)-2)/(1+Ah(kh)/Ah(kh+1));
    f0h=(kh-1-deltfh)*Fs/N
    amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)
    phh=(phmaxh+pi*deltfh)*180/pi
end
Ah(kh)=amh;f(kh)=f0h;
subplot(224);stem(f,Ah(1:N/2));title('加汉宁窗校正后');

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

eguang8116

经过校正后的频谱图，是不是只有峰值处是校正过的，其他谱线是不变的，那么校正后的谱图就不是等间隔了吧

liyubo

为什么需要校正？

tanke19860219

那你有相位差校正法的程序吗？谢谢了啊，我现在正在做，遇到些困难啊

zhaoyixu

原帖由 tanke19860219 于 2008-6-1 15:24 发表
那你有相位差校正法的程序吗？谢谢了啊，我现在正在做，遇到些困难啊

过几天会写的，呵呵，共同学习，mgfj

zhwang554

2楼 提出  经过校正后的频谱图，是不是只有峰值处是校正过的，其他谱线是不变的
    能量中心法需从5条谱线求出频率偏离,比值法需从2条谱线求出频率偏离,它们都是从振幅谱最大值处找到频率点,再用能量中心法或比值法,所以只有峰值处是校正过的，其他谱线是不变的.
    这儿提供一种FFT和apFFT组合校正法,它可以将整个振幅校正谱,频率校正谱,相位校正谱图示出,
    图中全部谱线都校正过, 从图中可找任何频率成分的振幅校正值,频率校正值,相位校正值.

     FFT须N个样点,apFFT须2N-1个样点.FFT和apFFT组合校正法须2N-1个样点,对前N个样点作FFT,得FFT振幅谱a1及相位谱p1,对全部2N-1个样点作apFFT,得振幅谱a2及相位谱p2,

apFFT相位谱                             p2                               即初相位校正谱
FFT和apFFT相位差                  (p1-p2) /180/(1-1/N)    即频率校正谱
FFT功率谱/apFFT振幅谱          a1.^2/a2                       即振幅校正谱

close all;clc;clear all;
N=1024;
w=2*pi;
t=-N+1:N-1;
y=1.0*exp(j*(w*t*49.1/N+50.0*pi/180))+0.8*exp(j*(w*t*149.2/N+100*pi/180))+0.6*exp(j*(w*t*249.3/N+150*pi/180))+...
+0.4*exp(j*(w*t*349.4/N+200*pi/180))+0.2*exp(j*(w*t*449.5/N+250*pi/180));
y1 = y(N:end);
win =  hanning(N)';;
win1 = win/sum(win);
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);
p1 = mod(phase(y11_fft)*180/pi,360);
y2 = y(1:2*N-1);
win =  hanning(N)';;
winn =  conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);
a2 = abs(y2_fft);
p2=mod( phase(y2_fft)*180/pi,360);
      ee=(p1-p2)/180/(1-1/N);
      aa=(a1.^2)./a2;
subplot(4,1,1),stem(a2,'.');title('apFFT振幅谱');ylim([0,1]);xlim([0 N/2]);grid
subplot(4,1,2),stem(p2,'.');title('初相位校正谱');ylim([0,400]);xlim([0 N/2]);grid
subplot(4,1,3);stem(ee,'.');title('频率校正谱');ylim([-1,1]);xlim([0 N/2]);grid
subplot(4,1,4);stem(aa,'.');title('振幅校正谱');ylim([0,1.5]);xlim([0 N/2]);grid
      disp('相位校正值')
     p2(50:100:450)
      disp('频率校正值')
    ee(50:100:450)+(49:100:450)
      disp('振幅校正值')
      aa(50:100:450)

运行结果:
初相位校正值   50.0000000000024          100.000000000003         150.000000000002           200                                     249.999999999991
   频率校正值   49.0999999419078712     149.19999993660561     249.299999913428547    349.399999883880428      449.499999834899008   
   振幅校正值   0.999999755117187         0.799999777122636       0.599999946932585         0.400000089303518          0.200000165320486

[ 本帖最后由 zhwang554 于 2008-6-12 11:31 编辑 ]

zhaoyixu

http://forum.vibunion.com/forum/thread-66278-1-1.html
已发出，相位差法之时域平移法。
相位差法之改变窗长法的程序随后发出

zhaoyixu

相位差法之改变窗长法：
http://forum.vibunion.com/forum/ ... e%3D2&frombbs=1

yangzj

原帖由 zhwang554 于 2008-6-12 08:07 发表
2楼 提出  经过校正后的频谱图，是不是只有峰值处是校正过的，其他谱线是不变的
    能量中心法需从5条谱线求出频率偏离,比值法需从2条谱线求出频率偏离,它们都是从振幅谱最大值处找到频率点,再用能量中心法或比值法 ...

全相谱不校正的相位是2N-1点序列中第N点的相位,即中点的相位,而不是起始点的相位.

zhwang554

6楼程序中
t=-N+1:N-1
中间点 t=0, 即信号y中各频率成份的初相位

[ 本帖最后由 zhwang554 于 2008-6-13 21:22 编辑 ]

eguang8116

非常感谢zhwang554 和yangzj！

eguang8116

在比值程序中
amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)
这句sinc（）里为什么不是deltfh*pi？
在文献里好像都乘了pi，但是乘了以后计算结果就不对了

yangzj

原帖由 eguang8116 于 2008-7-30 18:29 发表
在比值程序中
amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)
这句sinc（）里为什么不是deltfh*pi？
在文献里好像都乘了pi，但是乘了以后计算结果就不对了

matlab里sinc(x)=sin(pi*x)/(pi*x)

eguang8116

以前一直没有注意到这个问题

gongzuobiji

我想知道    修正量是怎么计算出来的   也就是比值法的原理是什么