在 Matlab 中设计 Butterworth 滤波器并获得滤波器 [ab] 系数作为在线 Verilog HDL 代码生成器的整数

信息处理 matlab
2022-01-06 04:20:00

我使用 Matlab 设计了一个非常简单的低通巴特沃斯滤波器。下面的代码片段演示了我所做的。

fs = 2.1e6;
flow = 44 * 1000;
fNorm =  flow / (fs / 2);
[b,a] = butter(10, fNorm, 'low');

在[b,a]中存储了滤波器系数。我想将 [b,a] 获取为整数，以便我可以使用在线 HDL 代码生成器在 Verilog 中生成代码。

Matlab [b,a] 值似乎太小而无法与在线代码生成器一起使用（服务器端 Perl 脚本拒绝生成带有系数的代码），我想知道是否可以获得 [b, a] 以可用作正确输入的形式。

我在 Matlab 中得到的 a 系数是：

1.0000
-9.1585
37.7780
-92.4225
148.5066
-163.7596
125.5009
-66.0030
22.7969
-4.6694
0.4307

我在 Matlab 中得到的 b 系数是：

1.0167e-012
1.0167e-011
4.5752e-011
1.2201e-010
2.1351e-010
2.5621e-010
2.1351e-010
1.2201e-010
4.5752e-011
1.0167e-011
1.0167e-012

使用在线生成器，我想设计一个具有 12 位位宽和 I 或 II 滤波器形式的滤波器。我不知道上面链接中的“小数位”是什么意思。

