基于跨帧相位差法实现混合信号分离

前言

本文的问题背景是设计一个信号分离装置，可以分离由两个正弦波信号合成的混合信号，也就是说要测量两个原始信号的频率、幅度和相位，并尽可能提高精度。原始信号的频率范围为 $100Hz \sim 10kHz$ ，且信号的频率差最小为 $100Hz$ 。
设计的难点在于测量精度，对于分析混合信号，我们会首先想到快速傅里叶变换，即 FFT ，但是 FFT 要做到 5% 甚至 1% 的精度，必须将采样点数取得足够大，根据奈奎斯特采样定律，本设计中采样率最低应为 $20kHz$ ，为了留有一些容差，实际采样率应大于 $25kHz$ ，此时要做到 $1Hz$ 的测量精度，采样点数就应取 $25kHz / 1Hz = 25000$ ，而 FFT 一般要求采样点数为 $2$ 的幂次方，所以必须取到 $32768$ ，显然对于单片机来说这是一个超大的数组，一般 STM32 又以 uint16_t 类型存储 ADC 采样值，这时采样数组所需要的内存就有 $2 \times 32768Bytes = 64KB$ ，而进行 FFT 时又需要将数据进行零填充或转化为复数存储，所需要的空间又会成倍数增长，对于单片机来说这样的大小是无法做到的，所以我们需要采取其他方法来提高测量精度。
如标题所言，本文将介绍一种所谓跨帧相位差法的算法，理论上能够使测量频率无限逼近于真实频率，当频率测量足够准确时，幅度和相位就可以使用很简单的算法得到足够高的精度。在使用该算法之前，首先要得到频谱图中与信号真实频率最接近的频率 bin 索引。

在 FFT 之前

由傅里叶变换衍生而来的 快速傅里叶变换(Fast Fourier Transform, FFT) 能够将时域信息转化为频域信息，从而提供了分析信号的工具。理想的傅里叶变换需要求和整个时域的数据，即从0到无穷，而我们采样的点数是有限的，因此直接用采样得到的数据做 FFT 时，就会出现第一个问题：频谱泄露。

频谱泄露

专业地来说，当我们对一个时域信号进行离散傅里叶变换 (Discrete Fourier Transform, DFT) 以分析其频率成分时，如果这个信号不是周期性的，或者它的周期与我们进行DFT的窗口长度不完全匹配，那么信号的能量就会“泄露”到相邻的频率 bin (频率分辨率的最小单位) 中，原本应该集中在某个频率上的能量会扩散开来，形成一些不属于信号真实频率成分的“旁瓣”。

想象一下，你用一个固定长度的“窗户”去观察一段无限长的波浪。如果这个波浪的完整周期正好能被这个窗户的长度完美包含，那么你看到的频率成分就会很干净。但是，如果波浪的周期和窗户的长度不匹配，那么在窗户的边缘就会出现“截断”效应，这个截断就相当于给原始信号乘以了一个非理想的窗函数。

这个非理想的窗函数在频域上也有其对应的形状，它不是一个完美的冲击函数，而是具有一定的宽度和旁瓣。原始信号的真实频谱会和这个窗函数的频谱进行卷积，导致原始频谱的能量扩散到周围的频率上，形成了我们所说的频谱泄露。

从数学角度来讲，给定一个原始连续信号 $y(t)$ ，采样的时间为 $[0, T]$ ，那么采样得到的离散时间信号 $x[n]$ 就可以看做是 $y(t)$ 乘以一个矩形窗函数 $w(t)$ ：

$$ w(t) = \begin{cases} 1, \quad 0 \leq t \leq T \\ 0, \quad otherwises \end{cases} $$

然后对 $y(t) w(t)$ 采样得到 $x[n]$ 。不妨先令 $x[n] = y[n]w[n]$ ，根据傅里叶变换相关定理，时域的乘积的傅里叶变换就是频域的卷积，即 $X(e^{jw}) = Y(e^{jw}) * W(e^{jw})$ ，显然这不是我们想要的结果，事实上，矩形窗函数的频域函数为 $sinc$ 函数，使用 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)

