ZOOM-FFT详解

FFT虽然迅速方便地实现了数字频谱分析的功能，但是对于希望在某些感兴趣的频率范围内以较高的分辨率进行频域分析，传统的 FFT 就做不到了。因为普通 FFT在频域的分辨率是有一定限度的，其分析频率范围总是从0Hz 开始到 1/2倍的采样频率为最高分析频率。这种分析方法称为基带变换分析法。设采样频率为fs，FFT点数为NFFT，频率分辨率△f的关系式为△f =fs/NFFT。可以看出，一旦采样频率fs和NFFT确定后，频率分辨率是一个固定值。

ZOOM-FFT称为细化的快速傅立叶变换，又称为选带快速傅立叶变换。ZOOM-FFT的功能是对信号的频率进行局部细化放大，使感兴趣的频带获得较高的频率分辨率。实现FFT细化功能的算法有几种，如频移法-、相位补偿法和最大熵谱法等，目前应用最关的是频移法。

频移细化FFT的方法是基于离散傅立叶变幻的频移原理。设fk为所需要细化的频带的中心频率，对信号乘以exp(-j2pifkt)进行数字化平移后，原fk处的谱线移至频率轴的0处，用低通滤波器滤除所需细化的频段外的频率成分，并以细化倍数为间隔进行重采样，最后对重采样后的数据做FFT变换，并对变换结果重新排序。细化FFT的具体步骤如下：

1，遵循采样定理的原则，为防止采样信号的频率混淆，首先需要通过模拟低通抗混叠滤波器滤波或这顶足够高的采样频率fs。然后需要采集足够长度的信号数据，数据的长度为细化倍数D与NFFT的乘积，即为DNFFT。

2，将采样信号进行频移（复调制），即乘以单位旋转因子exp(-j2pifkt)，这样就把频率原点由0处移到所需要细化的频率fk处，频率分量fk停留在频率为零出的位置上，形成了一个以fk为频率零点的新的信号xk(t)。

3，用低通数字滤波器对频移后的数据进行滤波，去除信号所需要细化频带外的频率成分。

4，对滤波后数据进行重采样，重采样的采样频率为fd/D，也就是每隔细化倍数（D-1）个点取一个数据。

5，对重采样后的数据进行长度为NFFT的FFT计算。

6，对FFT计算结果重新排序。对于重采样所得的是一个复值序列，在进行FFT计算时，全部数据都是有用的信息。因为它是以调制频率fk为新的零频率，实际上不存在对称性，所以负频率处的一般复数结果都是有效的。在MATLAB函数FFT的变化结果中，由于对应正频率的数据是从头开始按正向排列存放在结果向量的前半部分，对饮复频率的数据则是从尾开始按反向排列存放在数据的后半部分。如果过用户只需要分析细化FFT的正频率域的，需要将细化FFT结果排列在（3NFFT/4+1）–NFFT的数据赋值到一个新向量1–NFFT/4的单元，将排列在1–NFFT/4的数据赋值到这个新向量的（NFFT/4+1）–NFFT/2的单元，于是便可的到频率从0到折叠频率的正频率域的细化FFT变换结果。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 细化FFT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear               % Remove items from workspace, freeing up system memory

clc                 % Clear Command Window

clf                 % Clear current figure window

close all hidden   % removal draws only those lines that are not obscured

                    % by other objects in the field of view. 

format long

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 加载数据

fun = @(t) sin(5*2*pi*t) +sin(20*2*pi*t) + randn(size(t));

fs = 300;                       % 采样频率

n = 10000;                       % 数据长度

t = 0:1/fs:(n-1)/fs;             % 建立离散时间列向量

x = fun(t);                     % 产生时间序列数据

fi = 15;                         % 最小细化截止频率

np = 8;                         % 放大倍数

nfft = 1024;                     % FFT长度

nt = length(x);                 % 计算输入数据长度 

fa = fi+0.5*fs/np;               % 最大细化截止频率  

nf = 2^nextpow2(nt);             % 取大于nt且最接近的2的整数次方为FFT长度

na = round(0.5*nf/np+1);         % 确定细化带宽的数据长度 

% 频移

n = 0:nt-1;                     % 建立一个按1递增的向量

b = n*pi*(fi+fa)/fs;             % 确定单位旋转因子向量 

y = x.*exp(-1i*b);               % 乘以单位旋转因子进行频移

% 频移后的低通滤波（频域滤波）

b = fft(y,nf);                   % FFT变换

a(1:na) = b(1:na);               % 正频率带通内的元素赋值

a(nf-na+1:nf) = b(nf-na+1:nf);   % 负频率带通内的元素赋值

b = ifft(a,nf);                 % IFFT变换

c = b(1:np:nt);                 % 重采样 

a = fft(c,nfft)*2/nfft;         % 进行细化FFT变换  

y2 = zeros(1,nfft/2);           % 变换结果重新排序

y2(1:nfft/4) = a(nfft-nfft/4+1:nfft); % 排列负频率段的数据

y2(nfft/4+1:nfft/2) = a(1:nfft/4);     % 排列正频率段的数据 

n = 0:(nfft/2-1);

f2 = fi + n*2*(fa-fi)/nfft;           % 定义细化FFT频率向量    

% 对输入数据做FFT用来变换比较  

y1 = fft(x,nfft)*2/nfft;               % FFT变换

f1 = n*fs/nfft; 

ni = round(fi*nfft/fs+1);             % 定义与细化FFT频率向量相同的频率范围

na = round(fa*nfft/fs+1);

% 绘制输入时程曲线图形

subplot(2,1,1); plot(t,x);

xlabel('时间(s)');ylabel('幅值');grid on;

title('输入时程曲线图形');axis([min(t) max(t) min(x) max(x)])

% 在相同的频率范围下绘制幅频曲线图

nn = ni:na;

subplot(212);plot(f1(nn),abs(y1(nn)),':',f2,abs(y2));

xlabel('频率(Hz)');ylabel('幅值');legend('普通','细化');grid on;    

title('相同的频率范围下绘制幅频曲线图');