频谱细化-----Zoom-FFT算法介绍及MATLAB实现

栅栏效应

若长度为$N=T/dt=T\cdot f_s$，其中$T$为采样时间，$dt$为采样间隔，$f_s$为采样频率，通过傅里叶变换后得到$X\left(f_i\right)$，$f_i=i\cdot f_s/N,i=0,1,2,...,N/2$，由于频率分辨率的限制，那么相邻两个离散点之间的频率无法获得，会导致部分有用的频率被漏掉，这就是所谓的栅栏效应。
解决这种问题的一种方法是可以增加数据的点数或者提高傅里叶变换的点数，但是这样会增加系统的硬件成本。不过可以通过某些方法在不改变硬件的基础上提高频率的分辨率，下面简单介绍一种经典的方法Zoom-FFT算法。

原理介绍

在实际的应用中，只需要对感兴趣的某段频率的信号进行分析，下面介绍一种基于复调制的频谱细化方法。

复调制移频

模拟信号$x\left(t\right)$经过A/D变换后转化为离散信号$x_0\left(n\right)$。如果感兴趣的频带为$f_1\sim f_2$，中心频率为$f_e=\left(f_1+f_2\right)/2$，对$x_0\left(n\right)$进行复调制，可以得到
$$\begin{array}{rl}& x\left(n\right)=x_0\left(n\right)e^{-j2\pi n{f}_e/f_s}=x_0\left(n\right)\cos \left(2\pi n{f}_e/f_s\right)-j{x}_0\left(n\right)\sin \left(2\pi n{f}_e/f_s\right)\\ & =x_0\left(n\right)\cos \left(2\pi n{L}_0\Delta f/N\Delta f\right)-j{x}_0\left(n\right)\sin \left(2\pi n{L}_0\Delta f/N\Delta f\right)\\ & =x_0\left(n\right)\cos \left(2\pi n{L}_0/N\right)-j{x}_0\left(n\right)\sin \left(2\pi n{L}_0/N\right)\end{array}$$
其中，$f_s=N\Delta f$是采样频率，$\Delta f$是谱线间隔，$L_0=f_e/\Delta f$为频率的中心移位，频谱显示对应中心频率的谱序号，即$f_e=L_0\Delta f$。同时可以看出，复调制使得$x_0\left(n\right)$的$f_e$搬移到$x\left(n\right)$的零频点，即$X_0\left(k\right)$的第$L_0$条谱线移到$X\left(k\right)$中零点频谱的位置.为了得到$X\left(k\right)$零点附近的部分细化频谱，可重新抽样把频率降到$f_s/D$，$D$为细化倍数。为了使得抽样后的频率不发生频谱混叠，需要在抽样前进行低通滤波。

低通滤波

为了保证重新采样后的信号在频谱分析时不发生频谱混叠，需进行抗混叠滤波，滤出需要分析的频段信号，则数字低通滤波器的截止频率 $f_c\le f_s/2D$

重新抽样

信号经过移频、低通滤波后，分析信号点数变少，再以较低的采样频率进行重新采样，在通过补零保证相同的采样点数时，样本的总长度加大，频谱的分辨率也就得到了提高。
设原采样频率为 ,采样点数为$N$，则频率分辨率为$f_s/N$，现重采样频率为$f_s/D$，当采样点数仍是$N$时，其分辨率为$f_s/\left(D\ast N\right)$，分辨率提高了$D$倍。这样就在原采样频率不变的情况下得到了更高的频率分辨率。

复数FFT

重新采样后的信号实部和虚部是分开的，需要对信号进行N点复FFT，从而得出$N$条谱线，此时分辨率为$\Delta f^{{}^{\prime }}=f{{}_s}^{{}^{\prime }}/N=f_s/ND=\Delta f/D$，可见分辨率提高了$D$倍。

频率调整

经过算法运行后的谱线不为实际频率的谱线，需要将其反向搬移，转换成实际频率，进而得出细化后的频率。
下面对该方法进行仿真，仿真参数如下

参数值	参数名称
采样频率	1000Hz
细化倍数	50
滤波器阶数	200
FFT点数	2048
频率	166.4Hz、165Hz

结果如下

从图中可以看出，直接进行FFT，频率与原始频率有些许误差，并且不是很能分辨，但是细化后的频谱能够很好地区分信号的频率，很准确地得到了信号的频率。
代码如下：

clear;
close all;
clc;
fs=1000;   %采样频率
N=2048;    %傅里叶变换点数
D=50;      %细化倍数
M=200;     %滤波器阶数

t=(0:N*D+2*M)/fs;    %时间轴
x=4*sin(2*pi*166.4*t)+2*sin(2*pi*165*t+pi/6)...
    +0.6*randn(1,N*D+2*M+1);  %信号

xf=fft(x,N);     %傅里叶变换
xf=abs(xf(1:N/2))/N*2;
subplot(211);
plot((0:N/2-1)*fs/N,xf);
axis([163.5 169 0 4])
grid on
xlabel('f/Hz');ylabel('Amplitude')
title('直接进行傅里叶变换结果')
fe=167;     %中心频率

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;

p=0:N-1;
w=0.5-0.5*cos(2*pi*p/N);
xrz=zeros(1,N/2);
xiz=zeros(1,N/2);
L=10;   %循环次数
for i=1:L
    for k=1:N
        kk=(k-1)*D+M;
        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+1j*xiz).*exp(-1j*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;
subplot(212);
plot(f,xz);
axis([163.5 169 0 4])
grid on
xlabel('f(Hz)');ylabel('Amplitude')
title('Zoom-FFT细化后的频谱')

参考文献：
[1]王旭刚. 基于FMCW体制K波段测距雷达的研究与实现[D].南京航空航天大学,2012.