学写程序－比值校正法
第一次发贴，正学习中。编写了一个单频率谐波函数的比值校正法校正程序。
也不知道写得对不对~~~~~~~~~~~~~~~~~请多多指教~~~~~~~~~~
参考自丁康老师的《离散频谱校正技术》一书和版主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倍是幅值恢复系数
%加矩形窗
=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('加矩形窗校正后');
%加汉宁窗
=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 编辑 ] 经过校正后的频谱图，是不是只有峰值处是校正过的，其他谱线是不变的，那么校正后的谱图就不是等间隔了吧
回复 楼主 的帖子
为什么需要校正？回复 楼主 的帖子
那你有相位差校正法的程序吗？谢谢了啊，我现在正在做，遇到些困难啊 原帖由 tanke19860219 于 2008-6-1 15:24 发表那你有相位差校正法的程序吗？谢谢了啊，我现在正在做，遇到些困难啊
过几天会写的，呵呵，共同学习，mgfj
FFT和apFFT组合校正法 (回复2楼)
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)+;%构成长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();xlim();grid
subplot(4,1,2),stem(p2,'.');title('初相位校正谱');ylim();xlim();grid
subplot(4,1,3);stem(ee,'.');title('频率校正谱');ylim([-1,1]);xlim();grid
subplot(4,1,4);stem(aa,'.');title('振幅校正谱');ylim();xlim();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 编辑 ]
回复 4楼 的帖子
http://forum.vibunion.com/forum/thread-66278-1-1.html已发出，相位差法之时域平移法。
相位差法之改变窗长法的程序随后发出
回复 4楼 的帖子
相位差法之改变窗长法：http://forum.vibunion.com/forum/viewthread.php?tid=66415&extra=page%3D2&frombbs=1 原帖由 zhwang554 于 2008-6-12 08:07 发表
2楼 提出经过校正后的频谱图，是不是只有峰值处是校正过的，其他谱线是不变的
能量中心法需从5条谱线求出频率偏离,比值法需从2条谱线求出频率偏离,它们都是从振幅谱最大值处找到频率点,再用能量中心法或比值法 ...
全相谱不校正的相位是2N-1点序列中第N点的相位,即中点的相位,而不是起始点的相位.
回复 9楼
6楼程序中t=-N+1:N-1
中间点 t=0, 即信号y中各频率成份的初相位
[ 本帖最后由 zhwang554 于 2008-6-13 21:22 编辑 ] 非常感谢zhwang554 和yangzj！
还有个问题
在比值程序中amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)
这句sinc（）里为什么不是deltfh*pi？
在文献里好像都乘了pi，但是乘了以后计算结果就不对了 原帖由 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)
yangzj兄，太感谢了
以前一直没有注意到这个问题 我想知道 修正量是怎么计算出来的 也就是比值法的原理是什么
页:
[1]
2