2
$\begingroup$

I've implemented a very simple first order recursive low pass filter in c using the algorithm:

double a, b, prevOutput;

double lowPassFilter(double input)
{
    double output = a * input + b * prevOutput;
    prevOutput = output; 
    return output;
}

where:

a = 1.0 - b

According to The Scientist and Engineer's Guide to Digital Signal Processing, the cutoff frequency of the filter is defined by the relation:

$$b = e^{-2πf_c}$$

$f_c$ is the -3db cutoff frequency as a fraction of the sampling frequency. Using $f_c = 0.25$, I calculated that $b = 0.20788$ (approximately) and plugged it into my filter. I wanted to verify my filter, so I passed an impulse function through it and recorded the output:

double result[128];
for (int i = 0; i < 128; ++i)
{
    if (0 == i) result[i] = lowPassFilter(1.0);
    else result[i] = lowPassFilter(0.0);
}

I then performed FFT on the result to find the frequency response:

complex_t fftResult* = fft(result);
for (int i = 0; i < 65; ++i)
{
    double amplitude = fftResult[i].real * fftResult[i].real;
    amplitude += fftResult[i].img * fftResult[i].img;
    amplitude = sqrt(magnitude);
    printf("f = %.3f: %f\n", i / 128.0, amplitude );
}

The result I saw was somewhat expected but also not exactly what I expected. In the DC bin ($f = 0$) I saw a amplitude very close to one and I saw the amplitude attenuate as it approached 0.5. However, I also expected that at $f = 0.25$ the amplitude should have been around $0.70795$ (i.e. -3db). but instead the amplitude was around $0.775$:

frequency response of the lp filter

 ...
f = 0.242: 0.782501
f = 0.250: 0.774971
f = 0.258: 0.767380
 ...

Aside from the inefficient code, is there something wrong with the method I am using to obtain the frequency response or am I interpreting the cut-off frequency relation incorrectly or something else?

$\endgroup$
6
$\begingroup$

That formula for the cut-off frequency is a very inaccurate approximation. In this answer I derived the exact relation between the coefficient of a first order recursive averaging filter and its 3-dB cut-off frequency. Note that in the quoted answer I used the constant $\alpha=1-b$. From formula $(3)$ in that answer we get for the coefficient $b$

$$b = 2-\cos(\omega_c)-\sqrt{(2-\cos(\omega_c))^2-1)}\tag{1}$$

With $\omega_c=2\pi f_c$ we get for $f_c=0.25$ a value $b=0.26795$.

The frequency response of the filter is given by

$$H(e^{j\omega})=\frac{1-b}{1-be^{-j\omega}}\tag{2}$$

The figure below shows the magnitude response $|H(e^{j\omega})|$, and it shows that the desired cut-off frequency is achieved with the chosen value of $b$:

enter image description here

Note that the 'EDIT' part of the answer quoted above also explains where the inaccurate formula for the cut-off frequency comes from. The following figure shows the resulting cut-off frequency of the filter as a function of the desired cut-off frequency when the exact formula $(1)$ is used and when the approximation is used:

enter image description here

Clearly, the approximation is quite bad and it only kind of works for very small cut-off frequencies. I see no reason why it should be used at all.

$\endgroup$
  • $\begingroup$ +1 I was looking for some mistakes in the code and missed the approximation part... $\endgroup$ – msm Oct 22 '16 at 8:28
  • $\begingroup$ @msm: Yes, I happen to have stumbled upon that strange formula before (see the other answer quoted in my answer). $\endgroup$ – Matt L. Oct 22 '16 at 8:34
  • $\begingroup$ I have sometimes used the approximation for parameter smoothing as it is easy to construct (up to a constant multiple of frequency) just by thinking how the slope of the impulse invariant exponentially decaying continuous impulse response is related to the cutoff frequency, and for that application it is accurate enough. $\endgroup$ – Olli Niemitalo Oct 22 '16 at 9:51
  • $\begingroup$ Tried various other fc values using the provided equation and it indeed works perfectly. The last chart in this answer is a really nice touch, it seems that for fc < 0.1 the approximation works much better. $\endgroup$ – DaveS Oct 24 '16 at 1:46
  • 1
    $\begingroup$ @gramm: It depends on how you define $f_c$. I use $\omega=2\pi f$, which means that $f=0.5$ corresponds to Nyquist, so $f_c\ge 0.5$ doesn't make sense. (And indeed, $f_c=1$ is the same as $f_c=0$ due to periodicity). $\endgroup$ – Matt L. Aug 20 at 12:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.