speccorr.m

这里给出另一组对谱线校正的方法

%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
    %fprintf('IP=%5.6f   %5.6f   %4d\n',IP,ShiftCorr);
    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