9种谱校正方法及matlab 程序代码采样间隔归一化成1T ∆=,采样长度为N .这样FFT 离散谱线为0,1)i X i N =-（,相应的频率分辨率2/(1/)N f N ωπ∆=∆=. 设FFT 离散谱线局部极高谱线为m (为了数学上简洁,假定从0开始,注意在MATLAB 环境下数组实际操作的是从1开始),记频偏量δωδω=∆. 我们需要使用谱线m 和与之相邻一条次高谱线，记这连续两条谱线中左边一条序号为M (当次高谱线在m 左侧时1M m =-,反之M m =).下面列出若干算法的δ计算公式1. 加矩形窗的精确谱校正[1]i i i X U jV =+111()sin()()cos()M M M M opt M MV V M U U M K U U ωω+++-∆+-∆=- 1211cos()sin()cos()sin()opt M M opt M M K M Z V U M K M Z V U M ωωωωωω++-∆⎡⎤=+⎢⎥∆⎣⎦-∆+∆⎡⎤=+⎢⎥∆+∆⎣⎦ 2121cos()cos()()Z M Z M M m Z Z ωωωδ∆+∆-∆=+-- 2. 加矩形窗情形,采用解析单频模型的幅值比校正[1, 2]11||()||||M M M X M m X X δ++=+-+ 3. 加汉宁窗情形,采用解析单频模型的幅值比校正[1, 2]112||||()||||M M M M X X M m X X δ++-=+-+ 4. 加矩形窗情形,采用解析单频模型的复比值校正[3]11Re ()M M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭5. 加汉宁窗情形,采用解析单频模型的复比值校正[3]112()M M M MX X M m X X δ+++=+-- 6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]11Re ()M m M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭ 11m R m mX X X δ++=-,1111m m L m m m m X X X X X X δ---=-=-- 0.5)0.5)m L m R δδδδδ=-++（（7. 加汉宁窗情形,采用解析单频模型的复合复比值校正[3]112Re ()M M m M M X X M m X X δ++⎛⎫+=+-⎪-⎝⎭ 112m m R m mX X X X δ+++=-,1111221m m m m L m m m m X X X X X X X X δ----++=-=-- 0.5)0.5)m L m R δδδδδ=-++（（8. 加矩形窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 11L R L R L Rααδδαα==---，, 00 R R L R δδδδδ>>⎧=⎨⎩当且其它9. 加汉宁窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 212111L RL R L Rααδδαα++==---，,00 R R L R δδδδδ>>⎧=⎨⎩当且其它References1. Schoukens, J., R. Pintelon,H. Van Hamme. The interpolated fast Fourier transform: Acomparative 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.IEEETransactions on Signal Processing. 1994, 42(5): 1264-1268.%========================这是调用调试==================DT=1;N=1024;PHI=pi/3;Ampl=1;CiR=11.9; %cycles in recordFreq=CiR/(DT*N); %frequencyTV=[0:N-1];DatVec=Ampl*cos(Freq*TV*2*pi+PHI);FV=fft(DatVec);figuresubplot(2,1,1);plot(TV,DatV ec);subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));grid on[MV,MI]=max(abs(FV));%加矩形窗的解析校正--1FreqShift=SpecCorr(FV,MI,N,1);%加矩形窗的解析单频模型校正--2FreqShift=SpecCorr(FV,MI,N,2);%加汉宁窗的解析单频模型校正--3HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,3);%加矩形窗的解析单频模型校正+复比值法--4FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,4);%加汉宁窗的解析单频模型校正+复比值法--5HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,5);%加矩形窗的解析单频模型校正+复比值法+左右平均--6FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,6);%加汉宁窗的解析单频模型校正+复比值法+左右平均--7HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,7);%加矩形窗的解析单频模型校正+Quinn算法--8FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,8);%加汉宁窗的解析单频模型校正+Quinn算法--9HanDat=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% spectrumfunction PeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)% picking up the second highest spectrum lineif(abs(SpecVec(PeakIdx-1))>abs(SpecVec(PeakIdx+1)))IP=[PeakIdx-1,PeakIdx];ShiftCorr=-1; %shift aligning with the PeakIdxelseIP=[PeakIdx,PeakIdx+1];ShiftCorr=0; %shift aligning with the PeakIdxendII=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 windowU=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 PeakIdxelseif(method==3) %based on the analytical-single-tone model for Hanning windowPeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecV ec(IP(1))))/(abs(SpecV ec(IP(2)))+abs(SpecVec(IP(1 ))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==4) %based on the analytical-single-tone model for rectangular window with complex correctionPeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==5) %based on the analytical-single-tone model for Hanning window with complex correctionPeakShift=(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+averagePeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));MaxPeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxLeftShift=real(SpecVec(PeakIdx)/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1)))-1;RightShift=real(SpecVec(PeakIdx+1)/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx)));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==7) %based on the analytical-single-tone model for Hanning window with complex correction+average,????PeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));MaxPeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdxLeftShift=(2*SpecVec(PeakIdx)+SpecV ec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;RightShift=(2*SpecVec(PeakIdx+1)+SpecV ec(PeakIdx))/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx ));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==8) %Quinn method for the rectangular windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = a1/(1-a1);RightShift = -a2/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endelseif(method==9) %Quinn method for the Hanning windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = (2*a1+1)/(1-a1);RightShift = -(2*a2+1)/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endendend。