zoomfft.m

FFT变换的一个东西

%ZoomFFT谱
%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;