谱校正方法

谱校正方法采样间隔归一化成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2121cos()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，,00RRLR当且其它9.加汉宁窗,Quin校正[4]11Re()Re(),Re()Re()mmLRmmXXXX212111LRLRLR，,00RRLR当且其它References1.Schoukens,J.,R.Pintelon,H.VanHamme.TheinterpolatedfastFouriertransform:Acomparativestudy.IEEETransactionsonInstrumentationandMeasurement.1992,41(2):226-232.2.谢明,丁康.频谱分析的校正方法.振动工程学报.1994,7(2):172-179.3.陈奎孚,王建立,张森文.频谱校正的复比值法.振动工程学报(已投).2007.4.Quinn,B.G.EstimatingfrequencybyinterpolationusingFouriercoefficients.IEEETransactionsonSignalProcessing.1994,42(5):1264-1268.%========================这是调用调试==================DT=1;N=1024;PHI=pi/3;Ampl=1;CiR=11.9;%cyclesinrecordFreq=CiR/(DT*N);%frequencyTV=[0:N-1];DatVec=Ampl*cos(Freq*TV*2*pi+PHI);FV=fft(DatVec);figuresubplot(2,1,1);plot(TV,DatVec);subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));gridon[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);===========这是子程序======================%spectrumcorrectionassemble%thesamplingintervalis1s(orunitary)%Input:SpecVec--DiscrteFourierSpectrumfromFFT%PeakIdx--thepeakindex,notingthematrixinMatLabstartfrom1%TL--thelength(orthepointnumber)oftheFFT%method--correctionmethod%output:PeakShift--thecorrectedpeakshiftingfromthepeakindiscrete%spectrumfunctionPeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)%pickingupthesecondhighestspectrumlineif(abs(SpecVec(PeakIdx-1))abs(SpecVec(PeakIdx+1)))IP=[PeakIdx-1,PeakIdx];ShiftCorr=-1;%shiftaligningwiththePeakIdxelseIP=[PeakIdx,PeakIdx+1];ShiftCorr=0;%shiftaligningwiththePeakIdxendII=IP(1)-1;%notingthattheindexofamatrixinMATLABstartsfrom1,notzeroif(method==1)%ananalyiticsolutionforrectangularwindowU=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)%basedontheanalytical-single-tonemodelforrectangularwindowPeakShift=abs(SpecVec(IP(2)))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr;%shiftaligningwiththePeakIdxelseif(method==3)%basedontheanalytical-single-tonemodelforHanningwindowPeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecVec(IP(1))))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr;%shiftaligningwiththePeakIdxelseif(method==4)%basedontheanalytical-single-tonemodelforrectangularwindowwithcomplexcorrectionPeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr;%shiftaligningwiththePeakIdxelseif(method==5)%basedontheanalytical-single-tonemodelforHanningwindowwithcomplexcorrectionPeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));PeakShift=real(PeakShift)+ShiftCorr;%shiftaligningwiththePeakIdxelseif(method==6)%basedontheanalytical-single-tonemodelforrectangularwindowwithcomplexcorrection+averagePeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));MaxPeakShift=PeakShift+ShiftCorr;%shiftaligningwiththePeakIdxLeftShift=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)%basedontheanalytical-single-tonemodelforHanningwindowwithcomplexcorrection+average,????PeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));MaxPeakShift=real(P