信息处理 - 如何对这个巴特沃斯滤波器的参数进行逆向工程？ - 吾爱随笔录

如何对这个巴特沃斯滤波器的参数进行逆向工程？

信息处理数字滤波器巴特沃思逆向工程

2022-02-23 08:27:46

我正在尝试开发一个需要带通滤波器的应用程序。采样率为 30Hz，我希望滤波器保留的频率范围是：0.5Hz-4Hz。

我使用 mkfilter 来设计我自己的过滤器，但它始终被我在 github 上的示例项目中找到的其他过滤器所超越。作者没有记录他使用的参数并且无法访问。

是否可以对产生此过滤器的 mkfilter 参数进行逆向工程？（见下面的快速代码）

xv[0] = xv[1]; xv[1] = xv[2];
xv[2] = xv[3]; xv[3] = xv[4];
xv[4] = xv[5]; xv[5] = xv[6];
xv[6] = xv[7]; xv[7] = xv[8];
xv[8] = xv[9]; xv[9] = xv[10];
xv[10] = value/gain

yv[0] = yv[1]; yv[1] = yv[2];
yv[2] = yv[3]; yv[3] = yv[4];
yv[4] = yv[5]; yv[5] = yv[6];
yv[6] = yv[7]; yv[7] = yv[8];
yv[8] = yv[9]; yv[9] = yv[10];

yv[10] = (xv[10] - xv[0]) + 5 * (xv[2] - xv[8]) +
         10 * (xv[6] - xv[4]) + (-0.0000000000 * yv[0]) +
         (0.0357796363 * yv[1]) + (-0.1476158522 * yv[2]) + 
         (0.3992561394 * yv[3]) + (-1.1743136181 * yv[4]) +
         (2.4692165842 * yv[5]) + (-3.3820859632 * yv[6]) +
         (3.9628972812 * yv[7]) + (-4.3832594900 * yv[8]) +
         (3.2101976096 * yv[9])

return yv[10]
```

1个回答

整个事情似乎是差分方程的直接实现，这当然不是最好的方法。就稳定性和数值噪声而言，级联（或并行）二阶部分更可取。

绘制整个图，我们可以看到转角频率为 1.5Hz 和 9Hz。它确实是巴特沃斯带通，但具有不同的截止频率。

根据评论中的问题更新

线性时不变滤波器可以用它的差分方程来描述

\sum_{k} a_{k} y [n - k] = \sum_{k} b_{k} x [n - k]

$\sum_k a_ky[n-k] = \sum_k b_kx[n-k]$

传递函数可以通过 Z 变换获得为

H (z) = \frac{\sum_{k} b_{k} z^{- k}}{\sum_{k} a_{k} z^{- k}}

$H(z) = \frac{\sum_k b_k z^{-k}}{\sum_k a_k z^{-k}}$

代码直接实现差分方程，即

y [n] = \sum_{k = 0}^{N - 1} b_{k} x [n - k] - \sum_{k = 1}^{N - 1} a_{k} y [n - k]

$y[n] = \sum_{k = 0}^{N-1} b_k x[n-k] - \sum_{k = 1}^{N-1} a_k y[n-k]$

查看我们得到的b = [1 0 -5 0 10 0 -10 0 5 0 -1]代码a = [1 -3.2101976096 4.3832594900 -3.9628972812 3.3820859632 -2.4692165842 1.1743136181 -0.3992561394 0.1476158522 -0.0357796363 0]

我们可以将其放入传递函数的方程中，并以任何归一化频率进行评估 $z = e^{j\omega}$

在 Matlab 中，这看起来像。

% filter coefficients
b = [1 0 -5 0 10 0 -10 0 5 0 -1];
a = [1 -3.2101976096  4.3832594900  -3.9628972812 3.3820859632 -2.4692165842 1.1743136181 -0.3992561394 0.1476158522 -0.0357796363  0];
% make a log frequency vector
fs = 30; % sampling rate
fr = logspace(log10(.3),log10(15),1000);
% calculate frequency response
H = freqz(b,a,fr,fs);
% Convert to enegery and normalize
E = H.*conj(H); E = E./max(E);
% and plot it
semilogx(fr,10*log10(E),'Linewidth',2);
% pretty up the graphs
grid on;
xlabel('Frequency in Hz');
ylabel('Level in dB');
set(gca,'ylim',[-10 1]);

其它你可能感兴趣的问题

上一篇如何知道滤波器是低通还是高通？下一篇如何在matplotlib中将正弦波的FFT居中？