基于FIR的高效插值器

在超声相控接收系统中，需要对采样的超声信号做超分辨率处理。考虑超声回波信号具有明显的窄带特性，就可以使用DSP课堂中的插0滤波法，选FIR的原因是可具有线性相移并且更加稳定。

FIR实现内插的数学原理

N阶FIR滤波器的时域输入输出关系为：

$$
y[n] = \sum^{N-1}_{k=0} h[k]x[n-k]
$$

式中$y$是输出，$x$是输入，$h$是滤波器的系数。本文中不对FIR的DSP知识进行回顾，直接使用Matlab工具设计滤波器。但是有必要分析一下，为什么可以使用补0+FIR方案实现超分辨率，以及FIR滤波器的参数如何确定。

首先，我们对原信号$x[n]$做有限长离散傅里叶变换，得到频谱$X$：

$$
X[k] = DFT(x[n])=\sum_{n=0}^{N-1} x[n]\cdot W_N^{kn}, 0\le k \le N-1
$$

同时逆变换有：

$$
x[n] = IDFT(X[k]) = \frac 1N \sum^{N-1}_{k=0} X[k]W_N^{-kn}, 0\le n\le N-1
$$

式中，$W_N=e^{-j2\pi /N}$是离散的旋转因子。如果在$x[i-1]$和$x[i]$之间插入$L-1$个0，那么插值后原本长度为$N$的序列$x[n]$变为了长度$(N-1)L+1$的序列$x_{ip}[n]$，注意$x_{ip}[jL]=x[j]$，做DFT有：

$$
X_{ip}[k]=DFT(x_{ip}[n]) = \sum_{n=0}^{(N-1)L} x_{ip}[n]\cdot W_{NL}^{kn}, 0\le k \le (N-1)L
$$

因为除了$n=jL, 0\le j \le N-1$之外其他都是0，因此有：

$$
X_{ip}[k] = \sum_{j=0}^{N-1} x_{ip} [jL] \cdot W_{NL}^{kjL}, 0 \le k \le (N-1)L
$$

$$
=\sum^{N-1}_{j=0} x[j]\cdot W_N^{kj}，0\le k \le (N-1)L
\\
=\begin{cases}
X[k], & 0 \le k \le N-1 \\
X[k-N], & N \le k \le 2N-1 \\
\vdots
\end{cases}
$$

因此我们可以发现原频谱$X$和内插0后频谱$X_{ip}$发生了对应内插倍数的频谱复制。从信号与系统的角度，只要我们的频谱没有发生混叠，那就一定可以恢复出原信号，这里可以看出频谱复制后是绝对不可能发生频谱重叠，但是如果原本信号的频谱两端并不为0，那么因为滤波器的不理想性也会导致还原信号错误。

为了演示，笔者使用高斯幅度调制正弦信号的超声脉冲回波信号作为原信号$x[n]$，进行插0、FIR操作，并在每一步做FFT，如下图读者就可以很简单地理解到具体原理。

对应的matlab分析代码：

figure()

%% 生成脉冲数据
F0 = 2e6;  % 谐振频率
bw_ratio = 0.9; % 带宽分数

Fs = 10e6; % 采样频率
T = 6*(1/F0);  % 信号采样时间范围，6个周期

% 创建时间向量
t = -T:1/Fs:T;
t_length = size(t);
t_length = t_length(2);

% 产生高斯脉冲
[yi, ~, ye] = gauspuls(t, F0, bw_ratio);


% FFT分析
[y_power,y_freq] = myEventFFT(yi,Fs/1e6);


% 绘制高斯脉冲
subplot(2,3,1)
plot(t+T, yi, t+T, ye);
title('Gaussian Modulated Ultrasonic Signal');
xlabel('Time (s)');
ylabel('x[n] Amplitude');

% 绘制频谱
subplot(2,3,4)
plot(y_freq, y_power, "LineWidth",2);
title('FFT(x[n])');
xlabel("f /MHz");
ylabel("|P(f)|" );

%% 插0

inter_mult = 10; % 数据插值倍率

y_after_inter = zeros(1,inter_mult*t_length);

for id = 1: t_length
    y_after_inter(id*inter_mult) = yi(id);
end

t_after_inter = (1:length(y_after_inter))/length(y_after_inter)*2*T;

% FFT分析
[y_inter_power,y_inter_freq] = myEventFFT(y_after_inter,Fs/1e6*inter_mult);

% 绘制插值后的图像和频谱
subplot(2,3,2)
plot(t_after_inter,y_after_inter)
title('after interpolate');
xlabel('Time (s)');
ylabel('x_{ip}[n] Amplitude');

