探讨滤波器优化算法RAG-N和BHM算法的详细解读: 基于 MATLAB 和 Simulink 的递归算法(RAG-N OR BHM) FIR 滤波器量化和低成本结构优化设计-含MATLAB代码

探讨滤波器优化算法RAG-N和BHM算法的详细解读：Recursive Algorithm FIR Filter Quantization and Low-Cost Structure Optimization Design Based on MATLAB and Simulink（基于 MATLAB 和 Simulink 的递归算法(RAG-N OR BHM) FIR 滤波器量化和低成本结构优化设计）

IEEE引用：

Z. Shang, “Recursive Algorithm FIR Filter Quantization and Low-Cost Structure Optimization Design Based on MATLAB and Simulink,” 2023 IEEE 16th International Symposium on Embedded Multicore/Many-core Systems-on-Chip (MCSoC), Singapore, 2023, pp. 530-535, doi: 10.1109/MCSoC60832.2023.00084.

Bib TeX: ``

@INPROCEEDINGS{10387811,
  author={Shang, Zixia},
  booktitle={2023 IEEE 16th International Symposium on Embedded Multicore/Many-core Systems-on-Chip (MCSoC)}, 
  title={Recursive Algorithm FIR Filter Quantization and Low-Cost Structure Optimization Design Based on MATLAB and Simulink}, 
  year={2023},
  volume={},
  number={},
  pages={530-535},
  doi={10.1109/MCSoC60832.2023.00084}}

摘要

在数字滤波领域，有限脉冲响应（FIR）滤波器因其稳定的结构和线性相位特性而备受青睐。然而，由于传统的直接型滤波器运算速度较慢，且需要大量的硬件资源，因此许多研究人员探索使用递归方法来减少电路中加法器的使用，从而降低数字滤波器的硬件成本。其中，RAG-N 算法和 BHM 算法是该领域较为著名的算法。本文研究的重点是有限脉冲响应滤波器系数的量化性能以及递归算法的系数优化分析。本文首先设计并论证了经过 8 位、10 位和 12 位量化后采用递归算法优化的有限脉冲响应滤波器的结构，并在此基础上利用 Simulink 建立了相应的仿真模型。最终结果表明，虽然 RAG-N 和 BHM 等启发式优化算法可以显著减少硬件资源的使用，但随着量化位数的增加，这些算法的计算时间可能会增加，这可能会由于逻辑深度的增加和路径延迟的增加而对电路设计产生影响。因此，未来的研究应更多地关注这些问题，以提高 FIR 滤波器的性能和效率。Keywords- FIR filter; quantization; recursive algorithms; RAG-N; BHM; Matlab; Simulink

简介

随着数字通信技术的飞速发展，高质量信号处理的需求对滤波器的性能和资源消耗提出了更高的要求。在数字信号处理中，有限冲激响应（FIR）滤波器因其结构稳定简单、具有线性相位特性等优点，在实时信号处理中得到了广泛应用。
目前，广泛使用的 FIR 滤波器结构主要有直接结构、转置结构和分布式算法结构。直接结构将滤波运算转化为积的累加；转置结构在直接结构的基础上翻转加法方向，使输入输出方向一致，从而获得更高的工作频率限制。
FIR 滤波器的阶数和量化位宽是决定数字滤波器硬件资源消耗的两个重要因素。通常，阶数越高，量化位宽越大，硬件资源消耗就越大。由于 FIR 滤波器的参数是固定的，因此引入了固定乘法器优化算法来优化处理过程，即通过将乘法分解为加法和位移操作来降低计算复杂度。
在 FIR 滤波器硬件中，常量乘法通常采用少乘法器技术来实现，这种技术主要使用位移寄存器和加法器来代替乘法器。在滤波器架构中，加法器可进一步分为乘法器块（MB）中的乘法器块加法器（MBA）和延迟单元中的结构加法器（SA）。由于给定滤波器阶数的延迟单元和结构加法器数量相对固定，因此在现有的数字滤波器固定乘法器优化算法中，MB 数量一直是计算滤波器硬件复杂度的重要指标。
在本研究中，我们首次构建了一个 Simulink 模型，以进一步验证递归算法在优化低成本结构方面的实用性，并充分考虑量化位宽的影响。

滤波器结构和优化原理（RAG-N和BHM的工作原理）

有限脉冲响应（FIR）滤波器由一条分接延迟线和一组加法器和乘法器构成，其中每个乘法器接收的操作数都是 FIR 滤波器的系数。FIR 滤波过程本质上是一个逐级信号延迟过程，在此过程中，来自每一级延迟的输入信号经过加权和求和后生成 FIR 滤波器的输出，其核心操作是乘法累加操作。其核心操作是乘法累加操作。

