王老师您好！我在http://home.vibunion.com/blog-62061-17741.html
这里看到您的图1，我想把程序写出来，但是有两个问题：
1、（h）和（j）两个子图我画不出来，不知道怎么编写那个插值公式？
2、我本只想一次性显示3个点的频率、相位幅值，但是我的程序确实一次1个点，一次2个点，一次3个点这样显示，其实我想要的是最后一种，怎么避免前两种的出现呢？
下面是我的程序；
close all;clc;clear all;
N=128;
t=-N+1:N-1;
f=[9.5 29.3 49.1];
pha=[200 200 200];
A=[1.0 0.6 0.2];
y=1.0*cos(2*pi*t*9.5/N+200*pi/180)+0.6*cos(2*pi*t*29.3/N+200*pi/180)+0.2*cos(2*pi*t*49.1/N+200*pi/180);
y1 = y(N:2*N-1);%后N个输入数据
win = hanning(N)';;
win1 = win/sum(win);%窗归1
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);%FFT振幅谱
p1 = mod(phase(y11_fft)*180/pi,360);%FFT相位谱
y2 = y(1:2*N-1);%2N-1个输入数据
%win = hanning(N)';;
winn = conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);%窗归1
y22= y2.*win2;
y222=y22(N:end)+[0 y22(1:N-1)];%构成长N的apFFT输入数据
y2_fft = fft(y222,N);;
a2 = abs(y2_fft);%apFFT振幅谱
p2= mod(phase(y2_fft)*180/pi,360);%apFFT相位谱
ee=mod((p1-p2)/180/(1-1/N),1);%频率偏离校正值
aa=(a1.^2)./a2*2;%振幅校正值
a11=a2./a1;
a12=a11.^2;
for i=1:3
aa2=(A(i)/N).*(sin(pi*ee)/sin(pi*ee/N)).*exp(pha(i)+N/(N-1)*pi*ee);
aa3=(A(i)/N^2).*(sin(pi*ee)/sin(pi*ee/N)).^2.*exp(pha(i));
end
subplot(5,2,1);stem(2*a1,'.');title('（a）FFT振幅谱');ylim([0,1]);xlim([0 N/2]);grid
subplot(5,2,2),stem(2*a2,'.');title('（b）apFFT振幅谱');ylim([0,1]);xlim([0 N/2]);grid
subplot(5,2,3);stem(p1,'.');title('（c）FFT 相位谱');ylim([0,400]);xlim([0 N/2]);grid
subplot(5,2,4),stem(p2,'.');title('（d）apFFT 相位谱');ylim([0,400]);xlim([0 N/2]);grid
subplot(5,2,5);stem(ee,'.');title('（e）校正频偏');ylim([0,1]);xlim([0 N/2]);grid
subplot(5,2,6);stem(aa,'.');title('（f）校正振幅1');ylim([0,1]);xlim([0 N/2]);grid
subplot(5,2,7);stem(a11,'.');title('（g）归一振幅谱');ylim([0,1]);xlim([0 N/2]);grid
subplot(5,2,8);stem(aa2,'.');title('（h）校正振幅2');ylim([0,1]);xlim([0 N/2]);grid
subplot(5,2,9);stem(a12,'.');title('（i）归一功率谱');ylim([0,1]);xlim([0 N/2]);grid
subplot(5,2,10);stem(aa3,'.');title('（j）校正振幅3');ylim([0,1]);xlim([0 N/2]);grid
for i=1:3
r(i)=round(f(i));
disp('频率校正值')
fff(i)=floor(f(i))+ee(r(i)+1)
disp('振幅校正值')
aaa(i)=aa(r(i)+1)
disp('初相位校正值')
ppp(i)=p2(r(i)+1)
end
程序改动的地方编的有点烂，还望老师不吝赐教，谢谢！ |
发表于 2007-5-29 08:03
显身卡
