基于矢量水听器的潜标探测系统——信号处理部分-编程阁

一、项目概述

基于矢量水听器的潜标探测系统_带水听器潜标-CSDN博客

前篇如上。

依旧是毕设项目，依旧是基于AD采集信号，将信号搬运到PS端的DDR3上，在ARM端裸机开发进行信号处理（FIR、FFT、矢量测向），将处理后的结果发送到上位机，并可以存储到SD卡内。

调了小一周，终于调好了。（累！）

硬件选择依旧是ZYNQ7010+AD7606并行采集。

二、PL端设计调整

基于矢量水听器的潜标探测系统_带水听器潜标-CSDN博客

根据前篇，我设计的是AD输出四路通道，后面接了四路FIFO，然后再接了四路DMA，最终通过HP接口将数据传到PS端的DDR3中。

周四我们组的硬件老师来了，我顺便问了一下，老师说就算是四个DMA，最终也是从一个HP接口传输数据，存到DDR3中反倒不合适，不如从一开始AD7606就交织输出，用一个FIFO缓存，再通过一个DMA搬运到DDR3中。

老师也建议，可以保留四路FIFO，但水声信号的数据量并不是很大，FIFO的深度可以减小一些，节省一些逻辑资源。同时，老师建议，ZYNQ7010的资源比较少，在PL端做多做到FIR，其余的挪到PS端来做。

2.1问题1：

我根据老师的建议做了一下，也就是AD7606→FIFO转AXI Stream→FIR→DMA IP核→DDR3，但在测试阶段，出现了解析错误，经查询是设置FIR通道的问题，如果我使用FIR多通道，例如我设置了四个通道的FIR，那如果我只输入三通道的实际数据，会发生数据错位。

原因：FIR IP核内部是严格按照TDM（时分复用）来处理交织数据的，如果只输入三路通道，从第四个时钟开始，CH0的数据会被IP核误判成CH3，会造成所有通道的滤波完全乱序，输出波形毫无规律，会造成数据的交叉污染。

所以我思来想去，还是决定在ARM端做FIR、FFT与矢量测向，虽然没有办法利用PL端的数据加速，但相较于在PL端做，也更加的方便。

2.2问题2：

还有一处修改：之前我灵机一动，在FIFO转AXI Stream模块中将数据拼接成了高16位全是0，低16位是数据，也就是在内存中是这个样子的：

内存: [0x0000_CH1] [0x0000_CH2] [0x0000_CH3] [0x0000_CH4] ...

在实际对信号进行信号处理过程中，这部分的解交织代码总是写不好，于是我决定保留最初的“两个16bit样本打包成一个32bit字”的数据格式。

原因：原始FIFO是16bit写、32读。AD7606每采一个16bit样本写进FIFO，FIFO内部每攒够2个16bit样本，就拼接成一个32字输出。也就是dout是32bit，里面同时包含了两个16bit AD样本。

assign M_AXIS_TDATA = {dout[15:0], dout[31:16]};

{dout[15:0], dout[31:16]}是半字交换（调整字节序），确保DMA写到DDR后，PS端解交织时拿到的通道顺序是对的。

故，我现在的PL端的BD图如下所示：

具体内容没怎么变，大家参考之前发的两篇文章就可以捏！

三、PS端信号处理代码实现

在之前的文章里，我已经在PS端利用lwIP实现TCP数据传输，以及SD卡的读写，这一次主要实现信号处理算法。

3.1解交织

想要在ARM端进行信号处理，要先对存入DDR3的数据进行解交织，故代码逻辑为：

uint32_t word = src[word_idx + w]; uint16_t sample_hi = (uint16_t)(word >> 16); // CH1, CH3... uint16_t sample_lo = (uint16_t)(word & 0xFFFF); // CH2, CH4...

解交织后，4路数据被拆分到独立地址：

CH1 →0x02000000
CH2 →0x02100000
CH3 →0x02200000
CH4 →0x02300000

3.2信号处理部分

在DMA传输完成后，主循环依次调用：

