选带快速傅立叶变换ZOOM-FFT的matlab实现

up目录

一、理论基础

二、核心程序

三、测试结果

---

一、理论基础

      zoom-fft是一个信号传输过程，其中输入信号被向下混频到基带，然后被抽取，然后被传递到标准FFT。ZOOM-FFT称为细化的快速傅立叶变换，又称为选带快速傅立叶变换。ZOOM-FFT的功能是对信号的频率进行局部细化放大，使感兴趣的频带获得较高的频率分辨率。实现FFT细化功能的算法有几种，如频移法-、相位补偿法和最大熵谱法等，目前应用最关的是频移法。频移细化FFT的方法是基于离散傅立叶变幻的频移原理。设fk为所需要细化的频带的中心频率，对信号乘以exp(-j2pifkt)进行数字化平移后，原fk处的谱线移至频率轴的0处，用低通滤波器滤除所需细化的频段外的频率成分，并以细化倍数为间隔进行重采样，最后对重采样后的数据做FFT变换，并对变换结果重新排序。

       zoomFFT这种频率细化的方法主要步骤是：移频-抗混叠滤波-重采样-复FFT处理-频率调整，这种方法是利用降低采样频率Fs（但又不会产生频率混叠），而采样点数N不变的，来提高所关心频段的频率分辨率（单纯的FFT在降低Fs时受采样定理限制，不能太小，不然会产生混叠）。其实这种方法方法，在降低采样频率，保持采样点数不变的情况下，必然还是会导致采样时间的增加，这是我们不希望看到的。而减少计算量是对的，zoomFFT可以选择感兴趣的频段对其进行频谱分析（要经过滤波处理），相对来说所需的采样频率Fs1就比较低，而全频段的FFT分析，为了防止频率混叠，需要满足Fs2>2*Fc （Fc为信号中最大频率），可见Fs1<Fs2。而频谱分辨率为Fs/N = Fs/(Fs*t)=1/t 只与采样时间有关（t为采样时间），这时要使这两种算法频率分辨率相同，只需要使其采样时间相同，由于Fs1<Fs2，则ZoomFFT的采样点数N1<全频段FFT的采样点数N2，从而zoomFFT的计算量是N1*log2(N1)，FFT的计算量是N2*log2(N2)，相对来说计算量减小了。

细化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变换结果。

二、核心程序

主要步骤如下：

% 步骤一：乘以exp
zoom_fft_xx   = (x_real_zoom+j*x_imag_zoom).*exp(-j*2*pi*(0:N-1)*frequency_shift/Fs);
% 步骤二：数字低通+重采样
zoom_fft_xx   =  zoom_fft_xx.*w;
zoom_fft_xx   =  zoom_fft_xx-sum(zoom_fft_xx)/N;
% 步骤三：FFT变化
zoom_fft_xx   =  fft(zoom_fft_xx);
% 步骤四：频率调整
zoomfft_x     =  zoomfft_x+(abs(zoom_fft_xx(1:N/2))/N*2).^2;
zoomfft_x     =  sqrt(zoomfft_x);

%加载波形信号，用于处理
load data_sin_cos_square.mat

D     = str2num(get(handles.ZoomFactor,'string'));%放大倍速，最大为10，大于10.代码报错
f1    = str2num(get(handles.ZoomF1    ,'string'));%分辨频率下限
f2    = str2num(get(handles.ZoomF2    ,'string'));%分辨频率上限
ws    = str2num(get(handles.edit18    ,'string'));%通带
wp    = str2num(get(handles.edit19    ,'string'));%阻带
%进行原始的傅立叶变化
fft_x = fft(x,N);
fft_x = abs(fft_x(1:N/2))/N*2;
k     = 1:M;                          
w     = 0.5+0.5*cos(pi*k/M);          %Hanning窗

%移频量
frequency_shift = D*f1;                       %移频量
frequency_diff  = Fs/D/N;
f               = f1:frequency_diff:f1+(N/2-1)*frequency_diff;
zoomfft_x       = zeros(1,N/2);
%定义滤波器
%为了便于观察，滤波器的分实部和虚部
H_filter_real(1)     = (ws-wp)/pi;
H_filter_real(2:M+1) = (sin(ws*k)-sin(wp*k))./(pi*k).*w;
H_filter_imag(1)     =  0;
H_filter_imag(2:M+1) = (cos(ws*k)-cos(wp*k))./(pi*k).*w;

k=0:N-1;
w=0.5-0.5*cos(2*pi*k/N);%Hanning窗

for k=1:N
    kk=(k-1)*D+M;
    x_real_zoom(k)=x(kk+1)*H_filter_real(1)+sum(H_filter_real(2:M+1).*(x(kk+2:kk+M+1)+x(kk:-1:kk-M+1)));
    x_imag_zoom(k)=x(kk+1)*H_filter_imag(1)+sum(H_filter_imag(2:M+1).*(x(kk+2:kk+M+1)-x(kk:-1:kk-M+1)));
end

%这里开始正式处理
% 步骤一：乘以exp
zoom_fft_xx   = (x_real_zoom+j*x_imag_zoom).*exp(-j*2*pi*(0:N-1)*frequency_shift/Fs);
% 步骤二：数字低通+重采样
zoom_fft_xx   =  zoom_fft_xx.*w;
zoom_fft_xx   =  zoom_fft_xx-sum(zoom_fft_xx)/N;
% 步骤三：FFT变化
zoom_fft_xx   =  fft(zoom_fft_xx);
% 步骤四：频率调整
zoomfft_x     =  zoomfft_x+(abs(zoom_fft_xx(1:N/2))/N*2).^2;
zoomfft_x     =  sqrt(zoomfft_x);
%定位到第二个图上
axes(handles.axes2);
plot((0:N/2-1)*Fs/N,fft_x);hold on;%显示原始的频谱图
plot(f,zoomfft_x,'r');hold off;    %显示zoomfft处理后的频谱图
axis([f1 f2 0 1.4*max(fft_x)]);    %显示坐标范围
grid on;                           %显示方格
up24

三、测试结果