1.前言

在上一节中,我们为fft模块简单的添加了抛物线插值算法,并且衡量了插值算法对不同频率差的幅度与频率的估计情况。然而,在与AI的多轮对线之下,我了解到: 不同窗函数适用的插值算法不同。这就意味着,原来设计的架构(基于算法进行的分类)并不适用。我们需要基于窗函数,对幅值矫正算法重新分类。同时,为矫正算法模块加入相位估计算法,以应对2025-G题(电路模型探究装置)这种需要相位谱(复数谱)的情形。

2.整体架构

2.1输出规范

得益于写抛物线插值时就考虑到了将来需要加入其他频率估计算法的情形。所以之前已经有了统一接口:

@dataclass
class EstimateData:
    freq_estimated: float
    amplitude_estimated: float

加入相位数据只需要在下方加入

phase_estimated: float

即可。

2.2 整体架构

如同《Head First 设计模式》中提到的, 程序的设计模式,应该面向修改关闭,面向拓展开放。面向不变关闭,面向变化打开。这有点类似与《矛盾论》中提到的"抓主要矛盾"的哲学思想。在本程序中,我们注意到,其变化量有两个: 窗函数(Window function), 插值算法。因此,我们应该面向这两个变化量留出拓展的余地,将窗函数和插值算法解耦,而面向其他无关紧要的细节,我们可以写成紧耦合的算法。但现在有以下几个问题:

原有的抛物线算法适用于所有窗函数,只不过估计精度不同。
所有窗函数适用的相位矫正算法都相同,但是频差与幅值估计算法均不同
一个窗函数可能有好几种不同的频差估计算法(如矩形窗就有rife和cadan两种频差估计算法),但不是所有算法都适用与某个窗函数。
原有的no_interp 和 parabolic_interp函数是紧耦合的,输入复数谱,频率bin,采样率,采样点数,相干增益后直接输出EstimateData,和新算法(将频率,幅值,相位估计算法拆分成不同函数)的写法有很大不同。首先,对于问题4,我们采取了简单直接的做法: 构造一个抽象的(通用的)函数estimate_freq_amplitude_phase(),并将其他频率估计函数设置为私有函数,以免产生错误调用。对于问题4,我在抽象函数estimate_freq_amplitude_phase()中,选择了使用选择分支来解决这个问题,因为旧的模式不会再增加,之后的代码架构一定是新的算法模式,因此我们直接采用了这种紧耦合的方法,不过这对我们后续的拓展毫无影响。其次,对于问题2-3,我们采用了较为复杂的策略模式来实现。选择策略模式是因为其能较为简单的在嵌入式C代码中实现(使用结构体或swich-case方式达到类似的效果) 我的办法是使用了四个字典来满足我的调用要求:每种窗函数只能调用对应的算法,不能调用适用于其他窗函数的算法。以下为字典列表由于幅值估计算法每一个窗只有一个,所以只用维护一个算法列表即可。

_AMP_METHODS = {
        "rect": _amp_rect,
        "hann": _amp_hann,
        "hamming": _amp_hamming,
        "blackman": _amp_blackman,
        "blackmanharris": _amp_blackmanharris,
        "flattop": _amp_flattop
        }

频差估计算法,这个用了三个算法字典,每个字典的作用分别为auto模式下调用的插值算法,不同窗函数能调用的插值算法表,所有算法对应的函数指针(python字典的值能设置为函数,这里姑且不严谨的叫它指针)

_DELTA_METHODS = {
        "rife": _delta_rect_rife,
        "cadan": _delta_rect_candan,
        "granke": _delta_hann_grandke,
        "offelli": _delta_hamming_offelli,
        "andria_blackman": _delta_blackman_andria,
        "agrez": _delta_blackmanharris_agrez,
        "andria_flattop": _delta_flattop_andria
        }

所有的插值算法列表

_AUTO_DELTA_METHODS = { #自动窗函数选择表
        "rect": "cadan",
        "hann": "granke",
        "hamming": "offelli",
        "blackman": "andria_blackman",
        "blackmanharris": "agrez",
        "flattop": "andria_flattop"
        }