如图 1(a)所示，输入信号与滤波器的每个常数系数相乘后被送入延迟单元。这一操作通常被称为多常数乘法（MCM）问题，主要由移位寄存器和加法器网络实现。如图 1(b)所示，加法器可细分为用于延迟单元的结构加法器（SA）和用于常数乘法单元的乘法器块加法器（MBA）。
当滤波器阶数确定后，延迟单元和 SA 的数量就相对固定了（不过，如果某些系数为零，SA 的数量可能会减少）。因此，FIR 滤波器的实现复杂度主要取决于 MBA 的数量。换句话说，通过减少滤波器的 MBA 数量，可以有效降低滤波器的复杂度。

A. 乘法器子空间（MBA）中的传统二进制系数加法器

如图 2 所示，两个滤波器系数的二进制表示被设置为 10101。如果这两个系数独立实现，那么每个系数将需要 2 个乘法单元加法器 (MBA)，因此总共需要 4 个 MBA。然而，通过适当的递归优化算法，我们能够合并这些重复的滤波器系数，从而将加法器的总需求降至 2 个。这表明，适当使用递归优化算法（图形启发式优化算法）可以大大减少数字滤波器的硬件资源消耗[8]。

B. 递归优化算法（RAG-N 和 BHM）

以常数系数集 {1, 7, 16, 21, 33} 为例，为了对同一输入信号进行乘法运算，可以使用有向无环图（DAG）集中生成所有系数乘法运算[5]，如图 3 所示。图 3(b) 显示了递归算法图-N（简称 RAG-N）生成的结果。
递归优化算法的基本原理在于，已经生成的节点（或称为基本元素）可以用来生成尚未生成的系数。例如，只需添加一个加法器，就能从现有的 7 生成 21，而如果单独生成 21，则需要两个加法器（21=2⁴+2²+1）。因此，递归优化算法可以有效减少整个乘法模块所需的加法器总数。
以 RAG-N 为例，该算法包含两个部分：优化部分和启发式部分 [5]。在算法执行过程中，需要两个查找表：第一个查找表列出了每个系数单独实现所需的最小加法器数量（即单个系数的最优成本），第二个查找表则列出了实现每个系数最优成本的具体方法，这些方法可能有多种，例如，5 可以由 2 + 3 或 7 - 2 得到。

Matlab 实验设计(先量化滤波器系数再使用算法精简)

A. 滤波器参数设置

本研究设计的有限脉冲响应（FIR）滤波器具有以下特点：采样频率为 16 kHz，属于低通滤波器类型；通带纹波小于或等于 1 dB；在阻带区域的衰减至少为 30 dB；通带截止频率设为 3.4 kHz，阻带截止频率设为 4 kHz。这种设计使滤波器既能实现其功能，又能有效控制滤波过程中可能出现的误差和干扰。

B. 滤波器量化要求

在实际硬件实施中，处理器以有限精度（通常是固定位宽）表示数据。通过量化过程，理论上无限精度的数字被转换为有限精度的表示，从而可以在实际硬件中运行。这为后续优化和模拟滤波器系数提供了条件。
量化过程可分为四个步骤：

1. 找出滤波器系数的最大绝对值 M。
2. 用 M 对所有滤波器系数进行归一化处理，即用 M 除以所有系数。
3. 将所有滤波系数乘以 (2^(Q-1)-1)，并将处理后的系数四舍五入，形成整数系数。
4. 将整数滤波系数转换为二进制补码数据。
   本文将重点介绍量化过程及其对 8 位（8 位）、10 位（10 位）和 12 位（12 位）情况的分析。

Matlab 实验结果

A. 滤波器量化

图 4 和图 5 分别显示了不同量化位宽条件下滤波器的幅频响应图。


从图 4 中的归一化幅频图可以看出，滤波器系数的量化位宽越小，量化后的滤波性能就越差。相反，量化位宽越大，量化前后的滤波性能差异就越小。当量化位宽达到一定值（12 位，本例中为 12 位）时，滤波性能几乎与原始滤波器相同。
从图 5 的非规范化幅频图（（非规范化））可以看出，量化位宽越大，滤波器通带的增益就越大。例如，当量化位宽为 12 比特时，增益可达 72 dB。这一结果与我们之前的讨论是一致的。这是因为，根据上述滤波器量化方法，系数的量化过程相当于滤波器增益扩大了 (2^(Q-1)-1)/M 的系数。

B. 递归优化算法（RAG-N 和 BHM）的复杂性优化结果

