离散频谱校正的程序

期待着一天

Fs=1024;N=1024;L=100;
t =(0:N+L-1)/Fs;
tt=0:N-1;
hanning=0.5-0.5*cos(2*pi*tt/N);
windowtype=input('请选择加窗类型1.矩形窗2.汉宁窗');
x=10.343*cos(2*pi*298.30453*t+240*pi/180);%.*hanning;%L+N点时间序列
%取出两段时间序列
for i=1:N
    x1(i)=x(i);
    x2(i)=x(i+L);
end
if windowtype==1
    y1=fft(x1);%对第一段信号加矩形窗做FFT变换
    y2=fft(x2);%对第二段信号加矩形窗做FFT变换
elseif windowtype==2
    y1=fft(x1.*hanning);%对第一段信号加汉宁窗做FFT变换
    y2=fft(x2.*hanning);%对第二段信号加汉宁窗做FFT变换
else
    error('选择有误,请重新选择');
end
Y1=abs(y1(1:N/2)/N*2);%第一段信号幅值
Y2=abs(y2(1:N/2)/N*2);%第二段信号幅值
f=(1:N/2)*Fs/N;
subplot(211);
   if windowtype==1
        plot(f,Y1);
        xlabel('f');ylabel('A');title('加矩形窗校正前');grid on
    elseif windowtype==2
        plot(f,2*Y1);
        xlabel('f');ylabel('A');title('加汉宁窗校正前');grid on
   end
[Y1Amax,k1]=max(Y1);
[Y2Amax,k2]=max(Y2);
phase1=angle(y1(k1));%第一段峰值所对应的相位
phase2=angle(y2(k2));%第二段峰值所对应的相位
   if windowtype==1
        A_uncorrect=Y1Amax
    elseif windowtype==2
        A_uncorrect=Y1Amax*2
   end
f_uncorrect=(k1-1)*Fs/N
phase_uncorrect=phase1*180/pi
delt=mod(phase2-phase1-2*pi*(k1-1)*L/N,2*pi);
%将delt调整到（-pi,pi）之间
if delt<-pi
    delt1=delt+2*pi;
elseif delt>pi
    delt1=delt-2*pi;
else delt1=delt;
end
deltf=delt1/(2*pi*L/N);
f_correct=(k1-1+deltf)*Fs/N
phase_correct=(phase1-deltf*pi)*180/pi        
Y_correct=zeros(1,N/2);
if windowtype==1
    A_correct=Y1Amax/sinc(deltf)
    Y_correct(k1)=A_correct;
    f(k1)=f_correct;
    subplot(212);stem(f,Y_correct);grid on
    xlabel('f');ylabel('A');title('加矩形窗校正后');
elseif windowtype==2
    A_correct=2/sinc(deltf)*Y1Amax*(1-deltf^2)
    Y_correct(k2)=A_correct;
    f(k2)=f_correct;
    subplot(212);stem(f,Y_correct);grid on
    xlabel('f');ylabel('A');title('加汉宁窗校正后');
end
请帮忙看下这个程序哪个地方不对

songzy41

LZ编写的程序好象主要取自于“离散频谱校正--相位差法之改变窗长法（源程序）”：
http://forum.vibunion.com/thread-66415-1-14.html
做了小部分修改，应该没有什么大错吧。