在原始信号后添加零点(Zero Padding)，提高频率分辨率，使频谱更平滑，便于观察泄露现象。
选择整数周期采样

在截取信号时，尽可能让信号包含整数个周期，可显著减少频谱泄露。

其中加窗和增加采样点数是最常用的方法，但是如果要使用针对 arm 优化后的 CMSIS-DPS 库进行 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));

得到频率 bin 索引

DFT将时域信号转换到频域，频域上的输出也是离散的，分布在所谓的频率“桶”(frequency bins) 里。每个频率桶代表一个特定的频率。频率分辨率就是这些相邻频率桶之间的间隔。

想象一下，在一个固定的时间窗口 $T$ 内观察一个信号。在这个时间窗口内，能够被完整观察到的最低频率（除了直流分量 $0Hz$ 之外）就是其周期恰好等于 $T$ 的频率。这个频率就是 $f_{min} = 1/T$。

这个 $1/T$ 就是我们能分辨出来的最小频率变化量。任何频率变化如果小于这个值，在长度为 T 的观测窗口内就无法完整地展现出一个周期的差异，从而难以区分。

显然，这个时间窗口就是采样的总时长，即 $T = NT_S = N/f_s$ ，将其代入最小频率表达式，有

$$ \Delta f = \frac{f_s}{N} \tag 1 $$

跨帧相位差法测频

加窗并进行 FFT 之后，计算频域幅值得到频谱图，就可以得到每个信号的频率 bin 索引，至此完成了算法的前期准备工作，那么，接下来该怎么做呢？

首先，设信号为正弦波，那么时域中采样得到的数据可表示为

$$ x[n] = Asin(2\pi fnT_s + \varphi_0) \tag 2 $$

其中 $T_s = \frac{1}{f_s}$ 为采样周期，$f_s$ 为采样频率，我们知道频谱图中的频率分辨率为 $\Delta f = f_s/N$ ，那么对于真实频率 $f$ ，我们假定一个理想索引

$$ m = \frac{f}{\Delta f} \tag 3 $$

此处 $m$ 不一定为整数，因为无法保证所测频率一定卡在频谱图的栅格上。此时真实频率 $f$ 可以表示为

$$ f = m \Delta f = m \frac{f_s}{N} \tag 4 $$

因此我们只要找到这个 $m$ 的值，就可以精准得到真实频率。跨帧相位差法就是间接地计算 $m$ ，从而得到精度极高的测量频率。

我们先将正弦转换到复数域中：

$$ x[n] = \frac{A}{2j}[e^{j2\pi fnT_s + \varphi_0} + e^{-j2\pi fnT_s - \varphi_0}] \tag 5 $$

对一帧 $N$ 点实信号 $x[n]$ ，我们知道其第 $k$ 个 DFT 系数为

$$ X_k = \sum_{n = 0}^{N-1} x[n]e^{-j2\pi kn/N} \tag 6 $$

所以我们对正频率一侧做 DFT ，有

$$ \begin{align*} X_k &=\frac{A}{2j} \sum_{n = 0}^{N-1} e^{j2\pi fnT_s + \varphi_0}e^{-j2\pi kn/N} \\\\ &= \frac{A}{2j} \sum_{n = 0}^{N-1} e^{j\varphi_0}e^{j2\pi(n-\frac{N}{2})T_s+j2\pi f\frac{N}{2}T_s - j2\pi kn/N} \\\\ &= \frac{A}{2j} \sum_{n = 0}^{N-1} e^{j\varphi_0}e^{j2\pi(n-\frac{N}{2})(m-k)/N-j\pi (k-m)} \tag 7 \end{align*} $$

我们不需要关注求和部分，因此将求和无关项提出，这里需要注意如果加了窗函数，那么求和部分也会有窗函数的影响，

$$ X_k = AW(offset)e^{j\varphi_0}e^{-j\pi (k-m)} \tag 8 $$

其中 $W(offset)$ 是由窗函数和偏移量共同决定的复因子。

