谱校正方法

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

下面列出若干算法的计算公式

1. 加矩形窗的精确谱校正[1]

XiUijVi

Kopt(VM1VM)sin(M)(UM1UM)cos(M)

UM1UMKcos(M)Z1VMoptUMsin(M)

Koptcos(M)Z2VM1UM1sin(M)Z2cos(M)Z1cos(M)(Mm)

Z2Z12. 加矩形窗情形,采用解析单频模型的幅值比校正[1, 2]

|XM1|(Mm)

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

2|XM1||XM|(Mm)

|XM||XM1|4. 加矩形窗情形,采用解析单频模型的复比值校正[3]

ReXM1XM1XM(Mm) 5. 加汉宁窗情形,采用解析单频模型的复比值校正[3]

2XM1XM(Mm)

XM1XM6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]

XM1mReXM1XM(Mm) RXm1XmXm1,L 1Xm1XmXmXm1XmXm1（0.5m)L（0.5m)R

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

mReR2XM1XMXM1XM(Mm) 2Xm1Xm2XmXm12Xm1Xm,L 1Xm1XmXmXm1XmXm1（0.5m)L（0.5m)R

8. 加矩形窗,Quin校正[4]

LLRe(Xm1)Re(Xm1) ,RRe(Xm)Re(Xm)L，RR, 1L1RR 当R0且L0 

R 其它9. 加汉宁窗,Quin校正[4]

LLRe(Xm1)Re(Xm1) ,RRe(Xm)Re(Xm)2L121, ，RR1L1R

 当R0且L0 RR 其它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. 3. 4.

谢明,丁康. 频谱分析的校正方法. 振动工程学报. 1994, 7(2): 172-179. 陈奎孚, 王建立,张森文. 频谱校正的复比值法. 振动工程学报(已投). 2007.

Quinn, B.G. Estimating frequency by interpolation using Fourier coefficients. IEEE Transactions on Signal Processing. 1994, 42(5): 12-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

LeftShift=real(SpecVec(PeakIdx)/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1)))-1; RightShift=real(SpecVec(PeakIdx+1)/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx))); %average

PeakShift=(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 PeakIdx

LeftShift=(2*SpecVec(PeakIdx)+SpecVec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;

RightShift=(2*SpecVec(PeakIdx+1)+SpecVec(PeakIdx))/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx));

%average

PeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift; elseif(method==8) %Quinn method for the rectangular window a1 = real(SpecVec(PeakIdx-1)/SpecVec(PeakIdx)); %left a2 = real(SpecVec(PeakIdx+1)/SpecVec(PeakIdx)); %right LeftShift = a1/(1-a1); RightShift = -a2/(1-a2);

if (LeftShift>0 & RightShift>0) PeakShift = RightShift; else

PeakShift = LeftShift; end

elseif(method==9) %Quinn method for the Hanning window a1 = real(SpecVec(PeakIdx-1)/SpecVec(PeakIdx)); %left a2 = real(SpecVec(PeakIdx+1)/SpecVec(PeakIdx)); %right LeftShift = (2*a1+1)/(1-a1); RightShift = -(2*a2+1)/(1-a2); if (LeftShift>0 & RightShift>0) PeakShift = RightShift; else

PeakShift = LeftShift; end end end