SignalProcess_ProcessAllChannels(actual_fs); // 4路FFT VectorProcessing_Execute(actual_fs, result); // 互谱+DOA send_spectrum_data(); // 发送频谱 VectorProcessing_SendResult(); // 发送方位结果

首先，我们要逐通道进行FFT处理，主要使用CMSIS-DSP的实数FFT来处理。

3.2.1初始化

先对FFT进行初始化：

arm_status status = arm_rfft_fast_init_f32(&S, FFT_SIZE);

作用是预计算旋转因子，填充结构体S内部查找表，这部分只执行一次，之后所有的1024点FFT复用该实例，避免重读计算cos/sin。

再对FIR滤波器进行初始化：

arm_fir_init_f32(&S_fir[i], FIR_NUM_TAPS, fir_coeffs, fir_state[i], FFT_SIZE);

我设计的FIR滤波器参数如下（仅作为参考，请大家结合自己的实际来配置参数）：

20ksps四通道带通FIR滤波器： 参数设计： 响应类型：带通 设计方法：FIR等波纹 采样率：20000HZ 通带下边界：50HZ，通带上边界：5000HZ 阻带下边界：20HZ或30HZ，阻带上边界：5500HZ 通带纹波：1dB 阻带衰减：60dB 阶数：128阶

这么设计滤波器的好处在于：1）50HZ以下的噪声能压掉，2）5KHZ内有效噪声基本保留，3）四通道一致性比较好

其中，fir_state是延迟线缓冲区，大小 =FIR_NUM_TAPS + FFT_SIZE = 129 + 1024 = 1153。

3.2.2单通道FFT

在单通道FFT中，第一步，要将ADC原始值转换成电压。

int16_t raw = (int16_t)adc_data[i]; // 16bit 有符号 fft_input[i] = ((float32_t)raw / 32768.0f) * 5.0f;

AD7606是±5V、16bit双极性，故+32767对应＋5.0V，-32768对应-5.0V，故应该按照代码的方式进行转换。

再进行FIR滤波：

arm_fir_f32(&S_fir[ch_idx], fft_input, filtered_input, FFT_SIZE);

再进行RMS计算：

for (int i = 0; i < FFT_SIZE; i++) { sum_sq += filtered_input[i] * filtered_input[i]; } arm_sqrt_f32(sum_sq / FFT_SIZE, &result->rms_value);

这一步必须在FFT之前，我将这部分挪到FFT之后做，输出的值与真实值不符合，后来才知道原因：因为arm_rfft_fast_f32会原地破坏输入缓冲区的filtered_input，把它当成临时工作区。

再执行FFT：

arm_rfft_fast_f32(&S, filtered_input, fft_output, 0);

参数0代表着正变换（也就是时域→频域），输入1024点实数时域信号，输出打包数据fft_output[1024]。

其中，arm_rfft_fast_f32对1024点实数输入的输出格式如下：

RFFT 输出格式
数组索引	内容	物理意义
`fft_output[0]`	`R0`	DC 分量（实部）
`fft_output[1]`	`R512`	Nyquist 点（实部）
`fft_output[2]`	`R1`	频点 1 实部
`fft_output[3]`	`I1`	频点 1 虚部
`fft_output[4]`	`R2`	频点 2 实部
`fft_output[5]`	`I2`	频点 2 虚部
...	...	...
`fft_output[1022]`	`R511`	频点 511 实部
`fft_output[1023]`	`I511`	频点 511 虚部

DC和Nyquist两个纯实数点被CMSIS塞到了最前面，[0]是DC，[1]是Nyquist，从[2]开始才是正常的（R1，I1）交错。

再计算幅度谱，将其从双边谱转换为单边谱：

// DC 分量 fft_magnitude[0] = fabsf(fft_output[0]) / FFT_SIZE; // 中间频点（1 ~ 511） for (int i = 1; i < FFT_SIZE/2; i++) { float32_t real = fft_output[2*i]; float32_t imag = fft_output[2*i + 1]; fft_magnitude[i] = (sqrtf(real*real + imag*imag) * 2.0f) / FFT_SIZE; } // Nyquist 点 fft_magnitude[FFT_SIZE/2] = fabsf(fft_output[1]) / FFT_SIZE;