subplot(2,3,5)
plot(y_inter_freq, y_inter_power, "LineWidth",2);
title('FFT(x_{ip}[n])');
xlabel("f /MHz");
ylabel("|P(f)|" );

%% FIR带通滤波


% 设置滤波器的频率点
% 归一化： 实际频率/(采样率/2)
f_fir_lp = [0 (F0*bw_ratio)/(Fs*inter_mult/2) (F0*bw_ratio)/(Fs*inter_mult/2) ...
            (F0*(2-bw_ratio))/(Fs*inter_mult/2) (F0*(2-bw_ratio))/(Fs*inter_mult/2) 1];
m_fir_lp = [0 0 1 1 0 0];

b_fir_lp = fir2(80,f_fir_lp,m_fir_lp);

% 换带通方案更好
f_fir_bp = [(F0*bw_ratio)/(Fs*inter_mult/2) (F0*(2-bw_ratio))/(Fs*inter_mult/2)];
b_fir_bp = fir1(80,f_fir_bp);

%显示滤波器特性
% freqz(b_fir_bp,1)

% 滤波
y_after_filter = filter(b_fir_bp,1,y_after_inter);

% FFT分析
[y_fil_power,y_fil_freq] = myEventFFT(y_after_filter,Fs/1e6*inter_mult);

% 绘制滤波后的图像和频谱
subplot(2,3,3)
plot(t_after_inter,y_after_filter)
title('after filter');
xlabel('Time (s)');
ylabel('x_{f}[n] Amplitude');


subplot(2,3,6)
plot(y_fil_freq, y_fil_power, "LineWidth",2);
title('FFT(x_{f}[n])');
xlabel("f /MHz");
ylabel("|P(f)|" );

%% 定义辅助函数
function [vec_power,vec_freq] = myEventFFT(vec_smpData, F_s)
%求采样数据的FFT变换，绘图，返回格式 [单边功率, 频率]
%   采样率共用导入数据时候的全局参数 F_s
% 
% global F_s;

% 数据长度
L_smpData = length(vec_smpData);

% 求频谱
vec_freq = F_s/L_smpData*(0:(L_smpData/2));  % 频率轴

P2 = abs( fft(vec_smpData) / L_smpData);  %FFT
P1 = P2(1 : L_smpData/2+1); % 取单边
P1(2:end-1) = 2*P1(2:end-1);

vec_power = P1;

end

内插器RTL框图

用于内插的FIR滤波器都可以做运算优化，因为插0后有很多运算都是0乘操作。比如上例中设计了N=80阶滤波器，数据倍率为L=10，若此时推入滤波器的$x_{ip}[n]$为有效采样数据，则：

$$
y[n] = \sum_{k=0}^{79} h[k]x_{ip}[n-k] \\
= \sum_{k}^{0,10,20,\cdots,70} h[k]x_{ip}[n-k] \\
=\sum_{i=0}^{7} h[i\cdot 10] x_{ip}[n-i \cdot 10]
$$

也就是可以优化为只有N/L=8次乘法操作。类似的，考虑前一时钟推入的为有效采样数据，此时时钟推入的为插值的0，那么：

$$
y[n]= \sum_{i=0}^{7} h[i\cdot 10-1]x_{ip}[n-i \cdot 10-1]
$$

还是只进行8次乘法操作，并且会发现在L=10次运算中，有效的$x_{ip}$是同一组数据。因此，作如下图设计，首先每隔L次时钟，滑动更新smpData寄存器数据；而系数全部存在RAM中，总共有N/L=8个RAM读端口，RAM地址由一个计数器提供，从而实现滑动系数；计数器给每个RAM的地址依次相差L。

根据上述优化分析，可以设计如下插值器的框图。

具体实现

使用DSP48E2原语

Ultrascale系列的DSP单元，如下图支持加法、乘法、逻辑判断。

直接基于DSP48E2原语进行设计，能够获得很多硬件层面的加速优势和功耗优势，比如：

• 实现更低能耗和高性能的流水线设计
• 使用BRAM或者UltraRAM来存储滤波器参数，从而优化CLB的使用
• 如果不用乘法可以USE_MUL=NONE来减小功耗
• 当多个DSP间有数据传递的关系，可以使用专用的级联（Cascade）资源从而减小对通用布线资源的占用。专用资源也能够获得更好的运算性能和更低的功耗
• 可以使用时间复用的方法，在低速应用时

参考文档

• PG323 DSP Macro LogiCORE IP Product Guide
• UG579 UltraScale Architecture DSP Slice User Guide

Donate