论文阅读1：分段FFT算法实现

1.导言

        海量数据的信号频率估计是工程实践中常见的难题，其庞大的运算代价给数据处理带来了困难。 而常规的频率估计算法为了保证估计精度，通常会进一步增加运算代价。对此，提出了一种用于过采样场景下海量数据的基于频谱聚焦的频率估计方法。该方法将数据分段进行频谱聚焦处理，降低了运算代价， 解决了海量数据处理困难的问题，同时提高了频率估计精度。仿真结果表明，与传统算法相比，该算法具有运算代价较小，估计精度较高的优势。

2.论文简介

论文题目： 基于频谱聚焦的海量数据信号频率估计方法

关键词：频率估计；过采样；大数据量；运算代价；高精度；频谱聚焦

论文作者： 蔡博文 ，马武舰   ，陈新亮 ，梁振楠   ，刘泉华 

论文链接：

https://kns.cnki.net/kcms2/article/abstract?v=VcTOyLYtvEw5IXsm6moTqurIkrFfsxrPgjYAmPbKdZ5wLa30CILW8gaYa0vlUGrIc6uQrorFXmyt-5JC89qdIgCvRVdfZtl8RfVVk-kQabMF_S9PoP4GQmibCD1aYQMGuzRfqjFNBWOPWr32URpUbahWd7feNJ_RcFzo_8UhyDpMSG6sZpFJvQf4SmC5nbfD5tP7KuNeu14=&uniplatform=NZKPT&language=CHS

3.频谱聚焦算法原理

        为了缩短单次 FFT 的数据长度，并提高频谱 精度，提出频谱聚焦的处理方法。 考虑实际接收的数据流，每次截取一个 PRT 进行频谱聚焦处理，处理流程如图 1 所示。

        处理步骤包括：

        （1 ）对接收数据流截取一个 PRT ，然后对截取 的 PRT 间隔 M 位抽值 ，分为 M 个子通道。

        （2 ）各子通道做 L倍补零傅里叶变换 ，得到其 频谱。

        （3 ）各子通道频谱， 补偿相位差 后进行 叠加 ，得到叠加后的频谱，即为频谱聚焦结果。

4.频谱聚焦算法MATLAB复现

% 参数设置
fs = 1024; % 采样频率
N = 16384; % 信号长度
f1 = 50; % 第一个正弦波频率
f2 = 223; % 第二个正弦波频率
f3 = 483; % 第三个正弦波频率
f4 = 51; % 第四个正弦波频率
f5 = 52; % 第五个正弦波频率
A1 = 1; % 第一个正弦波幅度
A2 = 0.7; % 第二个正弦波幅度
A3 = 0.5; % 第三个正弦波幅度
A4 = 0.9; % 第二个正弦波幅度
A5 = 0.8; % 第三个正弦波幅度

% 生成复合信号
t = (1:N)/fs;
x0 = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t) + A3*sin(2*pi*f3*t)+A4*sin(2*pi*f4*t)+A5*sin(2*pi*f5*t);

% 直接进行FFT
X_direct = fft(x0, N);
f_direct = (0:N-1)*(fs/N);

% 绘制直接FFT结果
figure;
subplot(2, 1, 1);
plot(f_direct, abs(X_direct));
title('Direct FFT Result');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;

% 频谱聚焦方法参数
M = 8; % 子通道数
L = 1; % 补零倍数
PRT = N; % 假设PRT等于整个信号长度

% 数据分段和抽样
x = zeros(M, N/M);
for q = 1:M
    x(q, :) = x0((q-1)*N/M + 1:q*N/M);
end

% 傅里叶变换
Xq = zeros(M, L*N/M);
for q = 1:M
    Xq(q, :) = fft(x(q, :), L*N/M);
end

% 相位补偿和叠加（这里为了简化，没有进行实际的相位补偿计算）
X = zeros(1, L*PRT/M);
for q = 1:M
    % 相位补偿，这里假设信号频率为50Hz，实际情况需要根据信号特性来确定
    phase_compensation = exp(-1j*2*pi*(0:L*PRT/M-1)*q/fs/128);
    X = X + Xq(q, :) .* phase_compensation;
end

% 频率轴
f_focus = (0:L*N/M-1)*(fs/L/N);

% 绘制频谱聚焦结果
subplot(2, 1, 2);
plot(f_focus*8, abs(X));
title('Spectrum Focusing Result');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;

% 输出峰值频率
[~, idx] = max(abs(X));
fprintf('Peak frequencies: ');
for i = 1:length(idx)
    fprintf('%f Hz ', f_focus(idx(i)));
end
fprintf('\n');

        代码中比较了16384长度数据16384点直接FFT和分成8个子通道即8段2048点FFT后聚焦之后的结果对比，代码中绘制频谱聚焦结果时为了方便与直接FFT比较，乘以M即通道数，最后显示直接FFT和分段FFT之后聚焦乘以M的频率估计结果基本相同。其结果如下图所示：

              