计算其频率分辨率：

float32_t freq_resolution = sample_rate / FFT_SIZE;

假设采样率Fs = 2000 Hz，则分辨率为2000/1024 ≈ 1.953 Hz/ bin，峰值频率为peak_idx × 1.953 Hz。

对峰值进行检测

arm_max_f32(&fft_magnitude[1], FFT_SIZE/2 - 1, &peak_mag_temp, &peak_idx); peak_idx += 1;

这一步跳过了DC，避免0HZ直流干扰。

互谱计算的代码如下：

// 互谱实部：Re{ P*(f) · V(f) } = p_real*v_real + p_imag*v_imag float32_t cross_real = p_real * v_real + p_imag * v_imag;

这一步是矢量水听器DOA估计的数学基础，p是声压通道（标量），V是振速通道（矢量，分为Vx/Vy/Vz），P* · V的实部正比于声强的流向，Ix = Re{P* · Vx}表示X方向的声能流，Iy = Re{P* · Vy}表示Y方向的声能流。

DOA估计代码：

azimuth_rad = atan2f(Iy_val, Ix_val);

声能流再水平面的投影方向就是声源方位角。

总结下来，FFT的数据流如下所示：

ADC原始值(uint16_t)
↓ 强制转 int16_t，映射到 ±5V
时域电压(float32_t[1024])
↓ arm_fir_f32 滤波
滤波后信号(float32_t[1024])
↓ arm_rfft_fast_f32 (正变换)
打包频域数据(float32_t[1024])
├─ [0] = DC (R0)
├─ [1] = Nyquist (R512) ← 你代码里复数提取漏了这步
└─ [2..1023] = (R1,I1)...(R511,I511)
↓ 解析 + 幅度计算
单边幅度谱(float32_t[513])
├─ [0] = DC 幅度
├─ [1..511] = 频点1~511幅度 (×2归一化)
└─ [512] = Nyquist 幅度
↓ arm_max_f32 搜峰
峰值频率/幅度 + RMS
↓ 存入 ch_results[4]
TCP发送 / 互谱计算 / SD存储

3.3矢量处理与DOA估计

首先统一目标频率，以Ix.peak_idx（互谱峰值频点）作为统一的目标频率bin，读取该频点处的互谱值：

float32_t Ix_val = Ix.cross_spectrum[target_bin]; float32_t Iy_val = Iy.cross_spectrum[target_bin]; float32_t Iz_val = Iz.cross_spectrum[target_bin];

根据公式tan(θ) = Iy / Ix来计算声能流在X-Y平面的夹角，也就是水平方向的方位角：

azimuth_rad = atan2f(Iy_val, Ix_val); azimuth_deg = azimuth_rad * 180.0f / PI; if (azimuth_deg < 0) azimuth_deg += 360.0f;

再根据垂直分量与水平合成幅度的夹角计算垂直方向的俯仰角：

horizontal_mag = sqrtf(Ix_val*Ix_val + Iy_val*Iy_val); elevation_rad = atan2f(fabsf(Iz_val), horizontal_mag);

3.3数据输出

3.3.1频谱数据发送

send_spectrum_data()为每个有效通道发送一帧：

帧头：0x55 0xA5+ 命令字0x10+ 通道掩码
通道头：0xAA+ 通道索引
参数：peak_freq,peak_magnitude,rms_value（3 个 float）
频谱数据：513 个 float（0~Nyquist 频率的幅度谱）

3.3.2矢量结果发送

VectorProcessing_SendResult()发送：

帧头：命令字0x20
DOA 数据：azimuth_deg,elevation_deg,confidence,target_freq
互谱峰值频率：Ix/Iy/Iz.peak_freq
互谱数据：128 个 float（Ix.cross_spectrum的前 128 点）

4.总结

目前测试下来，这套代码是可以跑通的，且我利用信号发射器模拟信号发射输出，得到的结果也是符合预期的，今天太晚了先到这里，明天多用几种方案再来验证一下我的这套系统

还有一些在实现期间遇到的问题，等我有时间整理一下下......

累！

基于矢量水听器的潜标探测系统——信号处理部分