使用上面列出的 [b,a] 系数运行代码生成器 (http://www.spiral.net/hardware/filter.html)，小数位设置为 20，位宽为 12，我收到以下运行错误：

Integer A constants: 1048576 -9603383 39613104 -96912015 155720456 -171714386 131597231 -69209161 23904282 -4896220 451621
Integer B constants: 0 0 0 0 0 0 0 0 0 0 0

Error: constants wider than 26 bits are not allowed, offending constant = -69209161, effective bitwidth = 7 mantissa + 20 fractional = 27 total.

An error has occurred - please revise the input parameters.

我该如何更改我的设计，以免发生此错误？

更新： 使用 Matlab 生成 6 阶巴特沃斯滤波器，我得到以下系数：

为一个：

1.0000
-5.4914
12.5848
-15.4051
10.6225
-3.9118
0.6010

对于 b：

0.0064e-005
0.0382e-005
0.0954e-005
0.1272e-005
0.0954e-005
0.0382e-005
0.0064e-005

运行在线代码生成器（http://www.spiral.net/hardware/filter.html），我现在收到以下错误（小数位为 8，位宽为 20）：

./iirGen.pl -A 256  '-1405' '3221' '-3943' '2719' '-1001' '153' -B  '0' '0' '0' '0' '0' '0' '0' -moduleName acm_filter -fractionalBits 8 -bitWidth 20 -inData inData  -inReg   -outReg  -outData outData -clk clk -reset reset -reset_edge negedge -filterForm 1  -debug  -outFile ../outputs/filter_1330617505.v 2>&1 
At least 1 non-zero-valued constant is required.  Please check the inputs and try again.

也许 b 系数太小，或者代码生成器 (http://www.spiral.net/hardware/filter.html) 想要另一种格式的 [b,a]？

更新：

也许我需要做的是将 [b,a] 系数按小数位数缩放，以获得整数形式的系数。

a .* 2^12
b .* 2^12

但是，我仍然认为 b 系数非常小。我在这里做错了什么？

也许另一种类型的过滤器（或过滤器设计方法）会更合适？有人可以提出建议吗？

更新： 正如 Jason R 和 Christopher Felton 在下面的评论中所建议的，SOS 过滤器会更合适。我现在已经编写了一些 Matlab 代码来获取 SOS 过滤器。

fs = 2.1e6;
flow = 44 * 1000;      
fNorm =  flow / (fs / 2);
[A,B,C,D] = butter(10, fNorm, 'low');
[sos,g] = ss2sos(A,B,C,D);

我得到的 SOS 矩阵是：

1.0000    3.4724    3.1253    1.0000   -1.7551    0.7705
1.0000    2.5057    1.9919    1.0000   -1.7751    0.7906
1.0000    1.6873    1.0267    1.0000   -1.8143    0.8301
1.0000    1.2550    0.5137    1.0000   -1.8712    0.8875
1.0000    1.0795    0.3046    1.0000   -1.9428    0.9598

是否仍然可以使用 Verilog 代码生成工具 (http://www.spiral.net/hardware/filter.html) 来实现这个 SOS 过滤器，或者我应该简单地手动编写 Verilog？有没有好的参考资料？

我想知道在这种情况下使用 FIR 滤波器是否会更好。

此外：递归 IIR 滤波器可以通过将系数表示为分数来使用整数数学来实现。（有关详细信息，请参阅 Smith 出色的 DSP 信号处理书：http ://www.dspguide.com/ch19/5.htm ）

以下 Matlab 程序使用 Matlab rat() 函数将巴特沃斯滤波器系数转换为小数部分。然后如评论中所述，二阶部分可用于数字实现过滤器（http://en.wikipedia.org/wiki/Digital_biquad_filter）。

% variables
% variables
fs = 2.1e6;                     % sampling frequency           
flow = 44 * 1000;               % lowpass filter


% pre-calculations
fNorm =  flow / (fs / 2);       % normalized freq for lowpass filter

% uncomment this to look at the coefficients in fvtool
% compute [b,a] coefficients
% [b,a] = butter(7, fNorm, 'low');
% fvtool(b,a)  

% compute SOS coefficients (7th order filter)
[z,p,k] = butter(7, fNorm, 'low');

% NOTE that we might have to scale things to make sure
% that everything works out well (see zp2sos help for 'up' and 'inf' options)
sos = zp2sos(z,p,k, 'up', 'inf'); 
[n,d] = rat(sos); 
sos_check = n ./ d;  % this should be the same as SOS matrix

% by here, n is the numerator and d is the denominator coefficients
% as an example, write the the coefficients into a C code header file
% for prototyping the implementation

 % write the numerator and denominator matices into a file
[rownum, colnum] = size(n);  % d should be the same
sections = rownum;           % the number of sections is the same as the number of rows
fid = fopen('IIR_coeff.h', 'w');

fprintf(fid, '#ifndef IIR_COEFF_H\n');
fprintf(fid, '#define IIR_COEFF_H\n\n\n');
for i = 1:rownum
   for j = 1:colnum

       if(j <= 3)  % b coefficients
            bn = ['b' num2str(j-1) num2str(i) 'n' ' = ' num2str(n(i,j))];
            bd = ['b' num2str(j-1) num2str(i) 'd' ' = ' num2str(d(i,j))];
            fprintf(fid, 'const int32_t %s;\n', bn);
            fprintf(fid, 'const int32_t %s;\n', bd);

       end
       if(j >= 5)  % a coefficients
            if(j == 5) 
                colstr = '1'; 
            end
            if(j == 6) 
                colstr = '2'; 
            end
            an = ['a' colstr num2str(i) 'n' ' = ' num2str(n(i,j))];
            ad = ['a' colstr num2str(i) 'd' ' = ' num2str(d(i,j))];
            fprintf(fid, 'const int32_t %s;\n', an);
            fprintf(fid, 'const int32_t %s;\n', ad);
       end
   end
end

% write the end of the file
fprintf(fid, '\n\n\n#endif');
fclose(fid);
3个回答

“小数位”是总线中您专用于表示数字的小数部分的位数（例如，3.75 中的 0.75）。

假设你有一个 4 位宽的数字总线，1001代表什么数字？如果将其视为正整数 (2^3 + 2^0 = 8 + 1 = 9)，它可能表示“9”。或者它可能意味着 -7 在二进制补码表示法中：(-2^3 + 2^0 = -8 + 1 = -7)。

那些带有一些分数的数字，即“实数”数字呢？实数可以在硬件中表示为“定点”或“浮点”。看起来那些过滤器生成器使用固定点。

回到我们的 4 位总线 ( 1001)。让我们引入一个二进制点，所以我们得到1.001. 这意味着现在使用点的 RHS 上的位来构建整数，并使用 LHS 上的位来构建分数。数字总线表示的数字设置1.001为 1.125 ( 1*2^0 + 0*2^-1 + 0*2^-2 + 1*2^-3 = 1 + 0.125 = 1.125)。在这种情况下，在总线的 4 位中，我们使用其中的 3 位来表示数字的小数部分。或者，我们有 3 个小数位。

因此，如果您有一个像上面那样的实数列表，您现在必须决定要表示它们的小数位数。这是权衡：您使用的小数位越多，您可以表示您想要的数字越接近，但您的电路需要更大。更重要的是，您使用的小数位数越少，滤波器的实际频率响应将偏离您一开始设计的频率响应越远！

更糟糕的是，您正在寻求构建一个无限脉冲响应 (IIR) 滤波器。如果您没有足够的小数位和整数位，这些实际上可能会变得不稳定！

所以Marty很好地处理了比特问题。关于滤波器本身，我认为您可能会从 matlab 收到有关缩放系数不佳的警告或投诉？当我绘制过滤器时，来自 scipy 而不是 matlab，但它可能非常相似。

回复

在通带处下降 100 dB！因此，您可能需要确保您需要一个更小的阶过滤器，这将有助于您的实施。当我使用 6 阶滤波器时，我不再收到有关不良系数的抱怨。也许尝试减少订单，看看它是否仍然满足您的要求。

正如所讨论的，最好使用部分的总和，即将高阶滤波器分解为级联的二阶滤波器。更新后的问题有 SOS 矩阵。使用此代码和此处的示例，Python 对象可用于生成各个部分。

在matlab中

save SOS

在 Python 中

import shutil
import numpy
from scipy.io import loadmat
from siir import SIIR

matfile = loadmat('SOS.mat')  
SOS = matfile['SOS']
b = numpy.zeros((3,3))
a = numpy.zeros((3,3))
section = [None for ii in range(3)]
for ii in xrange(3):
    b[ii] = SOS[ii,0:3]
    a[ii] = SOS[ii,3:6]

    section[ii] = SIIR(b=b[ii], a=a[ii], W=(24,0))
    section[ii].Convert()  # Create the Verilog for the section
    shutil.copyfile('siir_hdl.v', 'iir_sos_section%d.v'%(ii))

可以在此处找到有关定点的其他信息

其它你可能感兴趣的问题
上一篇谁能解释一下为什么 IIR 滤波器没有线性相位？ 下一篇特征检测器和特征描述符有什么区别？