基于跨帧相位差法实现混合信号分离
前言
本文的背景，是设计一个信号分离装置，用来处理由两个正弦波叠加而成的混合信号。更准确地说，这里所说的“分离”，本质上是要估计两个原始信号各自的频率、幅度和相位，并尽可能提高测量精度。原始信号的频率范围为 $100Hz \sim 10kHz$，且两信号的最小频差为 $100Hz$。
这个问题最先想到的工具通常是快速傅里叶变换（Fast Fourier Transform, FFT）。但如果直接依赖 FFT 的频率栅格来读数，精度往往会受到频率分辨率 $\Delta f = f_s/N$ 的限制。举个简单的例子：若采样率取 $25kHz$，想让频率分辨率达到 $1Hz$，就至少需要 $N=25000$ 个采样点；而为了适配嵌入式上常见的 FFT 实现，点数通常还会进一步取成 $2$ 的整数次幂，也就是说往往要取到 $32768$ 点。仅仅存储一帧 uint16_t 采样值，就已经要占用约 $64KB$ RAM；若再算上 FFT 输入输出缓冲区、复数数组以及其他中间变量，对 STM32 这类单片机来说压力就很大了。
因此，这篇文章的核心思路并不是“单纯把 FFT 做得更大”，而是先用 FFT 找到信号所在的粗略频率位置，再利用跨帧相位差法（phase difference across frames）对频率进行精细修正。当频率估计足够准确后，再使用最小二乘法（Least Squares, LS）去估计幅度与相位，整体精度就能明显提高。
需要先说明一点：这里的“跨帧相位差法”并不是跳过 FFT，而是建立在“已经知道信号大致落在哪个频率 bin 附近”的前提之上的。也就是说，在进入相位差测频之前，我们仍然需要先完成频谱粗定位。
在 FFT 之前
由傅里叶变换衍生而来的快速傅里叶变换（Fast Fourier Transform, FFT），能够把时域信息转换为频域信息，因此是分析信号最常用的工具之一。理想的傅里叶变换默认信号在整个时间轴上都存在，而实际采样只能截取一小段有限长度的数据。也正因为如此，直接对有限长度的采样序列做离散傅里叶变换（Discrete Fourier Transform, DFT）时，首先会遇到的问题就是频谱泄露。
频谱泄露
更专业地说，当我们对一个时域信号做离散傅里叶变换（Discrete Fourier Transform, DFT）时，如果该信号在当前观测窗口内不是整数周期，或者说信号周期与 DFT 窗口长度不匹配，那么信号能量就不会只集中在某一个频率 bin（频率“桶”）上，而会扩散到相邻 bin 中，形成旁瓣。这种“能量扩散”的现象，就是频谱泄露。
可以把它理解成这样：我们是拿一个有限长度的“窗户”去观察一段原本可以无限延伸的波形。如果窗内正好包含了整数个周期，那么截断相对干净；但如果只截到了半个周期或者若干个不完整周期，那么在窗边界处就会出现突变。这个突变并不是信号本身的频率成分，而是截断操作带来的额外频谱扩展。
从数学角度看，采样本质上相当于将原始信号乘上一个有限长度窗函数。若给定原始连续信号 $y(t)$，采样时间区间为 $[0,T]$，那么采样得到的离散信号可以看成是 $y(t)$ 乘以一个矩形窗函数 $w(t)$：
$$ w(t) = \begin{cases} 1, \quad 0 \leq t \leq T \\ 0, \quad \text{otherwise} \end{cases} $$然后对 $y(t)w(t)$ 进行采样，得到离散序列 $x[n]$。若记离散形式为 $x[n]=y[n]w[n]$，根据傅里叶变换相关定理，时域乘积对应频域卷积，即
$$ X(e^{j\omega}) = Y(e^{j\omega}) * W(e^{j\omega}) $$ 问题就出在这里：矩形窗的频域响应不是理想冲激，而是一个主瓣加若干旁瓣的 sinc 型结构。真实频谱与它卷积之后，本来应当集中在某一点的能量就会扩散到附近的频率位置。使用 Matlab 等工具可以直观看到矩形窗在频域上的形状：
下面这段是用于说明原理的 Matlab 演示脚本，用来观察矩形窗的时域与频域形状：
N = 64;
rectangular_window = ones(1, N);
rectangular_window(1, 1:16) = 0;
rectangular_window(1, 49:64) = 0;
fft_rectangular_window = fft(rectangular_window);
frequency_axis = linspace(-pi, pi - 2*pi/N, N);
figure;
subplot(2, 1, 1);
stem(0:N-1, rectangular_window);
title('矩形窗时域波形');
xlabel('采样点 (n)');
ylabel('幅度');
grid on;
subplot(2, 1, 2);
plot(frequency_axis, abs(fftshift(fft_rectangular_window))/N);
title('矩形窗频域幅度谱 (线性尺度)');
xlabel('归一化角频率 (弧度/采样)');
ylabel('幅度');
grid on;
xlim([-pi, pi]);
频谱泄露会带来几类实际问题：
- 频率分辨困难： 很接近的两个频率分量可能会因为能量泄露而混在一起，难以区分。
- 幅度估计偏差： 本来应落在主峰上的能量被分散到周围频点，导致读出来的主峰幅值偏小。
- 相位估计不稳定： 一旦目标频点附近混入较多泄露能量，相位也会跟着抖动。
- 出现虚假频率成分： 旁瓣有时会被误判成新的频率分量。
处理采样数据
为了减轻频谱泄露，常见的做法有以下几种：
加窗（Windowing）
对时域信号乘以窗函数，以减小截断效应。
- 常用窗函数包括：汉宁窗（Hanning）、汉明窗（Hamming）、布莱克曼窗（Blackman）、凯泽窗（Kaiser）、布莱克曼-哈里斯窗（Blackman-Harris）。
- 不同窗函数在“主瓣宽度”和“旁瓣抑制”之间会有不同权衡。
零填充（Zero Padding）
在原始信号后补零，可以让 FFT 频谱显示得更平滑，更容易观察峰值位置。需要注意的是，零填充并不会真正提高系统的物理分辨率，它更多起到频谱插值的作用，方便后续做峰值搜索和频率估计。
尽量截取整数周期
如果在截取信号时，能尽量保证采样窗口包含整数个周期，那么频谱泄露会显著减轻。不过在未知频率的测量场景里，这件事通常很难做到。
去除直流分量
实际采样中，ADC 偏置、前级电路失调等都可能引入直流分量。若不先去掉均值，频谱低频端往往会被抬高，影响峰值搜索。因此在 FFT 之前，通常应先减去均值。
在嵌入式场景下，加窗通常是最稳妥也最通用的做法。尤其是使用针对 Arm 平台优化后的 CMSIS-DSP 库进行 FFT 时，FFT 点数一般受到固定接口的限制，零填充不一定方便实现，因此实际工程里往往更依赖加窗。这里给出几种常用窗函数的生成代码。为了提高运行速度，下面使用 arm_math 中的三角函数接口：
下面这段是说明性示意代码，用于展示几种常见窗函数的生成方式：
// FFT_SIZE 是要进行 FFT 的点数
// 汉宁窗
for(i = 0; i < FFT_SIZE; ++i)
window[i] = 0.5f * (1.0f - arm_cos_f32(2.0f * M_PI * i / (FFT_SIZE - 1)));
// 汉明窗
for(i = 0; i < FFT_SIZE; ++i)
window[i] = 0.54f - 0.46f * arm_cos_f32(2.0f * M_PI * i / (FFT_SIZE - 1));
// 布莱克曼窗
for(i = 0; i < FFT_SIZE; ++i)
window[i] = 0.42323f - 0.49755f * arm_cos_f32(2.0f * M_PI * i / (FFT_SIZE - 1)) + 0.07922f * arm_cos_f32(4.0f * M_PI * i / (FFT_SIZE - 1));
// 布莱克曼-哈里斯窗
for(i = 0; i < FFT_SIZE; ++i)
window[i] = 0.35875f - 0.48829f * arm_cos_f32(2.0f * M_PI * i / (FFT_SIZE - 1)) + 0.14128f * arm_cos_f32(4.0f * M_PI * i / (FFT_SIZE - 1)) - 0.01168f * arm_cos_f32(6.0f * M_PI * i / (FFT_SIZE - 1));对于本文这个问题，窗函数的选择也有实际含义。由于两信号最小频差是 $100Hz$，而我们的目标并不是仅靠 FFT 直接读出最终频率，而是先找到两个主峰附近的 bin，再交给后续算法细化。因此更关心的是“峰值不要被旁瓣污染得太严重”。在这种情况下，像 Blackman-Harris 这类旁瓣抑制能力较强的窗往往更适合做前端粗定位。
在实际工程里，仅仅“生成窗函数”还不够。更关键的是在 FFT 前先完成减均值和乘窗，否则直流偏置和边界截断会同时干扰频谱。仓库中的 FFT_start() 实现大致如下：
下面这段是节选自仓库 Drivers/FFT/FFT.c 的实际工程代码，对应 FFT 前处理与频谱计算流程：
void FFT_start(uint8_t window_type)
{
Init_window(window_type);
float32_t sum_avr = 0.0f;
for(int i = 0; i < FFT_SIZE; ++i)
{
ADC_ConvData[i] = (float32_t)ADCbuff[i] * 3.3f / 4096.0f;
sum_avr += ADC_ConvData[i];
}
sum_avr /= FFT_SIZE;
for(int i = 0; i < FFT_SIZE; ++i)
{
ADC_ConvData[i] -= sum_avr; // 去直流
fft_inputbuf[i] = ADC_ConvData[i] * window[i]; // 乘窗
}
arm_rfft_fast_instance_f32 Rfft;
arm_rfft_fast_init_f32(&Rfft, FFT_SIZE);
arm_rfft_fast_f32(&Rfft, fft_inputbuf, fft_outputbuf, 0);
arm_cmplx_mag_f32(fft_outputbuf, mag, FFT_SIZE);
}这段代码对应的其实正是前面提到的那两件事：先把采样值转成电压并减去均值，再乘窗后送进 FFT。这样得到的频谱图会干净很多，后面的峰值搜索和相位测量也会更稳定。
得到频率 bin 索引
DFT 将时域信号转换到频域后，输出仍然是离散的，这些离散点就分布在所谓的频率“桶”（frequency bins）中。每一个 bin 对应一个离散频率，而相邻两个 bin 之间的间隔就是频率分辨率。
想象一下，在一个固定时间窗口 $T$ 内观察信号。这个时间窗口内，除了直流分量 $0Hz$ 之外，能够完整放下的最低频率，就是周期恰好等于 $T$ 的那个频率，也就是
$$ f_{min} = \frac{1}{T} $$这个 $1/T$，就是 FFT 频率栅格的间距。也就是说，如果只靠 FFT 的 bin 中心来读数，那么频率只能落在这些离散栅格上，这正是单帧 FFT 精度受限的根源。
显然，这个时间窗口就是采样总时长，即 $T = NT_s = N/f_s$，将其代入后可得
$$ \Delta f = \frac{f_s}{N} \tag 1 $$在后续算法中，这个公式会非常重要：FFT 负责给出“离真实频率最近的那个 bin”，而跨帧相位差法则负责把这个离散 bin 进一步细化为连续频率。
在仓库实现中，粗定位采用的是一个非常直接但够用的策略：先找到最大峰，再在它附近留出一小段保护区，随后搜索第二大峰。对应代码如下：
下面这段是节选自仓库 Drivers/FFT/FFT.c 的实际工程代码，对应双峰粗定位逻辑：
void find_peaks(uint32_t *k1, uint32_t *k2)
{
float32_t m = 10;
int i = 0;
for(int k = 0; k < FFT_SIZE / 2; ++k)
{
if(mag[k] > m)
{
m = mag[k];
i = k;
}
}
*k1 = i;
i = 0;
m = 10;
for(int k = 0; k < FFT_SIZE / 2; ++k)
{
if(k > *k1 - 4 && k < *k1 + 4)
continue;
if(mag[k] > m)
{
m = mag[k];
i = k;
}
}
*k2 = i;
}这种写法不算“通用峰值检测器”，但对于本文这个已知只有两路正弦、且两者频差不太小的任务来说，反而很实用。它的优点是开销小、逻辑直接，也便于在 STM32 上实时运行。
跨帧相位差法测频
完成加窗和 FFT 之后，我们可以从频谱图中找到每个信号对应的主峰 bin 索引。到这里，FFT 完成的是“粗定位”，而真正把频率精度继续往上提的，是跨帧相位差法。
为了推导方便，先只考虑单个正弦分量。设采样得到的信号为
$$ x[n] = A\sin(2\pi fnT_s + \varphi_0) \tag 2 $$其中 $T_s = \frac{1}{f_s}$ 为采样周期，$f_s$ 为采样频率。频谱图中的频率分辨率为 $\Delta f = f_s/N$，因此对于真实频率 $f$，我们定义它对应的理想 bin 索引为
$$ m = \frac{f}{\Delta f} \tag 3 $$这里的 $m$ 一般不是整数，因为真实频率几乎不可能恰好落在 FFT 的某个栅格中心上。于是有
$$ f = m \Delta f = m \frac{f_s}{N} \tag 4 $$因此，只要能估计出 $m$ 的小数部分，就能在粗定位 bin 的基础上把频率进一步细化。跨帧相位差法，本质上就是用相邻两帧在同一个 bin 处的相位变化来估计这部分细小偏差。
为了推导更清楚，不妨先用复指数形式来观察这个问题。考虑信号的正频率分量
$$ s[n] = Ce^{j(2\pi fnT_s + \varphi_0)} \tag 5 $$其中 $C$ 为复常数。对第 $r$ 帧信号（起始采样点为 $rN$）做 $N$ 点 DFT，其第 $k$ 个系数可写为
$$ X_r[k] = \sum_{n = 0}^{N-1} s[n+rN]e^{-j2\pi kn/N} \tag 6 $$将式（5）代入，有
$$ \begin{align*} X_r[k] &= \sum_{n=0}^{N-1} Ce^{j(2\pi f(n+rN)T_s+\varphi_0)}e^{-j2\pi kn/N} \\ &= Ce^{j(2\pi frNT_s+\varphi_0)} \sum_{n=0}^{N-1} e^{j2\pi (f/f_s-k/N)n} \tag 7 \end{align*} $$如果对每一帧都施加同一个窗函数，那么上式中的求和部分还会乘上相同的窗权重。关键在于：对于相邻帧来说，这一项基本不变，真正发生变化的是前面的整体相位因子。因此可以把它写成
$$ X_r[k] = B(k,f)e^{j(2\pi frT+\varphi_0)} \tag 8 $$其中 $T=NT_s$ 是一帧时长，$B(k,f)$ 是由窗函数和频率偏移共同决定的复系数。它会影响幅值，也会给相位带来一个固定偏置，但对于相邻两帧而言，只要我们观察的是同一个信号、同一个 bin，这个偏置基本保持不变。
因此，相邻两帧在同一 bin 上的相位差满足
$$ \Delta\varphi = \angle X_{r+1}[k] - \angle X_r[k] \approx 2\pi fT \pmod{2\pi} \tag 9 $$也就是说，一帧时间里信号本身转过了多少角度，频域系数的相位就会跟着转过多少角度。于是理想情况下有
$$ f = \frac{\Delta\varphi}{2\pi T} \tag {10} $$但真正测出来的相位差是“包裹”在区间 $[-\pi,\pi]$ 内的，也就是说我们拿到的是
$$ \Delta \varphi_{mea} = \Delta \varphi \pmod{2\pi}, \qquad \Delta \varphi_{mea} \in [-\pi,\pi] \tag {11} $$这看起来似乎丢掉了整圈数信息，但由于 FFT 已经给出了粗略 bin 索引 $k$，我们其实已经知道真实频率大约就在 $k\Delta f$ 附近。于是可写成
$$ f = k\Delta f + \delta f = \frac{k}{T} + \delta f,\qquad |\delta f| < \frac{1}{2T} \tag {12} $$于是
$$ \begin{align*} \Delta \varphi &= 2\pi fT \\ &= 2\pi\left(\frac{k}{T}+\delta f\right)T \\ &= 2\pi k + 2\pi \delta f T \end{align*} $$对它取模 $2\pi$ 之后，整数圈 $2\pi k$ 会被消去，于是得到
$$ \Delta \varphi_{mea} = 2\pi \delta f T \tag {13} $$从而可解出细化频差
$$ \delta f = \frac{\Delta \varphi_{mea}}{2\pi T} \tag {14} $$因此真实频率可以写成
$$ f = \frac{k}{T} + \frac{\Delta \varphi_{mea}}{2\pi T} \tag {15} $$这就是跨帧相位差法最核心的结论。它的直觉其实很简单：FFT 负责告诉我们“频率大概落在哪个格子”，而相位差负责告诉我们“它在这个格子里偏左还是偏右、偏了多少”。因此，虽然单帧 FFT 的频率分辨率是有限的，但借助相邻帧相位变化，我们可以把频率估计进一步细化到远小于一个 bin 间隔的尺度。
当然，这个结论成立也有前提：
- 同一个信号在相邻两帧中不能发生剧烈变化；
- 选取的 $k$ 应当是该信号对应的最近主峰 bin；
- 两帧之间的相位差在展开后应与该 bin 对应的频率邻域匹配；
- 信噪比不能太低，否则相位本身会抖动得很厉害。
在本文所面对的场景中，两路信号频率稳定、频差不小于 $100Hz$，因此这一方法效果很好。经实际测试，测频平均误差可以做到约 $0.002\%$。
下面给出仓库中更完整的主流程代码。需要注意的是，要连续采两帧数据并运行两次，这样才能得到跨帧的相位差：
下面这段是节选自仓库 Drivers/FFT/FFT.c 的实际工程代码，对应“FFT 粗定位 + 跨帧相位差修正”的主流程：
void process_signal(void)
{
FFT_start(BLACKMAN_HARRIS);
uint32_t k1 = 0, k2 = 0;
find_peaks(&k1, &k2);
float32_t f1 = (float32_t)k1 * (float32_t)SAMPLE_RATE / (float32_t)FFT_SIZE;
float32_t f2 = (float32_t)k2 * (float32_t)SAMPLE_RATE / (float32_t)FFT_SIZE;
float32_t *c1 = &fft_outputbuf[k1 * 2U];
float32_t *c2 = &fft_outputbuf[k2 * 2U];
float32_t phi1_now, phi2_now;
arm_atan2_f32(c1[1], c1[0], &phi1_now);
arm_atan2_f32(c2[1], c2[0], &phi2_now);
float32_t frameT = (float32_t)FFT_SIZE / (float32_t)SAMPLE_RATE;
if (prev[0].k == k1)
{
float32_t phi_prev;
arm_atan2_f32(prev[0].im, prev[0].re, &phi_prev);
float32_t delta_phi = phi1_now - phi_prev;
if (delta_phi > M_PI) delta_phi -= 2.0f * M_PI;
if (delta_phi < -M_PI) delta_phi += 2.0f * M_PI;
f1 += delta_phi / (2.0f * M_PI * frameT);
}
if (prev[1].k == k2)
{
float32_t phi_prev;
arm_atan2_f32(prev[1].im, prev[1].re, &phi_prev);
float32_t delta_phi = phi2_now - phi_prev;
if (delta_phi > M_PI) delta_phi -= 2.0f * M_PI;
if (delta_phi < -M_PI) delta_phi += 2.0f * M_PI;
f2 += delta_phi / (2.0f * M_PI * frameT);
}
prev[0] = (bin_prev_t){c1[0], c1[1], k1};
prev[1] = (bin_prev_t){c2[0], c2[1], k2};
if(!k2)
f2 = 0.001f;
tones[0].f = f1;
tones[1].f = f2;
} 整个流程是这样串起来的：FFT_start() 负责前处理和频谱计算，find_peaks() 给出两个粗略主峰，随后主流程在对应 bin 上提取复数频谱值，计算跨帧相位差，并将修正后的频率写回 tones[] 结构体中。
上面这段代码还有几个值得注意的工程细节。
第一，只有当上一帧和当前帧识别到的是同一个 bin 时，才能直接做相位差修正，否则就要重新建立跟踪关系。第二，相位差本身必须做包裹处理，也就是把结果限制在 $[-\pi,\pi]$，否则一旦跨越相位边界，就会把本来很小的相位变化误判成接近 $2\pi$ 的大跳变。第三，仓库里还保留了 if(!k2) f2 = 0.001f; 这样的保护逻辑，本质上是在第二峰未可靠识别时，避免后续最小二乘求解直接落到退化情形。
最小二乘法测幅测相
在测频阶段已经得到了精度较高的频率估计，因此接下来就可以把频率视为已知量，转而用最小二乘法去估计幅度和相位。这个思路本质上是“先测频，再做参数拟合”，相比直接在频率、幅度、相位三个维度同时搜索，会轻量得多，也更适合单片机实现。
对于双正弦混合信号，我们建立模型
$$ x[n] = A_1\cos(\omega_1 n+\varphi_1)+A_2\cos(\omega_2 n+\varphi_2)+\varepsilon[n] \tag 16 $$其中 $\omega_1,\omega_2$ 已通过 FFT 和跨帧相位差法标定，$\varepsilon[n]$ 为残差项。利用三角恒等变换，有
$$ \begin{align*} A_k\cos(\omega_k n+\varphi_k) &= A_k[\cos(\varphi_k)\cos(\omega_k n)-\sin(\varphi_k)\sin(\omega_k n)] \\ &= I_k\cos(\omega_k n)-Q_k\sin(\omega_k n) \tag 17 \end{align*} $$其中 $I_k = A_k\cos(\varphi_k)$，$Q_k = A_k\sin(\varphi_k)$。于是原式可以改写为
$$ \begin{align*} x[n] &= I_1C_{1n} - Q_1S_{1n} + I_2C_{2n} - Q_2S_{2n} + \varepsilon [n] \\ &= \boldsymbol{h}_n^T\boldsymbol{\theta} + \varepsilon[n] \tag 18 \end{align*} $$其中
$$ \begin{align*} &\boldsymbol{h}_n = [C_{1n}, -S_{1n}, C_{2n}, -S_{2n}]^T \\ &\boldsymbol{\theta} = [I_1,Q_1,I_2,Q_2]^T \\ &C_{kn} = \cos(\omega_k n), \qquad S_{kn} = \sin(\omega_k n) \end{align*} $$把所有采样点堆叠起来，就得到矩阵形式
$$ \boldsymbol{X} = \boldsymbol{H\theta} + \boldsymbol{\varepsilon}, \qquad \boldsymbol{H}\in \mathbb{R}^{N \times 4} \tag 19 $$其中
$$ \begin{align*} &\boldsymbol{X} = \Big[x[0], x[1], \dots, x[N-1]\Big]^T \\ &\boldsymbol{H}= \begin{bmatrix} C_{10} &-S_{10} &C_{20} &-S_{20} \\ C_{11} &-S_{11} &C_{21} &-S_{21} \\ {} &{} \vdots \\ C_{1,N-1} &-S_{1,N-1} &C_{2,N-1} &-S_{2,N-1} \end{bmatrix}_{N \times 4} \\ &\boldsymbol{\varepsilon} = \Big[\varepsilon [0],\varepsilon [1], \dots,\varepsilon [N-1] \Big]^T \end{align*} $$最小二乘估计定义为
$$ \hat\theta = \arg\min_\theta ||\boldsymbol{X}-\boldsymbol{H\theta}||^2=\arg\min_\theta \sum_n \varepsilon[n]^2 \tag 20 $$对其求导并令导数为零，可得正规方程
$$ \boldsymbol{H}^T\boldsymbol{H}\boldsymbol{\hat\theta}=\boldsymbol{H}^T\boldsymbol{X} \tag 21 $$这就是最小二乘法对应的标准线性方程组。由于本文中只有两路正弦，因此未知参数只有 $I_1,Q_1,I_2,Q_2$ 四个，最终只需求解一个 $4\times4$ 线性系统。求解之后即可得到
$$ \begin{align*} \boldsymbol{\hat\theta} &= [I_1,Q_1,I_2,Q_2]^T \\ &= [A_1\cos{\varphi_1}, A_1\sin{\varphi_1},A_2\cos{\varphi_2},A_2\sin{\varphi_2}]^T \tag 22 \end{align*} $$于是幅度和相位分别为
$$ A_k = \sqrt{I_k^2 + Q_k^2} $$$$ \varphi_k = \operatorname{atan2}(Q_k, I_k) $$ 这里使用 atan2 而不是普通反正切，是为了保留象限信息，避免相位落在错误的象限中。下面给出对应的最小二乘代码实现：
下面这段是节选自仓库 Drivers/FFT/FFT.c 的实际工程代码，对应双频最小二乘拟合：
/**
* @brief 最小二乘法计算幅度与相位
* @note 无
* @param f1, f2: 已确定频率
* @param x: 采样序列
* @param I1, Q1, I2, Q2: 求解参数组
* @retval 无
*/
void least_square(float32_t f1, float32_t f2, const float32_t *x, float32_t *I1, float32_t *Q1, float32_t *I2, float32_t *Q2)
{
// 计算一次角增量
float32_t w1 = 2.0f * M_PI * f1 / SAMPLE_RATE;
float32_t w2 = 2.0f * M_PI * f2 / SAMPLE_RATE;
// 迭代生成sin, cos序列
float32_t c1 = arm_cos_f32(w1), s1 = arm_sin_f32(w1);
float32_t c2 = arm_cos_f32(w2), s2 = arm_sin_f32(w2);
float32_t cn1 = 1, sn1 = 0, cn2 = 1, sn2 = 0;
float32_t Scc1 = 0,Sss1 = 0,Scc2 = 0,Sss2 = 0;
float32_t Scc12 = 0,Scs12 = 0,Ssc12 = 0,Sss12 = 0;
float32_t SxC1 = 0,SxS1 = 0,SxC2 = 0,SxS2 = 0;
for (uint32_t n = 0; n < FFT_SIZE; ++n)
{
float32_t xn = x[n];
// 同频能量
Scc1 += cn1 * cn1;
Sss1 += sn1 * sn1;
Scc2 += cn2 * cn2;
Sss2 += sn2 * sn2;
// 跨频交叉项
Scc12 += cn1 * cn2;
Scs12 += cn1 * sn2;
Ssc12 += sn1 * cn2;
Sss12 += sn1 * sn2;
// 数据投影
SxC1 += xn * cn1;
SxS1 += xn * sn1;
SxC2 += xn * cn2;
SxS2 += xn * sn2;
// 迭代sin/cos
float ncn1 = cn1 * c1 - sn1 * s1;
float nsn1 = sn1 * c1 + cn1 * s1;
float ncn2 = cn2 * c2 - sn2 * s2;
float nsn2 = sn2 * c2 + cn2 * s2;
cn1 = ncn1;
sn1 = nsn1;
cn2 = ncn2;
sn2 = nsn2;
}
// 构建正交矩阵和投影向量，部分交叉项可近似为0
float32_t H[4][4] = {
{Scc1, 0, Scc12, Scs12},
{0, Sss1, Ssc12, Sss12},
{Scc12, Ssc12, Scc2, 0 },
{Scs12, Sss12, 0, Sss2 }
};
float32_t b[4] = {SxC1, SxS1, SxC2, SxS2};
// 高斯消元
for (int k = 0; k < 4; k++)
{
// 主元选择
int m = k;
float_t maxv = fabsf(H[k][k]);
for (int i = k + 1; i < 4; ++i)
{
if(fabsf(H[i][k]) > maxv)
{
m = i;
maxv = fabsf(H[i][k]);
}
}
// 交换行
for(int j = k; j < 4; ++j)
{
float32_t t = H[k][j];
H[k][j] = H[m][j];
H[m][j] = t;
}
float32_t tb = b[k];
b[k] = b[m];
b[m] = tb;
// 标准化主行
float32_t diag = H[k][k];
for (int j = k; j < 4; ++j)
H[k][j] /= diag;
b[k] /= diag;
// 消去下方得到上三角式
for (int i = k + 1; i < 4; ++i)
{
float_t factor = H[i][k];
for(int j = k; j < 4; ++j)
H[i][j] -= factor * H[k][j];
b[i] -= factor * b[k];
}
}
// 回代求解
float32_t theta[4];
for(int i = 3; i >= 0; --i)
{
float32_t sum = b[i];
for(int j = i + 1; j < 4; ++j)
sum -= H[i][j] * theta[j];
theta[i] = sum;
}
*I1 = theta[0];
*Q1 = theta[1];
*I2 = theta[2];
*Q2 = theta[3];
}这段实现里有一个隐含前提：两路频率已经估计得足够准确，否则基函数 $\cos(\omega_kn)$ 和 $\sin(\omega_kn)$ 本身就会发生失配，最终会把一部分频率误差“吃进”幅度和相位估计里。这也是为什么本文要先做高精度测频，再做最小二乘拟合，而不反过来。
虽然跨帧相位差法可以得到很高的测频精度，但幅值和相位的最小二乘估计仍然会受到采样噪声、ADC 量化误差、模拟前端失真，以及频率估计残差的共同影响。经笔者测试，其平均误差约为 $1.7\%$。这也比较符合直觉：测频主要依赖“相位随时间的累计变化”，对单点采样误差不算太敏感；而最小二乘拟合则直接使用整段波形上的每一个采样点，因此更容易受到采样误差和模型失配的影响。
结语
到这里，本文实际上已经完成了“参数意义上的分离”：也就是从混合信号中估计出两个分量各自的频率、幅度和相位。但如果目标是进一步实现某一路信号的实时重建与输出，让它在示波器上看起来像是被“数字滤波”出来的一样，那么还会遇到新的问题。
不论使用单片机内部 DAC，还是外接 DDS 模块，只要信号发生端与输入信号源不共用同一个时钟基准，就不可避免地会存在频率偏差和相位漂移。短时间内这种偏差可能并不明显，但随着时间增长，采样时刻与信号真实相位之间的相对关系会持续变化，于是就会出现相位差累积。其结果是：即使你已经在参数上“识别”出了某一路信号，当你尝试把它重新输出并和原信号一起观察时，仍然会发现两者在时间轴上逐渐滑移，无法长期稳定重合。
要解决这个问题，最直接的方法其实是让采样系统与信号发生系统共用时钟源；如果这一点做不到，那么就需要引入闭环控制，实时追踪输入信号相位，并动态调整输出信号，使两者保持锁定。笔者曾尝试使用 STM32 的 DAC 来模拟 DDS 输出，但受限于硬件性能和整体系统结构，引入闭环之后效果仍不够理想。若后续要继续往“实时同步重建”这个方向推进，可能还是需要更专门的 DDS 或锁相实现方案。
🍃本文中所涉及的算法已经使用STM32F407VET6单片机实现，工程文件已上传至GitHub，欢迎访问。