auto模式下,不同窗函数选择的频率差估计算法

_WINDOW_CHECK_LIST = { #查看用户是否选择错了窗函数算法
        "rect": ["rife", "cadan"],
        "hann": ["granke"],
        "hamming": ["offelki"],
        "blackman": ["andria_blackman"],
        "blackmanharris": ["agrez"],
        "flattop": ["andria_flattop"]
        }

用于检测用户是否选择错了窗函数算法,如果选择错了,就报错接着是这些算法的实现,我说实话大部分都没有什么意思,基本上都是对着公式写代码,没有半点数据结构。其中δ为频率差,A为幅值。为了简便,我们做以下符号约定。 X为FFT后得到的复数谱,$\omega$ 为窗函数

符号	含义
$k_0$	幅度谱峰值索引
$r$	$\begin{equation}\begin{cases}1\quad X[k_0 + 1] ≥ X[K_0-1]\\-1\quad X[k_0 + 1] ≥ X[K_0-1]\end{cases}\end{equation}$
$\alpha$	$\frac{\\|X[k_0+r]\\|}{\\|X[k_0]\\|}$
$\delta$	频偏因子 ∈ [-0.5,0.5]
CG	相干增益 Coherent Gain,可以通过查窗函数表,或者算窗函数均值获得

定义函数$sinc(\delta) = \frac{\pi \delta}{sin(\pi \delta)}$

注:sinc函数就是矩形窗/冲激在频域上的幅频响应曲线

窗函数	算法名称	δ计算公式	A计算公式
rect	rife(幅度谱)	$\delta = r \cdot \frac{\alpha}{1+\alpha}$	$A = \frac {2\\|X[k_0]\\|}{N} \cdot sinc(\delta)$
rect	cadan(复数谱)	$\delta = \cal{R}\left\{\frac{X[k_0-1] - X[k_0+1]}{2X[k_0] - X[k_0-1] - X[k_0 + 1]}\right\}$	同上
Hann	Grandke	$\delta = r \cdot \frac{2\alpha - 1}{\alpha + 1}$	$A = \frac{2\\|X[k_0]\\|}{N \cdot CG} \cdot sinc(\delta) \cdot (1 - \delta^2)$
Hamming	Offelli & Petri	$\delta = r \cdot \frac{\alpha - 0.54}{0.46\alpha + 0.54}$	$A = \frac{2\\|X[k_0]\\|}{N \cdot CG} \cdot sinc(\delta) \cdot \frac{1 - \delta^2}{1 - 0.148 \delta^2} \approx \frac{3.704\\|X[k_0]}{N} \cdot sinc(\delta) \cdot (1 - \delta^2)$ ($\frac{2}{N \cdot CG} \approx 3.704$)
Blackman	Andira	$\delta \approx r \cdot \frac{3\alpha - 1.5}{\alpha + 1.5}$ ($\\|\delta\\| < 0.4$时$误差 < 0.001bin$)	$A = \frac{2\\|X[k_0]\\|}{N \cdot CG} \cdot sinc(\delta) \frac{(1 - \delta^2)(4 - \delta^2)}{4 - 1.69\delta^2} \approx \frac{4.762\\|X[k_0]\\|}{N} \cdot sinc(\delta) \cdot (1 - \delta^2)(1 - 0.25\delta^2)$ ($\frac{2}{N \cdot CG} \approx 4.762$)
Blackman-Harris	Agrez	$\delta \approx r \cdot \frac{4\alpha - 2}{\alpha + 2}$	$A \approx \frac{5.575\\|X[k_0]\\|}{N} \cdot sinc(\delta) \cdot (1 - \delta^2)(1 - 0.25\delta^2)(1 - 0.111\delta^2)$ ($\frac{2}{N \cdot CG} \approx 5.575$)
Flattop	Andria	$\delta \approx \frac{5.5\alpha - 4.5}{\alpha + 1}$ (粗近似)	$A \approx \frac{2\\|X[k_0]\\|}{N \cdot CG} \cdot (1 + 0.02\delta^2 + 0.0005\delta^4)$