鉴于 BHM 算法和 RAG-N 算法的经典性，本文将使用这两种算法作为递归优化算法的代表进行分析。

图 6 显示了从 4 位量化到 12 位量化的原始滤波器复杂度与经过 RAG-N 和 BHM 优化后的滤波器复杂度的对比图。从图 7 可以看出，滤波器的量化位宽越大，滤波器的复杂度就越高，优化效果就越差。这说明在选择优化策略时，需要根据实际情况权衡量化位宽与硬件资源消耗之间的关系。

SIMULINK仿真模拟验证

在模拟设计中，信号输入由两个数字信号处理器（DSP）和一个随机信号源组成。图 7 模拟了有限脉冲响应（FIR）滤波器的量化频率图，其中 "低通 "模块作为基准，是根据实验数据设计的。Myfilter "等模块表示使用直接型 FIR 滤波器结构的初始和量化滤波器。


图 8 验证了用递归算法优化的 FIR 滤波器结构的正确性，其系数量化为 8、10 和 12 比特。由于实验中的 RAG-N 算法和 BHM 算法的简化系数仅在排列方式上有所不同，因此本文仅以优化后的 RAG-N 仿真图作为递归算法优化的代表。图 9、图 10 和图 11 则分别设计了基于优化系数的滤波器结构图。



在实际应用中，需要根据具体要求和资源限制选择合适的优化策略和量化位宽，以获得最佳的滤波器性能和硬件资源使用效率。

A. 仿真结果

图 12、13、14 和 15 展示了 Simulink 仿真结果，与 MATLAB 实验结果一致，验证了滤波器结构设计的准确性。图 12 比较了直接型 FIR 滤波器的量化频谱。图 13、图 14 和图 15 则计算并比较了递归算法优化后的频谱图。




图 16 显示了各 Simulink 模块的运行时间。使用 RAGN 和 BHW 优化算法的 8 位量化滤波器结构速度最快，而 "低通 "模块速度最慢。值得注意的是，与直接结构相比，10 位和 12 位优化结构的速度并没有显著提高。这可能是由于优化滤波器结构的逻辑深度相对较高，原因是多次重复使用了常用术语，这也带来了明显的路径延迟。

这些结果表明，虽然优化算法可以减少硬件资源的使用，但也可能增加逻辑深度和路径延迟，从而影响运行速度。

结论

根据本文对递归优化算法Simulink的仿真，补充了现阶段证明递归优化算法优化优缺点的理论不足。递归优化算法（图启发式优化算法），如BHM、RAG-N、HCUB等，可以***大大节省硬件资源。不过，这些算法的计算时间会随着量化位数的增加而增加。此外，由于共同项的多重复用，这些算法可能会导致逻辑深度增加，并带来显著的路径延迟，所有这些都可能对电路设计产生不利影响。***因此，在未来的有限脉冲响应（FIR）滤波器优化研究中，我们需要深入解决这些问题，并对其进行优化，以进一步提高 FIR 滤波器的性能和效率。同时，在实际应用中，我们需要根据具体要求和资源限制，考虑选择合适的优化策略和量化位宽。

部分代码

IEEE下载链接：*https://ieee-dataport.org/documents/recursive-algorithm-fir-filter-quantization-and-low-cost-structure-optimization-design*

clc
clear all
close all
%% Filter Design
fpass = 3400;           % Pass band frequency (Hz)
fstop = 4000;           % stop band frequency (Hz)
ripple = 1;           % Passband ripple in dB 
atten = 30;          % Stopband ripple in dB
fs = 16000;        % Sampling frequency(Hz)
f = [fpass fstop];    % Cutoff frequencies(Hz)
a = [1 0];        % Desired amplitudes
dev = [(10^(ripple/20)-1)/(10^(ripple/20)+1) 10^(-atten/20)]; 
[n,fo,ao,wp] = firpmord(f,a,dev,fs);
b = firpm(n,fo,ao,wp);
% Quantization Filter Coefficients-Normalization
b2 = MYquantize(b, 2);
b4 = MYquantize(b, 4);
b6 = MYquantize(b, 6);
b8 = MYquantize(b, 8);
b10 = MYquantize(b, 10);
b12 = MYquantize(b, 12);
% Filter optimization coefficients-unnormalized
b01 = MYquantize1(b, 1);
b21 = MYquantize1(b, 2);
b41 = MYquantize1(b, 4);
b61 = MYquantize1(b, 6);
b81 = MYquantize1(b, 8);
b101 = MYquantize1(b, 10);
b121 = MYquantize1(b, 12);
% Parks-McClellan Bandpass Filter-Calculate the frequency response of the filter
[h,w] = freqz(b,1,1024);
[H,f] = freqz(b, 1, 1024, fs);
[H8,f8] = freqz(b8, 1, 1024, fs);
[H10,f10] = freqz(b10, 1, 1024, fs);
[H12,f12] = freqz(b12, 1, 1024, fs);
[H81,f81] = freqz(b81, 1, 1024, fs);
[H101,f101] = freqz(b101, 1, 1024, fs);
[H121,f121] = freqz(b121, 1, 1024, fs);


