离散频谱校正的程序
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
=max(Y1);
=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
请帮忙看下这个程序哪个地方不对 LZ编写的程序好象主要取自于“离散频谱校正--相位差法之改变窗长法（源程序）”：
http://forum.vibunion.com/thread-66415-1-14.html
做了小部分修改，应该没有什么大错吧。
页:
[1]