频域补零上采样_AURIX 学习笔记（12）频域法互相关实现超声测距

概述

上一篇 AURIX 学习笔记（11）外部中断以及基于时域互相关的超声测距介绍了 TC264D 实现超声波测距应用的结构、原理等，也提到时域法的不足：计算复杂度高、滤波手段少、没有充分利用单片机外设资源。频域法互相关利用 FFT 模块，补上了时域法的短板，可以取得比较良好的效果。本文介绍频域法相对频域法的改动——主要是相关计算的改动——以及在 TC264DA 上实现频域互相关的一些细节（浪费调试时间的坑）。

频域互相关原理

数学原理

上一次写过，离散互相关运算的时域表达：

并且，我们知道离散卷积运算的表达式是：

发现两个表达式唯一的不通是原信号

取

还是

，换句话说，离散互相关和离散卷积的唯一区别是原信号的一次

反褶。

所以还可以从卷积运算满足交换律推出互相关运算不满足交换律。并且很容易记住交换互相关的两个信号导致相关谱的反褶。

我们还知道：

时域卷积等价于频域乘积；
时域反褶等价于频域共轭；

所以，可以写出互相关的频域表达：

其中，

是离散傅里叶变换，

是离散傅里叶反变换，

是同维向量的对应项相乘，

是对复信号取共轭。

性能分析

原本

的互相关运算，将被转化为 2 次傅里叶变换、1 次傅里叶反变换，以及

的取共轭再相乘（注意，这里是复数相乘，一个复数乘是 4 个实数乘和 2 个实数加）。并且，很容易想到更多优化思路：

「参考信号频谱的共轭

」是个常量

，没必要反复算，算一次存起来即可；
实信号的频谱，正频率和负频率成分呈共轭对称关系，所以乘法只需要算一半，
也只需要存一半；
相对于 123.76 kHz 采样率（根据奈奎斯特采样定理，有效频率是一半），超声收发系统通带很窄，不到有效测量频率的
。而通带以外算不算都是 0。

于是，实际上每次只要算 1 次傅里叶变换、1 次傅里叶反变换，以及最多

频谱长度的共轭和复数乘法（乘，然后共轭生成负频率成分）。对比时域法的性能：

-O3 下，200 点参考的互相关谱 1 个点约 600 周期
1024 点 FFT/IFFT 5000 周期，算上拷贝，置零之类的，差不多 20000 周期能完成一次普通的互相关运算，但这相当于是 256 点参考的 1024 点互相关谱，
，是

的差不多 40 倍

为啥 1024 点 FFT 实现的是 256 点参考的 1024 点互相关谱呢？

理解这个需要在有限长信号下考虑互相关运算。互相关本质上是在计算两个信号在不同对齐方式下的内积，或者说，两个无限维向量在方向上的接近程度（假设它们的长度/幅度已经归一化）。所以对于实际上有限长，或者是窗口外全为零的信号来说，它们的窗口至少有一个公共非零点，才能得到非零的内积，才能得到非零的相关谱线。当

的

挪动参考信号时，互相关谱的第一点是参考信号的尾重合于信号的头，即

，把参考信号左移了参考信号长度少一个。所以教科书上写道 m 点参考和 n 点信号的互相关得到的是 m+n-1 点互相关谱。那么 1024 点的 DFT 自然得到 1024 点的互相关谱，m+n-1 不能大于 1024，否则由于 FFT 假设周期延拓，得到的是

混叠的互相关谱。而 m+n-1 较小则会在互相关谱里补零，有点浪费计算量。

所以 1024 点 FFT 是肯定不能计算 1024 点参考的互相关的，选择 256 的其他理由在后面。

知乎的分隔线好大！

分帧问题

当然，这种加速也是有代价的，代价就是「窗」的存在。时域互相关可以每一点都算一次，完全从定义出发。而频域则必须截取一段时域信号，才能享受这个加速的优势。截的越长，优势越大。

这个在软件互相关才是完全正确的，按理说 FFT 使得

加速到

，加速效果随 n 指数增长。但 AURIX 学习笔记（10）硬件 FFT 写过，TC264DA 的 FFT 模块差不多是线性的。但是（相同采样率下）更长的窗就是更多的点数，通常会使用更长的 FFT，由于频谱点数也变多，相当于频域分辨率加大，还是会获得一些精度提升。

有窗，就会遇到窗里有没有我要的信号？有多少？如果只有一半在窗里怎么办？等等问题。这些问题非常现实。为了解决窗里截到部分信号的问题，就需要重叠分帧。如图所示：

重叠分帧，保证至少有一帧包含完整的目标信号

只要满足这些条件：

帧长度不小于信号长度
两帧重叠部分长度不小于信号长度

就能保证至少有一帧包含完整的信号。

总结一下，现在我们的（帧、重叠部分、参考信号）长度选择受到以下限制：

而我们回头去考察时域互相关，会发现时域互相关其实就是

、

的而已，所以一次只能算出 1 个点也没什么奇怪的。这也告诉我们重叠长度越小，重复计算就越少，当然相当于性能就越高。