注:Flattop由于mainlobe较宽所以算法只能粗略估计频差δ 注: 幅度矫正公式 $\frac{2|X[k_0]|}{N \cdot CG}$是实信号的情况(由于实信号的共轭对称性),如果是复信号,这一项应该变为 $\frac{|X[k_0]|}{N \cdot CG}$

相位矫正公式 $\phi = \angle X[k_0] - \pi \delta \cdot \frac{N - 1}{N}$

注:这里统一的抽象插值函数选择了输入复数谱,首先是因为对于矩形窗,Candan算法需要用到复数谱。其次是在MCU中,这样的好处是我们不需要另外一个幅度谱数组(只需要一个原始FFT复数谱数组),只需要在计算的时候计算复数模长即可。

2.3 对复频谱的处理

在调试过程中,发现了一个新问题:原本处理实信号的算法只取幅度谱$[\frac{N}{2},N-1]$这个区间内的点(索引从零开始),但是引入复数谱之后,就开始带来一个问题:以1024点为例,原始复数谱的点数为1024点(按numpy数组而非cmsis数组),而幅度谱只取了其中的512-1023这部分。其实如果只看实信号,这是一个非常简单的问题,只需要将输入幅度谱的bin索引+512即可解决问题。但是,之前所有FFT后处理模块,如查找频谱峰等,都具有查询双边谱的条件。为此,我决定:让插值模块输出原始的双边频率(正频率或负频率), 以添加对复信号处理的兼容性。首先要处理的是幅度矫正函数,由于实信号的共轭对称性,所以实信号的幅度矫正因子和复信号不一样(是复信号的两倍,具体看上面的注解) 因此,我们将之前公式中的所有系数2去掉。同时,因为原始的频率谱(未经过fftshift)时是这样一个状态,假设是1024点的fft,那么经过fft处理后,我们的频率正半轴$[0,\frac{f_s}{2}]$,在$[0:513]$ 这个序列当中,而$[513:]$ 为负频率序列。因此,当信号索引$k_0$大于512时,频率bin的计算公式变成 $f_s = (\frac{2}{N} - k_0) \cdot \frac{f_s}{N}$,带上频差估计$\delta$变为: $f_s = (\frac{2}{N} - k_0 + \delta) \cdot \frac{f_s}{N}$. 这里还是$+\delta$的原因是负半轴的坐标顺序是从$(-\frac{f_s}{2N},0)$,因此,在正半轴时,如果频偏是向正方向移动的,对应的就是频偏向负半方向移动,对应的$k'_{0} = k_0 - \delta$ ,代入上面的换算公式就能得到以上的式子原始的用于处理实信号(只有正频率)的freq_estimated如

estimated_freq = (peak_bin + estimated_freq_delta) * (fs_hz / n_pts)

那么这样,我们在此加入一个条件判定即可

if peak_bin > n_pts // 2:
	peak_bin = (n_pts // 2) - peak_bin
estimated_freq = (peak_bin + estimated_freq_delta) * (fs_hz / n_pts)

这样,我们就能正确的输出负频率了。对于正负频率的处理,我们应该将这个工作放到其他模块,而不是这个矫正算法模块中。

3.测试结果

仿照抛物线插值算法的测试流程,以下为测试项目:

扫频测试:绘制频偏-RMSE曲线,反应频率估计值,幅度估计值随频率偏差的变化,扫频测试中,生成的测试信号采样率为1Mhz,采样点数1024点,量化精度16bit,SNR=10,相位随机变化(因为根据相位矫正公式可以知道相位的误差跟原始相位无关)$\phi ∈ [0,2\pi]$测试范围为$\delta \in [0.0 , 0.5)$ 步长0.01(单位:bin)
SNR测试: 绘制SNR-RMSE曲线,反应频率估计值和幅度估计值随SNR的变化,生成的信号频偏$\delta = 0.3 bin$ ,采样率1Mhz,采样点数1024点,量化精度16bit,测试范围$SNR \in [0.0,60)$ ,步长5dB。