# How can I reverse engineer the parameters of this butterworth filter?

I am trying to develop an application which will need a bandpass filter. The sampling rate is 30Hz, and the frequency range I'd like the filter to preserve is: 0.5Hz-4Hz.

I've used mkfilter to design my own filter but it is consistently being outperformed by this other filter I found in a sample project on github. The author didn't document the parameters he used and cannot be reached.

Would it be possible to reverse engineer the mkfilter parameters which produced this filter? (see swift code below)

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]
$$$$


The whole thing appears to be direct implementation of the difference equation, which is certainly not the best way to do this. Cascaded (or parallel) second order sections are much preferable in terms of stability and numerical noise.

Plotting the whole thing, we can see that the corner frequencies are at 1.5Hz and 9 Hz. It is indeed a Butterworth bandpass but with different cutoff frequencies.

Update based on the question in the comments

A linear time invariant filter can be described by it's difference equation

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

The transfer function can be obtained by the Z-transform as

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

The code implements the difference equation directly,i.e.

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

Looking at the code we get b = [1 0 -5 0 10 0 -10 0 5 0 -1] and a = [1 -3.2101976096 4.3832594900 -3.9628972812 3.3820859632 -2.4692165842 1.1743136181 -0.3992561394 0.1476158522 -0.0357796363 0]

We can just pop this into the equation for the transfer function and evaluate at any normalized frequency $$z = e^{j\omega}$$

In Matlab this would look like.

% 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]);
`
• Thanks, that totally answers the question, but allow me ask two more. 1. How did you figure that out (what tools did you use etc...?) 2. How would I implement it using cascaded or second order sections. Commented Feb 6, 2022 at 0:52
• I added the Matlab code to graph the filter. There is plenty of literature and existing code on second order sections. Every filter design class would cover this as well. Commented Feb 6, 2022 at 14:02
• Thanks so much. very helpful. Commented Feb 6, 2022 at 15:42