5.实验结果及分析：最大分辨率试验

% 参数设置

fs = 1024; % 采样频率
N = 2048; % 信号长度
f1 = 50; % 第一个正弦波频率
f2 = 224; % 第二个正弦波频率
f3 = 484; % 第三个正弦波频率
f4 = 51; % 第四个正弦波频率
f5 = 52; % 第五个正弦波频率
A1 = 1; % 第一个正弦波幅度
A2 = 0.7; % 第二个正弦波幅度
A3 = 0.5; % 第三个正弦波幅度
A4 = 0.9; % 第二个正弦波幅度
A5 = 0.8; % 第三个正弦波幅度
M = 8; % 子通道数
L = 1; % 补零倍数

                               

如图所示为上面参数的直接fft和分段fft对比图，共五组频率，分别为：

              f1 = 50; % 第一个正弦波频率

              f2 = 224; % 第二个正弦波频率

              f3 = 484; % 第三个正弦波频率

              f4 = 51; % 第四个正弦波频率

              f5 = 52; % 第五个正弦波频率

       f1，f4，f5为相近频率，可以看出算法无法分辨。

% 参数设置

fs = 1024; % 采样频率
N = 2048; % 信号长度
f1 = 50; % 第一个正弦波频率
f2 = 223; % 第二个正弦波频率
f3 = 481; % 第三个正弦波频率
f4 = 51; % 第四个正弦波频率
f5 = 52; % 第五个正弦波频率
A1 = 1; % 第一个正弦波幅度
A2 = 0.7; % 第二个正弦波幅度
A3 = 0.5; % 第三个正弦波幅度
A4 = 0.9; % 第二个正弦波幅度
A5 = 0.8; % 第三个正弦波幅度
M = 8; % 子通道数
L = 1; % 补零倍数

          

如图所示为上面参数的直接fft和分段fft对比图，共五组频率，分别为：

         f1 = 50; % 第一个正弦波频率

         f2 = 223; % 第二个正弦波频率

         f3 = 481; % 第三个正弦波频率

         f4 = 51; % 第四个正弦波频率

         f5 = 52; % 第五个正弦波频率

f1，f4，f5为相近频率，可以看出算法无法分辨。f2,f3不是M/2的倍数无法分辨。

% 参数设置

fs = 1024; % 采样频率
N = 2048; % 信号长度
f1 = 48; % 第一个正弦波频率
f2 = 224; % 第二个正弦波频率
f3 = 484; % 第三个正弦波频率
f4 = 56; % 第四个正弦波频率
f5 = 64; % 第五个正弦波频率
A1 = 1; % 第一个正弦波幅度
A2 = 0.7; % 第二个正弦波幅度
A3 = 0.5; % 第三个正弦波幅度
A4 = 0.9; % 第二个正弦波幅度
A5 = 0.8; % 第三个正弦波幅度
M = 8; % 子通道数
L = 1; % 补零倍数

                          

  

如图所示为上面参数的直接fft和分段fft对比图，共五组频率，分别为：

        f1 = 48; % 第一个正弦波频率

        f2 = 224; % 第二个正弦波频率

        f3 = 484; % 第三个正弦波频率

        f4 = 56; % 第四个正弦波频率

        f5 = 64; % 第五个正弦波频率

f1，f4，f5为相近频率，可以看出算法相近频率必须为M的倍数。f2,f3为不相近频率，不是M/2的倍数则算法无法分辨。

% 参数设置

fs = 1024; % 采样频率
N = 16384; % 信号长度
f1 = 50; % 第一个正弦波频率
f2 = 223; % 第二个正弦波频率
f3 = 483; % 第三个正弦波频率
f4 = 51; % 第四个正弦波频率
f5 = 52; % 第五个正弦波频率
A1 = 1; % 第一个正弦波幅度
A2 = 0.7; % 第二个正弦波幅度
A3 = 0.5; % 第三个正弦波幅度
A4 = 0.9; % 第二个正弦波幅度
A5 = 0.8; % 第三个正弦波幅度

   

如图所示为上面参数的直接fft和分段fft对比图，共五组频率，分别为：

          f1 = 50; % 第一个正弦波频率

          f2 = 223; % 第二个正弦波频率

          f3 = 483; % 第三个正弦波频率

          f4 = 51; % 第四个正弦波频率

          f5 = 52; % 第五个正弦波频率

采样点数扩大M倍之后，频率分辨率达到1Hz.

    结论：要达到与直接fft相同的分辨率，分段fft算法需采样点数扩大到原来的M（子通道数）倍。

6.总结

     在过采样场景下，本文所提出的基于频谱聚焦 的信号频率估计算法，减小了单次 FFT 点数，解决了海量数据处理困难的问题。同时，其估计精度与补零 FFT 相近，但要达到与直接fft相同的分辨率，分段fft算法需采样点数扩大到原来的M（子通道数）倍。该算法具有一定的工程实践意义，在射电天文的雷达校正等领域中有良好的应用前景。