已知一帧长度为 $N$ ，采样周期为 $T_s$ ，可以得到一帧的时长

$$ T = NT_s \tag 9 $$

把同一信号往后拖一帧，即 $n' = n + N$ ，那么有

$$ \begin{align*} x[n'] &= Asin(2\pi fnT_s + 2\pi fNT_s + \varphi_0) \\\\ &= Asin(2\pi fnT_s + 2\pi fT + \varphi_0) \tag {10} \end{align*} $$

再次做 DFT ，容易得到

$$ X_k' \propto e^{j(\varphi_0 + 2\pi fT)} $$

可以看出，同一点的 DFT 系数，在时域信号延后了一帧之后，相位仅多了转过的角度 $2\pi fT$ ，不妨设 $\Delta \varphi = 2\pi fT$ ，即两点绝对相位差，则

$$ f = \frac{\Delta \varphi}{2\pi T} \tag {11} $$

显然，如果我们能得到这个绝对相位差，就能完全复原出真实频率。但是我们计算相位差的过程实际上是两点测量的相位作差，得到的结果 $\Delta \varphi_{mea} \in [-\pi, \pi]$ ，这个结果实际上没有加上转过的整圈数，也就是说

$$ \begin{align*} \Delta \varphi_{mea} &= \Delta \varphi \ \ mod \ \ 2\pi \\\\ &= 2\pi fT \ \ mod \ \ 2\pi \tag {12} \end{align*} $$

注意到我们之前标定的理想索引 $m$ ，

$$ 2\pi fT = \frac{2\pi fN}{f_s} = 2\pi m $$

设 $\Delta bi n = m-k$ ，则

$$ \begin{align*} \tag {13} 2\pi fT &=2k\pi + 2\pi \Delta bin \\\\ \frac{\Delta \varphi_{mea}}{2\pi T} &= \frac{2\pi fT \ \ mod \ \ 2\pi}{2\pi T} \\\\ &= f -\frac{k}{T} \tag {14} \end{align*} $$

所以真实频率为

$$ f = \frac{k}{T} + \frac{\Delta \varphi_{mea}}{2\pi T} \tag {15} $$

根据此算法，理论上可测得频率无限接近于真实频率，即使考虑到采样过程中存在的噪声以及计算的误差，也可以达到极高的精度，经笔者测试，平均误差仅有 $0.002\%$ 。

在这里给出实现算法的主要程序，注意要连续采两帧数据运行两次：

// 开启FFT
FFT_start(BLACKMAN_HARRIS);
	
// 找到两信号粗估计bin下标
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;	// 频移周期，即每帧间隔 4096 / 40k = 0.1024 s
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};
		
tones[0].f = f1;  tones[1].f = f2;

最小二乘法测幅测相

在测频算法中已经得到了精度较高的频率，因此使用最小二乘法测量幅度与相位能达到很好的效果。下面介绍一下最小二乘法：

对于

$$ x[n] = A_1cos(w_1+\varphi_1)+A_2cos(w_2n+\varphi_2)+\varepsilon[n] \tag 1 $$

其中 $w_1$ ，$w_2$ 已通过 FFT 和相位差法标定，$\varepsilon[n]$ 为残差，有

$$ \begin{align*} A_kcos(w_kn+\varphi_k) &= A_k[cos(\varphi_k)cos(w_kn)-sin(\varphi_k)sin(w_kn)] \\\\ &= I_kcos(w_kn)-Q_ksin(w_kn) \tag 2 \end{align*} $$

其中 $I_k = A_kcos(\varphi_k)$ ，$Q_k = A_ksin(\varphi_k)$ ，于是有

$$ \begin{align*} x[n] &= I_1C_{1n} - Q_1S_{1n} + I_{2n} - Q_2S_{2n} + \varepsilon [n] \\\\ &= \boldsymbol{h}_n^T\boldsymbol{\theta} + \varepsilon[n] \tag 3 \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(w_kn), \ S_{kn} = sin(w_kn) \end{align*} $$

将其写成矩阵形式，