figure(1)
freqz(b,1,1024,fs)
figure(2)
hold on
grid on
plot([0 fpass fstop fs/2]/(fs/2),[1 1 0 0],w/pi,abs(h))
plot([0 fpass fstop fs/2]/(fs/2),[1 1 0 0],w/pi,abs(H))
plot([0 fpass fstop fs/2]/(fs/2),[1 1 0 0],w/pi,abs(H8))
plot([0 fpass fstop fs/2]/(fs/2),[1 1 0 0],w/pi,abs(H10))
plot([0 fpass fstop fs/2]/(fs/2),[1 1 0 0],w/pi,abs(H12))
% [0 fpass fstop fs/2]/(fs/2),[1 1 0 0],f12/pi,
legend('Ideal','firpm Design','8-bit', '10-bit', '12-bit')
xlabel 'Radian Frequency (\omega/\pi)', ylabel 'Magnitude'
figure(3)
hold on;
plot(f, 20*log10(abs(H)), 'b');
plot(f8, 20*log10(abs(H8)), 'r');
plot(f10, 20*log10(abs(H10)), 'g');
plot(f12, 20*log10(abs(H12)), 'k');
grid on
xline(fpass,'--b','Pass band frequency')
yline(-ripple/2,'--r','Pass band ripple')
yline(ripple/2,'--r')
xline(fstop,'--c','Stop band attenuation')
yline(-atten,'--c','Stop band attenuation') 
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
legend('Original', '8-bit', '10-bit', '12-bit');
title('Lowpass Filter Designed to Specifications')
figure(4)
hold on;
plot(f, 20*log10(abs(H)), 'b');
plot(f8, 20*log10(abs(H81)), 'r');
plot(f10, 20*log10(abs(H101)), 'g');
plot(f12, 20*log10(abs(H121)), 'k');
grid on
xline(fpass,'--b','Pass band frequency')
yline(-ripple/2,'--r','Pass band ripple')
yline(ripple/2,'--r')
xline(fstop,'--c','Stop band attenuation')
yline(-atten,'--c','Stop band attenuation') 
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
legend('Original', '8-bit', '10-bit', '12-bit');
title('Lowpass Filter Designed to Specifications')
%% Filter Performance
% Calculate the performance parameters of the filter
ripple_b = 20*log10(max(abs(H(f <= fpass)))) - 20*log10(min(abs(H(f <= fpass))));
atten_b = -20*log10(max(abs(H(f >= fstop))));
ripple_b8 = 20*log10(max(abs(H8(f8 <= fpass)))) - 20*log10(min(abs(H8(f8 <= fpass))));
atten_b8 = -20*log10(max(abs(H8(f8 >= fstop))));
ripple_b10 = 20*log10(max(abs(H10(f10 <= fpass)))) - 20*log10(min(abs(H10(f10 <= fpass))));
atten_b10 = -20*log10(max(abs(H10(f10 >= fstop))));
ripple_b12 = 20*log10(max(abs(H12(f12 <= fpass)))) - 20*log10(min(abs(H12(f12 <= fpass))));
atten_b12 = -20*log10(max(abs(H12(f12 >= fstop))));
% Print results
fprintf('Original filter: passband ripple = %f dB, stopband attenuation = %f dB\n', ripple_b, atten_b);
fprintf('8-bit quantized filter: passband ripple = %f dB, stopband attenuation = %f dB\n', ripple_b8, atten_b8);
fprintf('10-bit quantized filter: passband ripple = %f dB, stopband attenuation = %f dB\n', ripple_b10, atten_b10);
fprintf('12-bit quantized filter: passband ripple = %f dB, stopband attenuation = %f dB\n', ripple_b12, atten_b12);
% Calculate the order and coefficient complexity of the filter
% Order
order = length(b) - 1;
order8 = length(b8) - 1;
order10 = length(b10) - 1;
order12 = length(b12) - 1;
% Coefficient complexity- A multiplication operation is required for each non-zero coefficient. Complexity is defined as the number of nonzero elements in the filter coefficients.
complexity = sum(b ~= 0);
complexity8 = sum(b8 ~= 0);
complexity10 = sum(b10 ~= 0);
complexity12 = sum(b12 ~= 0);
% Filter Costs
filter = dfilt.dffirt(b);
filter2 = dfilt.dffirt(b2);
filter4 = dfilt.dffirt(b4);
filter6 = dfilt.dffirt(b6);
filter8 = dfilt.dffirt(b8);
filter10 = dfilt.dffirt(b10);
filter12 = dfilt.dffirt(b12);
% Calculating the complexity of filters
b_cost = cost(filter);
b2_cost = cost(filter2);
b4_cost = cost(filter4);
b6_cost = cost(filter6);
b8_cost = cost(filter8);
b10_cost = cost(filter10);
b12_cost = cost(filter12);
%Filter coefficients
% b
b_adds = b_cost.NAdd;
b_mults = b_cost.NMult;
b_states = b_cost.NStates;
% Number of multiplication operations per input sample每输入采样的乘法运算数
b_NMultInput= b_cost.MultPerInputSample;
% Number of addition operations per input sample每输入采样的加法运算数
b_NaddInput = b_cost.AddPerInputSample;
% 2
b2_adds = b2_cost.NAdd;
b2_mults = b2_cost.NMult;
b2_states = b2_cost.NStates;
% Number of multiplication operations per input sample每输入采样的乘法运算数
b2_NMultInput= b2_cost.MultPerInputSample;
% Number of addition operations per input sample每输入采样的加法运算数
b2_NaddInput = b2_cost.AddPerInputSample;
% 4
b4_adds = b4_cost.NAdd;
b4_mults = b4_cost.NMult;
b4_states = b4_cost.NStates;
% Number of multiplication operations per input sample每输入采样的乘法运算数
b4_NMultInput= b4_cost.MultPerInputSample;
% Number of addition operations per input sample每输入采样的加法运算数
b4_NaddInput = b4_cost.AddPerInputSample;
% 6
b6_adds = b6_cost.NAdd;
b6_mults = b6_cost.NMult;
b6_states = b6_cost.NStates;
% Number of multiplication operations per input sample每输入采样的乘法运算数
b6_NMultInput= b6_cost.MultPerInputSample;
% Number of addition operations per input sample每输入采样的加法运算数
b6_NaddInput = b6_cost.AddPerInputSample;
% 8
b8_adds = b8_cost.NAdd;
b8_mults = b8_cost.NMult;
b8_states = b8_cost.NStates;
% Number of multiplication operations per input sample每输入采样的乘法运算数
b8_NMultInput= b8_cost.MultPerInputSample;
% Number of addition operations per input sample每输入采样的加法运算数
b8_NaddInput = b8_cost.AddPerInputSample;
% 10
b10_adds = b10_cost.NAdd;
b10_mults = b10_cost.NMult;
b10_states = b10_cost.NStates;
% Number of multiplication operations per input sample每输入采样的乘法运算数
b10_NMultInput= b10_cost.MultPerInputSample;
% Number of addition operations per input sample每输入采样的加法运算数
b10_NaddInput = b10_cost.AddPerInputSample;
% 12
b12_adds = b12_cost.NAdd;
b12_mults = b12_cost.NMult;
b12_states = b12_cost.NStates;
% Number of multiplication operations per input sample每输入采样的乘法运算数
b12_NMultInput= b12_cost.MultPerInputSample;
% Number of addition operations per input sample每输入采样的加法运算数
b12_NaddInput = b12_cost.AddPerInputSample;

figure(5)
% 绘制条形图
subplot(2,1,1);
bar([order, order8, order10, order12]);
ylabel('Order');
set(gca, 'XTickLabel', {'Original', '8-bit', '10-bit', '12-bit'});
subplot(2,1,2);
bar([complexity, complexity8, complexity10, complexity12]);
ylabel('Complexity');
set(gca, 'XTickLabel', {'Original', '8-bit', '10-bit', '12-bit'});
figure(6)
bar([b_adds, b_mults,b_states,b_NMultInput,b_NaddInput;b8_adds, b8_mults,b8_states,b8_NMultInput,b8_NaddInput;b10_adds, b10_mults,b10_states,b10_NMultInput,b10_NaddInput; b12_adds, b12_mults,b12_states,b12_NMultInput,b12_NaddInput]);
set(gca, 'XTickLabel', {'Original', '8-bit', '10-bit', '12-bit'});
ylim([25 33]);
grid on
legend('Adds', 'Multiplications','states','multiplication operations per input sample','addition operations per input sample');
ylabel('Operation Count');