多相抽取器实现及matlab示例

1. 概要

        多相滤波器并不是指具有某种特性的时域或者频域特性的滤波器，它指的是一种滤波器实现的结构，常用于多速率数字信号处理（multi-rate digital signal processing）系统之中。

      以下采样处理为例。假设要对输入数据序列进行M倍下采样（downsampling or decimation），通常来说，为了避免混叠，在下采样之前需要先进行抗混叠低通滤波，然后再对滤波结果进行下采样。很显然，滤波输出中每M个数据有(M-1)个数据要被扔掉，这就意味着针对这(M-1)个数据的滤波处理完全被浪费了。多相滤波器就是用于解决这一计算浪费问题的一种有效的滤波器实现的结构。

        多相滤波器结构同样可以用于M倍上采样处理中，在上采样处理中，上采样之后需要进行滤波器处理以滤除镜像。常规的做法是先进行插零上采样（为什么是插零上采样而不是保持上采样呢？这个问题将在另外的文章中进行解释），然后再进行滤波。采用多相滤波器结构实现滤波器的话，同样可以将滤波和采样处理的顺序颠倒以大幅度降低所需要的运算量。

        本文主要说明多相抽取器（polyphase decimator，顾名思义，即以多相滤波器的方式实现下采样或者说抽取以及抗混叠滤波）的设计和实现。多相插值器的设计和实现将另外说明。

 

2. 诺贝尔恒等变换

        为什么可以进行上/下采样处理和滤波处理的顺序的交换？这个涉及到多速率信号处理中的基本的诺贝尔恒等变换（Nobel Identities，注意，不是高贵的身份！），如下图所示[1]。诺贝尔恒等变换所说明的就是在对滤波器进行适当地变形分解后，可以交换上/下采样处理和滤波处理的顺序而保持处理结果完全一致，此处暂不给出证明（以后补上或者另文说明），有兴趣的读者可以自行查阅相关文献。

图1 诺贝尔恒等变换（Nobel Identities）

 

3. 抗混叠滤波器的多相分解

        考虑输入信号为 , ,滤波器冲激响应为 ,则常规的先进行抗混叠滤波然后进行下采样的计算过程可以表示如下：

       

         经过适当的数学变换处理后可以得到：

        详细的（手撕）推导过程如下图所示： 

 

4. 多相抽取器实现结构

        根据上一节所得的数学表达式，可以很直观地映射得到多相滤波器的实现结构。

 

        经过分解以后，M个支路的子滤波器长度之和与分解之前的滤波器长度保持不变。但是这个M个滤波器是工作在下采样后的输出数据速率上的。考虑两种实现方式的处理量的比较，为了简便起见，这里仅考虑所需要的乘法运算。考虑原滤波器长度为L，按照直接形式（Direct Form）实现的话，滤波器每输出一个数据需要L次乘法。进一步，每M个输出中要扔掉（M-1）个数据只保留一个数据输出，则按照最终输出数据来算的话对应每个输出数据需要（L*M）个乘法运算。而以多相滤波器结构实现的话，对应每个输出数据只需要L个乘法运算。两者的运算量之比为(M:1)，M越大则运算量削减程度越大。即便考虑原滤波器的系数对称特性以对称折叠的形式实现也需要（L*M/2）个乘法运算，运算量之比也是(M/2:1)。 

 

5. Matlab示例

 

        在以下matlab示例程序中以两种方式实现下采样处理，并对比验证了两种实现方式所得到的处理结果完全一致。代码中也给出了多相滤波器结构实现中的滤波器的多相分解以及对应的输入数据序列的分解方式。

%2021-07-24 chenxy
%FIR decimator illustration. M = 4
clear; close all; clc

% Create input signal and filter
x = 1:128;
h = 1:1:48; % If the number of h is not integer times of M, should pad with zero.

% %%%%%% Direct Form (Inefficient): filtering followed by down-sampling %%%%%%
y =filter(h,1,x);  % Compute filter output first
y_ds = y(1:4:end); % 4:1 down-sampling

% %%%%%% Polyphase Form (Efficient):  %%%%%%
% Polyphase decomposition into 4 sub-filters 
sub_filt0 = h(1:4:end);
sub_filt1 = h(4:4:end);
sub_filt2 = h(3:4:end);
sub_filt3 = h(2:4:end);
% Generate input signals for 4 sub-filters
x0 = [  x(1:4:end)];
x1 = [0 x(2:4:end-4)];
x2 = [0 x(3:4:end-4)];
x3 = [0 x(4:4:end-4)];
% filter each polyphase component and add together
y_poly_ds = filter(sub_filt0,1,x0)+filter(sub_filt1,1,x1)+filter(sub_filt2,1,x2)+filter(sub_filt3,1,x3);

plot(y_ds - y_poly_ds); title('The difference between two implementation');

% assert(isequal(y_ds,y_poly_ds)); % May fail for floating-point signal processing
assert(max(abs(y_ds,y_poly_ds)) < 1e-10);

 

 

Reference:

[1] Multirate Noble Identities