谱校正方法

谱校正方法 采样间隔归一化成1T,采样长度为N.这样FFT离散谱线为0,1)iXiN（,相应的频率分辨率2/(1/)NfN. 设FFT离散谱线局部极高谱线为m(为了数学上简洁,假定从0开始,注意在MATLAB环境下数组实际操作的是从1开始),记频偏量. 我们需要使用谱线m和与之相邻一条次高谱线，记这连续两条谱线中左边一条序号为M(当次高谱线在m左侧时1Mm,反之Mm).

下面列出若干算法的计算公式 1. 加矩形窗的精确谱校正[1] iiiXUjV

111()sin()()cos()MMMMoptMMVVMUUMKUU





1211cos()sin()cos()sin()optMMoptMMKMZVUMKMZVUM













2121

cos()cos()()ZMZMMmZZ



2. 加矩形窗情形,采用解析单频模型的幅值比校正[1, 2] 11||()||||MMMXMmXX



3. 加汉宁窗情形,采用解析单频模型的幅值比校正[1, 2] 112||||()||||MMMMXXMmXX





4. 加矩形窗情形,采用解析单频模型的复比值校正[3] 11Re()MMMXMmXX





5. 加汉宁窗情形,采用解析单频模型的复比值校正[3]

112()MMMMXXMmXX





6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3] 11Re()MmMMXMmXX





11mRmmXXX,1111mmLmmmmXXXXXX



0.5)0.5)mLmR（（ 7. 加汉宁窗情形,采用解析单频模型的复合复比值校正[3] 112Re()MMmMMXXMmXX







112mmRmmXXXX,1111221mmmmLmmmmXXXXXXXX





0.5)0.5)mLmR（（ 8. 加矩形窗,Quin校正[4] 11Re()Re(),Re()Re()mmLRmmXXXX

11LRLRLR，,

00 RRLR当且其它

9. 加汉宁窗,Quin校正[4] 11Re()Re(),Re()Re()mmLRmmXXXX

212111LRLRLR，,

00 RRLR当且其它

References 1. Schoukens, J., R. Pintelon,H. Van Hamme. The interpolated fast Fourier transform: A comparative study. IEEE Transactions on Instrumentation and Measurement. 1992, 41(2): 226-232. 2. 谢明,丁康. 频谱分析的校正方法. 振动工程学报. 1994, 7(2): 172-179. 3. 陈奎孚, 王建立,张森文. 频谱校正的复比值法. 振动工程学报(已投). 2007. 4. Quinn, B.G. Estimating frequency by interpolation using Fourier coefficients. IEEE Transactions on Signal Processing. 1994, 42(5): 1264-1268.

%========================这是调用调试================== DT=1; N=1024; PHI=pi/3; Ampl=1; CiR=11.9; %cycles in record Freq=CiR/(DT*N); %frequency

TV=[0:N-1]; DatVec=Ampl*cos(Freq*TV*2*pi+PHI); FV=fft(DatVec); figure subplot(2,1,1);plot(TV,DatVec); subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));grid on

[MV,MI]=max(abs(FV)); %加矩形窗的解析校正--1 FreqShift=SpecCorr(FV,MI,N,1);

%加矩形窗的解析单频模型校正--2 FreqShift=SpecCorr(FV,MI,N,2);

%加汉宁窗的解析单频模型校正--3 HanDat=DatVec.*hanning(N,'periodic')'; FV=fft(HanDat); FreqShift=SpecCorr(FV,MI,N,3);

%加矩形窗的解析单频模型校正+复比值法--4 FV=fft(DatVec); FreqShift=SpecCorr(FV,MI,N,4);

%加汉宁窗的解析单频模型校正+复比值法--5 HanDat=DatVec.*hanning(N,'periodic')'; FV=fft(HanDat); FreqShift=SpecCorr(FV,MI,N,5);

%加矩形窗的解析单频模型校正+复比值法+左右平均--6 FV=fft(DatVec); FreqShift=SpecCorr(FV,MI,N,6);

%加汉宁窗的解析单频模型校正+复比值法+左右平均--7 HanDat=DatVec.*hanning(N,'periodic')'; FV=fft(HanDat); FreqShift=SpecCorr(FV,MI,N,7);

%加矩形窗的解析单频模型校正+Quinn算法--8 FV=fft(DatVec); FreqShift=SpecCorr(FV,MI,N,8);

%加汉宁窗的解析单频模型校正+Quinn算法--9 HanDat=DatVec.*hanning(N,'periodic')'; FV=fft(HanDat); FreqShift=SpecCorr(FV,MI,N,9);

===========这是子程序====================== %spectrum correction assemble % the sampling interval is 1 s (or unitary) %Input: SpecVec--Discrte Fourier Spectrum from FFT % PeakIdx--the peak index, noting the matrix in MatLab start from 1 % TL--the length (or the point number) of the FFT % method--correction method % output: PeakShift-- the corrected peak shifting from the peak in discrete % spectrum function PeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)

% picking up the second highest spectrum line if(abs(SpecVec(PeakIdx-1))>abs(SpecVec(PeakIdx+1))) IP=[PeakIdx-1,PeakIdx]; ShiftCorr=-1; %shift aligning with the PeakIdx else IP=[PeakIdx,PeakIdx+1]; ShiftCorr=0; %shift aligning with the PeakIdx end

II=IP(1)-1; % noting that the index of a matrix in MATLAB starts from 1, not zero if(method==1) %an analyitic solution for rectangular window U=real(SpecVec(IP)); V=imag(SpecVec(IP)); DW=2*pi/TL; KOPT=(sin(II*DW)*(V(2)-V(1))+cos(II*DW)*(U(2)-U(1)))/(U(2)-U(1)); Z=V.*(KOPT-cos((IP-1)*DW))./(sin(DW*(IP-1)))+U; Tmp1=(Z(2)*cos(DW*(II+1))-Z(1)*cos(DW*II))/(Z(2)-Z(1)); PeakPos=acos(Tmp1)/DW; PeakShift=PeakPos-(PeakIdx-1); elseif(method==2) %based on the analytical-single-tone model for rectangular window PeakShift=abs(SpecVec(IP(2)))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1)))); PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdx elseif(method==3) %based on the analytical-single-tone model for Hanning window

PeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecVec(IP(1))))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1)))); PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdx elseif(method==4) %based on the analytical-single-tone model for rectangular window with complex correction PeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1)))); PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdx elseif(method==5) %based on the analytical-single-tone model for Hanning window with complex correction PeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1))); PeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdx elseif(method==6) %based on the analytical-single-tone model for rectangular window with complex correction+average PeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1)))); MaxPeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdx