How does a low-pass filter programmatically work?

Question

I have been working on a simple low pass filter for < 100 Hz metering in my application. But so far, I am struggling with the theory behind it all. It's cool that I got it working, but I'd really enjoy it if I knew how/why it is working.

I found the following code:

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
    double PI = M_PI;
    double sqrt2 = sqrt(2);

    double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
    double QcWarp = tan(QcRaw); // Warp cutoff frequency

    double gain = 1 / ( 1 + sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) );

    by[2] = ( 1 - sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) ) * gain;
    by[1] = ( 2 - 2 * 2 / ( QcWarp * QcWarp ) ) * gain;
    by[0] = 1;

    ax[0] = 1 * gain;
    ax[1] = 2 * gain;
    ax[2] = 1 * gain;
}

To calculate the Coefficients. Then, in the audio samples, I 'low pass' them this way:

        xv[2] = xv[1];
        xv[1] = xv[0];

        xv[0] = pData[j];
        yv[2] = yv[1];
        yv[1] = yv[0];

        yv[0] = (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                   - by[1] * yv[0]
                   - by[2] * yv[1]);

        pData[j] = yv[0];

To get a lowpass design.

I am wondering a few things:

1. I receive the audio samples in a simple float* array. What is that float number? The only thing I see is a number, how is that a piece of sound?
2. The code is using past calculations (three of them) in the new calculation per sample. Does that mean that the first 2 data samples are not filtered correctly? (not that it would matter because it's just 2 samples, but just wondering)
3. Trying to learn everything, I found a couple of formulas for the Butterworth (2nd Pole) filter. How are those formulas reflected in this code? None of the formulas I found have these calculations that you can see in the 'getLPCoefficientsButterworth2Pole()' function.

Answer 1

1. The number float* array is a pointer to the array. It is a single number which contains the address of the first element of the array of float values.

2. Usually, the initial condition (i.e. the initial 'past' elements of x and y) are 0, but if their values are not equal to 0 it is no big problem either, because after a while the initial conditions have no effect on the output signal for any stable filter. And your filter is obviously stable.

3. A second order lowpass transfer function with Butterworth characteristic and cut-off frequency $\omega_a=2*\pi*f_a$ (with $f_a$ in Hertz) is given by

$$H(s) = \frac{\omega_a^2}{s^2+\sqrt{2}\omega_a s + \omega_a^2}\tag{1}$$ This is a standard result which you can easily find on the web. In order to get a discrete-time filter you can apply a so-called bilinear transformation:

$$s = 2f_s\frac{z-1}{z+1}\tag{2}$$

where $f_s$ is the sampling frequency. This is necessary because the frequency axis of the analog signals ($0\le f\le\infty$) needs to be mapped to the allowed frequency range for discrete-time signals ($0\le f\le f_s/2$). Since this transformation warps the frequencies, we need to compute the desired analog cut-off frequency from the given cut-off frequency $f_d$ in the discrete-time domain:

$$\omega_a = 2f_s \tan(\frac{\omega_d}{2})\text{ with }\omega_d = 2\pi f_d/f_s$$

If you insert (2) into (1) you get

$$H(z) = k\frac{z^2+b_1z+b_2}{z^2+a_1z+a_2}\tag{3}$$

with

$$k=\frac{\alpha^2}{1+\sqrt{2}\alpha+\alpha^2}\text{ (gain)}$$ $$a_1 = \frac{2(\alpha^2-1)}{1+\sqrt{2}\alpha+\alpha^2},\quad a_2=\frac{1-\sqrt{2}\alpha+\alpha^2}{1+\sqrt{2}\alpha+\alpha^2},\quad b_1=2,\quad b_2=1$$

where I've used $\alpha=\tan(\frac{\omega_d}{2})$. The transfer function (3) corresponds to the time domain equation

$$y(n) = kx(n) + kb_1x(n-1) + kb_2x(n-2) - a_1y(n-1) - a_2y(n-2)$$

As you can see, there's quite some resemblance with your code. However, there are some differences too and you should check the frequency response of your filter. I think that the above equations are correct, but it's up to you to check them and decide which version is the correct one.

Answer 2

You asked how a low pass filter works and mentioned that the filter uses past values of you're data. This is a non-technical discussion of what happens in a low pass filter.

The low pass filter takes differing views (shifted in time) of your signal, scales them and adds them together. You can imagine drawing your signal 3 times, one being current, the second being shifted by one sample time, the third being shifted by 2 sample times.

At low frequencies, all the views look very similar (shifting by a single sample barely changes where you are on the signal at any instant in time). In this case, the three versions will add together in a constructive (or at least non destructive) way, so the signal passes through the filter.

Now moving to higher frequencies, each shifted version of the signal becomes more distinct at any given instant (sample point) and in fact may even reverse in sign. At these higher frequencies, the three versions of your signal tend to cancel out (added destructively) so the signal becomes attenuated.

Different types of filters arrange for this kind of constructive/destructive interference to take place over appropriate frequency bands to create low-pass, band-pass or high-pass filters.

Answer 3

It depends on how you want it. For me I have implemented in C. So, What i did was that I generated filter coefficients in matlab using firls. Then I stored that filter coefficients in an array and then pass these filter coefficients to the filter function in C. In matlab it is quite easy to generated the filter coefficients, and then use the convolution operation or use a specific filter. But in C you have to implement the filter as well. For my part I have to do it for DSPs, so implemented in C and got the results in matlab using mex function. To calculate the coefficients use this command in matlab :

n=17                  %filter order
f=[0 0.167 0.333 1]   %Frequency band edges
a=[1 1 0 0]           %Desired amplitudes
fir= firls (n, f,a )