zoomfft的matlab程序(1)

%ZoomFFT谱的matlab程序
%x-信号序列
%fs-采样频率
%N-做谱点数
%fe-分析中心频率
%D-细化倍数
%L-平均段数
%M-滤波器半阶数
%f-返回频率向量
%xz-返回幅值谱

function [f xz]=ZoomFFT(x,fs,N,fe,D,L,M)
k=1:M;
w=0.5+0.5*cos(pi*k/M); %Hanning窗

fl=max(fe-fs/(4*D),-fs/2.2);%细化频率上线
fh=min(fe+fs/(4*D),fs/2.2);%细化频率下线

yf=D*fl; %移频量
df=fs/D/N;
f=fl:df:fl+(N/2-1)*df;
xz=zeros(1,N/2);

wl=2*pi*fl/fs;
wh=2*pi*fh/fs;
hr(1)=(wl-wh)/pi;
hr(2:M+1)=(sin(wl*k)-sin(wh*k))./(pi*k).*w;
hi(1)=0;
hi(2:M+1)=(cos(wl*k)-cos(wh*k))./(pi*k).*w;
k=0:N-1;
w=0.5-0.5*cos(2*pi*k/N);
for i=1:L
for k=1:N
kk=(k-1)*D+M+(i-1)*N;
xrz(k)=x(kk+1)*hr(1)+sum(hr(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
xiz(k)=x(kk+1)*hi(1)+sum(hi(2:M+1).*(x(kk+2:kk+M+1)-x(kk:-1:kk-M+1)));
end
xzt=(xrz+j*xiz).*exp(-j*2*pi*(0:N-1)*yf/fs);
xzt=xzt.*w;
xzt=xzt-sum(xzt)/N;
xzt=fft(xzt);
xz=xz+(abs(xzt(1:N/2))/N*2).^2;
end
xz=(xz/L).^0.5;


例子：
fs=10240;
N=1024;
D=100;
M=0;
x=cos(2*pi*256.4*t)+5*cos(2*pi*256.9*t);
[f xz]=ZoomFFT(x,fs,N,256,D,1,M);
plot(f,xz);



所谓细化，我认为定义应该是：（不增加信号时域长度的情况下），提高其频谱分析物理分辨率的操作。
ZoomFFT，实现时，如果细化为原来的D倍,变为Delta_f/D；则要求采样长度相应的也应该是原来的D倍，D*N；而既然如此，直接经FFT得到的频谱，其实分辨率本身就是Delta_f/D。 ZoomFFT所做的， 只不过直接FFT计算点数是D*N，ZoomFFT通过巧妙的办法，使计算量仍旧是N；并且，只选取我们所关心的“频带”去计算。 对我们而言，我们也可以直接FFT得到频谱后，只绘制我们关心的那个“频带”；相对该直接FFT方法，在获得信息上ZoomFFT没有任何增加；仅仅提高了效率（连这其实都值得商榷，如果我们同时有关心其他频率，ZoomFFT还得重算）。
所以，我认为ZoomFFT的更确切的称呼，应该是“选带分析”，而非“细化分析”；当然这是沿用我理解的前面的细化的定义。