频谱细化之ZOOMFFT算法

---

一、 算法简介

复调制 Zoom-FFT[64]的核心思想是对有限长连续的中频信号进行 FFT 变换得 到粗略的信号序列感兴趣区间，然后对粗略的信号序列执行第二次 FFT 分析， 进而得到精确的频谱。
若已知离散的中频信号序列 x(m),采样点数为 M，采样频率为 fs，放大倍数为 D，并做 M 点 FFT 变换，得到输出序列 x(k)。
算法的执行步骤如下：
(1)频移(复调制)
复调制是将感兴趣的频谱的起始点移动到零频点，确保感兴趣信号扩展到频 谱的显示范围内。若已知频谱的起止频率为 f1~f2，则信号的中频频率 f0可表示为：


式中， fs = M *f ,f 为谱线间隔， f0 / f 表示 f0的中心移位，即整个频谱中心 频率 f0所对应的频谱的序列。可以得知，复调制将频率分量 f0平移到 x1(m)的零频处。而为了将 X(k)零点附近的频谱细化，可以把频率降到 fs / D 重新抽样,D 是细化倍数。
(2)数字低通滤波
低通滤波可以保证频谱变换后的频率不会因频谱混叠而受到干扰，若选择的 细化倍数为 D，那么数字低通滤波器的截止频率要满足 fc  fs / 2D 。
(3)重采样
降低采样频率至 fs/D,然后重采样可以提高频率分辨率。根据采样定理，当采 样点保持 N 不变时，若采样频率降低至原来的 1/D，相当于采样窗口时间增长 D 倍，这时候分辨率变为原来的 1/D。
(4)复数 FFT
重采样后，信号的实部和虚部被分开，所以必须对信号做 N 点复 FFT，可以 得到 N 条谱线，便得到了细化的频谱。如果细化倍数为 D，重采样频率将下降 为原来的 1/ D，采样时间也增大到原来的 D 倍。需要注意的是，采样时间的增 大仅在第一次采样时增加了采样点数，而在后续重采样时，仍为 M 个采样点， 并没有增加之后的 FFT 运算时间。

二、算法仿真

clear;
close all;
clc;
fs=1024;   %采样频率
N=1024;    %傅里叶变换点数
D=50;      %细化倍数
M=200;     %滤波器阶数
t=(0:N*D+2*M)/fs;    %时间轴

% x = 1.5 * cos(2 * pi * 98 * t);

x = 1.5 * cos(2 * pi * 98 * t) + 2 * cos(2 * pi * 99 * t) + ...
3 * cos(2 * pi * 100. * t) + 3.5 * cos(2 * pi * 101 *t);


xf=fft(x,N);     %傅里叶变换
xf=abs(xf(1:N/2))/N*2;

fe=99.5;     %中心频率

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;

figure();
subplot(211);
plot((0:N/2-1)*fs/N, 2 * abs(xf(1:(N/2))) / N, 'b', 'LineWidth', 3);
grid on;
title('FFT和ZoomFFT仿真对比','FontSize',30);
ax1 = gca;  % 获取当前子图的坐标轴对象
ax1.FontSize = 25;  % 设置坐标轴上文字的字体大小
xlabel('Frequency (Hz)','FontSize',30);
ylabel('Amplitude','FontSize',30);
legend('FFT','FontSize',25);

subplot(212);
plot(f, xz, 'r', 'LineWidth', 3);
grid on;
ax1 = gca;  % 获取当前子图的坐标轴对象
ax1.FontSize = 25;  % 设置坐标轴上文字的字体大小
xlabel('Frequency (Hz)','FontSize',30);
ylabel('Amplitude','FontSize',30);
legend('ZoomFFT','FontSize',25);

三、仿真结果

---