进一步考虑，为了让 CPU 多做事，不应该让 ADC 中断这样的例行公事频繁打扰 CPU，那就要使用 DMA，而又由于 AURIX DMA 的奇怪设计，其地址循环必须是 2 的整数次幂。

当然地址循环是，中断不一定非要和地址循环对应，这里仍有优化空间。可以改但没必要。

所以选择帧长度=768，参考长度=重叠长度=256 就显得比较自然、顺理成章了。

实现

为了使用频域法计算互相关，进行下列配置：

分配内存
- 供 ADC 的 DMA 写的 uint16 vadc_buffer[256]
- 保存历史数据的 uint16 copy_buffer[512]
- 保存参考时域信号的 const sint16 reference_signal[256]（编译器懂，这个在 ROM 里）
- 保存参考共轭频谱的 sint32 ref_buffer[1024]（1024 是 512 个复数，实际上可以在外面算频谱，不存时域的，让这个在 ROM 里。但是反正内存够用，这样显得原汁原味）
- 用于 FFT 计算的 sint32 fft_buffer[2048]（这才是 1024 个复数）
初始化外设
- VADC 模数转换模块、及其 DMA
- ERU 外部中断模块
- ASC/UART 串口通信模块
- FFT 模块
计算：计算参考信号频谱，把一半的共轭保存到内存对应区域

计算过程仍分超声信号线和触发信号线，触发信号到来时总结最高峰到上一个触发信号的时延，减去帧长度并打印到串口。超声信号线计算步骤：

等待 DMA 中断
在 DMA 中断服务中将最新信号复制到 copy_buffer+256，也就是历史数据的后一半。这一步要第一时间做，免得下一个 VADC 完成冲了这次的结果。所以在中断服务里做就不错。拷贝用 memcpy，比手写循环快！即使手工优化到 64 位拷贝一次 memcpy 还快 20% 以上！
将 copy_buffer+256 拷贝到 fft_buffer+512，fft_buffer 的前 512 点是上次的 copy_buffer，要在上次计算结束时准备好
求 fft_buffer 中这 768 点的均值，减去均值。VADC 给出的是12 位无符号数，FFT 要补 0，不减去直流分量会导致引入一个阶跃信号，阶跃信号有极宽的频谱，会污染信号频带，导致准确性下降。之所以 768 点求一次，而不是 DMA 之后 256 点求一次，是为了避免直流分量的不够平稳在一次计算中每 256 点减的是不同的数，这个可以优化减少加法计算，但时间够用，没必要
FFT！
由于我们采样率仅有 123.76 kHz，高频几乎已经占满，不必再滤了。用 memset 将频谱 0~39 kHz 的成分干掉，然后将负频率成分全部干掉。这就是没有相位延迟的理想带通滤波，比什么 FIR、IIR 高到不知哪里去了
从 39 kHz 开始，乘上参考信号的频谱，这一步可以做某种白化滤波、频域加权、均衡之类的操作
IFFT！
在这个窗口的互相关谱上寻峰，保存

实现中的坑

注意，这个计算步骤里有一个大坑，就是 AURIX 的硬件 FFT 是定点运算，固定的 sint32。所以要时刻警惕溢出！参考信号先在线下做个 FFT 看一下，要保证其 FFT 之后幅度谱小于

，这样才能保证完全对齐时信号频谱与参考信号频谱的乘积不溢出，这个实际上要求参考信号幅值要比较小，最好先写程序验证一番。

此外，还有一些小坑。

虽然：

/** brief Configuration structure for the FFT job*/
typedef struct
{IfxFft_Fft      *fft;                     /**< brief Pointer To FFT.Used when (IFXFFT_FFT_OPTIMIZED == 0) */uint16           inputLength;             /**< brief Length of the input */uint16           outputLength;            /**< brief Length of the output */IfxFft_Length    fftLength;               /**< brief Length of the transform */IfxFft_Input     inputFormat;             /**< brief Input format */IfxFft_Output    outputFormat;            /**< brief Output format */IfxFft_Operation operation;               /**< brief Operation (FFT / IFFT) */boolean          useWindowFunction;       /**< brief Selection to use window function */IFX_CONST void  *inputPtr;                /**< brief Pointer to input data */void            *outputPtr;               /**< brief Pointer to output data */
} IfxFft_Fft_JobConfig;

里有 inputLength 和 outputLength 两个字段，但这两个字段不一样会导致 FFT 计算无法完成，没看到相关描述，原因不明。对我影响不大，我就没仔细研究。

不要让编译器自动分配大块内存，用 __at(p) 指定。亲测

#pragma section farbss bss_cpu0
#pragma section farbss restore

没有卵用，还是会从 0x60000000 开始分配。而 Cpu0 对应的内存地址从 0x70000000 开始，访问 0x60000000 段会慢 3 倍。

实测效果

知乎视频www.zhihu.com

代码在 https://github.com/autolaborcenter/test_aurix 。和本文描述不完全一样，因为还在探索之中。但不一样的主要是过滤、优化，外设使用基本不会有什么变化。要做到这个视频的效果还需要一些对峰值的后过滤，目前还在优化，这次就不写了。说实在话，这些优化都不如干脆发大点的信号有用，所谓力大砖飞是也。有什么问题欢迎讨论。