$$ \boldsymbol{X} = \boldsymbol{H\theta} + \boldsymbol{\varepsilon}, \ \boldsymbol{H}\in R^{N \times 4} \tag 4 $$

其中，$N$ 为采样点数，

$$ \begin{align*} &\boldsymbol{X} = \Big[x[0], x[1], \dots, x[N-1]\Big]^T \\\\ &\boldsymbol{H}= \begin{bmatrix} C_{11} &-S_{11} &C_{21} &-S_{21} \\ C_{12} &-S_{12} &C_{22} &-S_{22} \\ {} &{} \vdots \\ C_{1N} &-S_{1N} &C_{2N} &-S_{2N} \end{bmatrix}_{N \times 4} \\\\ &\boldsymbol{\varepsilon} = \Big[\varepsilon [0],\varepsilon [1], \dots,\varepsilon [N-1] \Big]^T \end{align*} $$

最小二乘估计定义为：

$$ \hat\theta = argmin_\theta ||\boldsymbol{X}-\boldsymbol{H\theta}||^2=argmin_\theta \sum_n \varepsilon[n]^2 \tag 5 $$

令

$$ \frac{\partial ||\boldsymbol{X}-\boldsymbol{H\theta}||^2}{\partial \boldsymbol{\theta}} = -2\boldsymbol{H}^T(\boldsymbol{X}-\boldsymbol{H\theta}) = \boldsymbol{0} $$

得到

$$ \boldsymbol{H}^T\boldsymbol{H}\boldsymbol{\hat\theta}=\boldsymbol{H}^T\boldsymbol{X} \tag 6 $$

使用高斯消元法求解该线性方程组，即可得到

$$ \begin{align*} \boldsymbol{\hat\theta} &= [I_1,Q_1,I_2,Q_2]^T \\\\ &= [A_1cos{\varphi_1}, A_1sin{\varphi_1},A_2cos{\varphi_2},A_2sin{\varphi_2}]^T \tag 7 \end{align*} $$

此时即可计算出幅度与相位

$$ A_k = \sqrt{I_k^2 + Q_k^2} \\\\ \varphi_k = tan^{-1} \frac{Q_k}{I_k} $$

最小二乘代码实现如下：

/**
 * @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; 
		
        int test = 0;
        if(k == 3)
            test = 1;
		
        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];
}

虽然使用跨帧相位差法可以得到精度极高的频率，但是最小二乘法仍有一定的误差，经笔者测试平均误差约为 $1.7\%$ 。其原因可能是单片机 ADC 采样存在一定的误差，且信号本身存在噪声，导致采得的电压值存在误差，这种误差对于测频算法影响较小，但最小二乘法直接用到了每个采样点的信息，因此容易受采样误差的影响。

结语

以上算法仅测量分离了两个信号的频率、幅度和相位，而要实现“数字滤波”的效果，还需要将其中一个信号输出，这里涉及到与原信号同步稳定显示的问题。

不论使用单片机的 DAC 模块还是外部 DDS 模块，由于内部模块与外部信号源采用了各自独立的时钟基准，这两个时钟源之间不可避免地存在频率偏差和相位漂移。这种时钟异步性会导致在长时间观测下，单片机采样时刻与信号源实际相位之间的相对关系发生持续变化，即产生相位差的累积效应。累积效应作用的结果是，当通过滤波器对处理后的信号进行观测时（例如在示波器上显示），会发现滤波输出信号与原始输入信号之间存在一个逐渐增大的相位差，表现为两个波形在时间轴上相对滑动或漂移。

要解决该问题，最粗暴的方法其实是统一信号源与信号发生模块的时钟源，但大多数情况无法实现；其次可以通过引入闭环控制，实时追踪信号源的相位并跟随，做到相对静止。笔者尝试使用了 DAC 模拟 DDS 输出信号，由于硬件和性能原因，引入闭环控制后未能达到理想的效果，后续笔者会尝试使用专门的 AD 芯片来解决该问题。

🍃本文中所涉及的算法已经使用STM32F407VET6单片机实现，工程文件已上传至GitHub，欢迎访问。