Fast Integer 8 Hz 2nd Order LP for Microcontroller

Question

I need an 8 Hz 2nd order LP filter. It does not have to be terribly accurate but it should be 2nd order-ish. Performance is much more important.

My sample rate at the moment is ~9.5 kHz. The filter should have a Q of 0.707 but again accuracy is not paramount.

I have been looking at biquad filters but this is for a little 8 bit microcontroller. It has floating point emulation but my understanding is that it's 10x slower. So I think 32 bit integers are about as far as I want to push it. Trying to do an 8 Hz biquad using 32 bit integers could prove difficult.

If I understand correctly these filters as just summing sample values multipled by some coefficient at different taps in a delay line so it seems to me there should be a way to implement a relatively simple filter by picking 3 taps and 3 integer coefficients and then use a divisor at the end of each step to give me a stable running value.

Does this make sense? If yes, can someone recommend a design procedure? I think I have GNU octave on my machine.

UPDATE 1:

While Dan's description of a moving average filter is very interesting, I believe it would be fobidden in my particular application.

In short, I'm creating an LED VU meter. Yes, it is trivial to make an LED meter that will give the operator an indication as to the level of a signal and it might even be easy to make it vaguely accurate. Currently my implementation uses a simple running average and it works.

However, being the obsessive compulsive type, I have taken an interest in how mechanical VU meters actually work and someone directed me to a paper on the topic "A model of the VU (volume-unit) meter, with speech applications" Lobdell 2006. According to the paper, the procedure is fairly simple:

full wave rectify -> 8 Hz 2nd order low pass -> scale -> log convert

That someone also specifically rejected the idea of using a moving average filter because it would not yield accurate results with all frequencies being metered.

So I think I will need to implement a proper biquad.

But it has been impressed apon me by the great answers posted here as to the importance of reducing the sampling rate (decimation). Unforunately again, the decimation method cannot be a running average. It must be a time domain calculation to properly represent the energy of all frequencies being metered.

UPDATE 2:

When I do this (Octave / Matlab):

clear all;
close all;

N = 4096;
Z = 16;
A = 256

m = [];

for mi = 1:A
    ip=round(rand(N,1)*512);
    op=[];
    rval = 0;
    for i = 1:N
        rval = (rval * (Z-1) + ip(i)) / Z;
        if mod(i, Z) == 0
            op = [op rval];
        end
    endfor

    op = fft(op,N/Z);
    m = [m; abs(op(2:N/Z/2))];
endfor

x = [2:N/Z/2];
y = mean(m);

semilogx(x,20*log10(y))

I get this:

There are no notches in the response as long as the decimated sample rate is the same as the moving average window size (meaning window size is also modulus in test to save sample). If it's not exactly the same, I do get horrible notches in the response.

So this is good. Is this what you (Dan) meant by how things would "fold" together? Near as I can tell the slope is about 3 dB / octave so if I do 4 passes that's a not terribly accurate fast integer 2nd order LP.

Can I depend on this behavior when translated to fixed-point math?

UPDATE 3:

I now have a 12db / octave filter. I decimated by a factor of 16 to reduce the sample rate to ~600 sps and used a basic 2nd order biquad:

clear all;
close all;

N = 4096;
A = 8;
Fs = 9500;
T = 1/Fs;
t = 0:T:(N-1)*T;
f = Fs*(2:(N/2))/N;

% decimated values
Z = 16;
dN = N/Z;
dFs = Fs/Z;
dT = 1/dFs;
df = dFs*(2:(dN/2))/dN;

m = [];

for mi = 1:A
    ip=round(randn(1,N)*512 + sin(2*pi*100*t)*256);
    op=[];

    n = N;

    rval = 0;
    v0 = 0;
    v1 = 0;
    v2 = 0;
    for i = 1:n
        rval = (rval * (Z-1) + ip(i)) / Z;
        if mod(i, Z) == 0
            v0 = v1;
            v1 = v2;
            v2 = (87 * rval - 27191 * v0 + 59613 * v1) / 32768;
            tmp = (v0 + v2) + 2 * v1;
            op = [op tmp];
        end
    endfor

    n = length(op);

    op = fft(op,n);
    m = [m; abs(op(2:n/2))];
endfor

y = mean(m);

semilogx(df,20*log10(y))

However, with 4096 samples I get this:

and with 65536 samples I get this:

as you can see the magnitude is drifting. Does this mean the filter is not stable?

How can I stablize it?

Answer 1

This answer provides a quick introduction to decimation concepts and CIC filters which I would consider as one solution given the description.

Bottom Line First

---

Given your use of a microcontroller, (implied emphasis on minimizing resources), and that you indicated you do not need a high performance filter- consider doing everything with Cascade-Integrator-Comb (CIC) structures: They have unity gain coefficients and require NO multipliers! They can also provide something very important you should consider doing: decimation.

Decimate by 64:

Skipping to a possible solution for you, first decimate your signal using a 2nd Order CIC decimator such as the structure I show below, which will take your 9.5KHz signal down to a $9.5KHz/8= 1187.5$ Hz sample rate. Follow this with an identical CIC decimator to reduce the sampling rate to 148.4 Hz.

The decimator above does low pass filter the input signal, with the equivalent frequency response plotted below, relative to the input sampling rate (1 on the Normalized Frequency domain is half the sampling rate ($F_s/2$ or 4.75 KHz.). This is a pre-filter but no-where near the 8Hz cutoff desired; The first null appears at 148.4 Hz and the -3dB rejection point is at 47.5 Hz. The upper envelope of the $\textrm{sinc}^2$ response shown approximates a $1/f^2$ roll-off or 2nd order response.

Simple Moving Average LPF:

With this lower rate signal, design your 8 Hz LPF. You could consider here a simplest solution of using just a cascade of 4 sample moving averages, with some bit shifting to scale back the level, as follows:

This structure will provide the response shown in the figure below, which is in addition to the broadband filtering rejection from the decimator plotted above. On the Normalized Frequency axis, 1 represents half the sampling rate, or 74.2 Hz, and 8 Hz is at 0.107. The attenuation at 8 Hz is 2.2 dB. The advantage is the simplicity (No multipliers! Unity gain coefficients!) If this filter is insufficient, then you can easily design higher performance filters, but that task is much easier at this lower sampling rate. If the concern is the passband droop, that can be compensated easily with a pre-emphasis inverse $\textrm{sinc}$ FIR filter using just 3 taps (I added a link at the bottom with details on that if interested).

Here is another interesting option, as it has two fewer taps and better overall stopband rejection! This is the same structure above, except replace two of the four-tap FIRs with three-tap FIRs. This positions a null in the region of the dominant stopband area.

Your performance requirements may be different than what is provided by the result above; the above is to provide an example of what can be done, quite simply, with a minimum of resources, and notably no multipliers. Even a high performance implementation would use the decimation portion of this approach, (or half-band decimating filters which are also quite efficient), followed by a shaping filter which would provide the final filtering requirements desired in addition to compensation for the droop from the $\textrm(sinc)^n$ shape in the passband.

---

Background

CIC

Brief notes on the CIC structure shown above. (The CIC filter is also referred to as the Hogenauer Filter based on Eugene Hogenauer's 1981 paper introducing the idea). This CIC structure is two "integrators" in cascade which are simple accumulators, followed by a downsampler, which throws away all but every 8th sample, followed by two "comb filters", which is a subtraction operation between two successive lower frequency samples as shown in the figure. This one is specifically a 2nd order CIC, $(CIC)^2$, and the order can be increased by adding more integrator and comb stages (5 in a row of each would be a $(CIC)^5$ which is what you often see used in high performance systems. The higher the order, the more rejection of frequency regions that would fold into the frequency of interest (see discussion on decimation below). The accumulators cycle delay block shown ($z^{-1}$) is one clock cycle at the higher input rate, and the comb cycle delay block shown is one clock cycle at the lower output rate. This decimator will work with any value of downsampling, but as that number increases the precision needed within the accumulator and comb elements increases. I used quick judgement without detailing an actual design that it would be easier to cascade two of these overall structures to get a total decimation of 64, rather than implement it in one down-sampler element, but the number of stages and the decimation value at each stage is a knob you can play with. Important to note that the accumulators are allowed to overflow by design, but they must "roll-over" on overflow for the CIC to work properly (as opposed to saturate), meaning $cnt_{max}+n = n$. Additionally they cannot be allowed to roll-over multiple times before the subsequent subtraction operations have completed. Basically the math for the subsequent subtraction needs to be consistent regardless of a previous accumulator overflow: If the difference between samples is allowed to grow larger than the precision, the an inconsistency will result due to modulo arithmetic. For this reason I suggest you start with the 2nd order (or even a first order) as I have sketched unless you have a thorough understanding of the CIC operation.

The CIC is just a mathematical transformation of the moving average filter! It is exactly equivalent in operation to a moving average filter, it is not an approximation. In fact, the structure below is the CIC equivalent of the 4 sample moving average filter shown as "FLT" in the second figure above. Cascading multiple stages of integrators and combs is the same result as cascading multiple moving average filters (we cascaded two stages for the $CIC^2$ implementation.) What also is clear by reviewing the structure below (after convincing yourself the input and outputs would be identical to that with a 4 sample scaled moving average filter, or divide the output by 4 to get a true moving average), is if we were to use this structure as the pre-filter for a decimate by 4 (see discussion on decimation below), our first inclination would be to filter the signal with the filter as provided, and then down sample by keeping every 4th sample. If you study the "comb" portion of the filter, you can see that the downsample can be moved prior to the comb and the comb can be reduced to one delay at the lower rate- with the same result! (taking advantage of the fact that we were ignoring all but every 4th sample anyway).

Why Decimate?

Regardless of final filter approach used, the fact that you want to make an 8 Hz low pass filter in a system with a 9.5 KHz sampling rate where resources are prime, should absolutely be decimated (meaning reduce the sampling rate digitally). This will significantly simplify your filter implementation for a given performance metric. To give some insight into this, I included a slide I had on estimating the number of taps (complexity) needed for a FIR filter. The key take-away is to notice how the number of taps is driven by the transition band, as a fraction of the sampling rate. Given the same transition in Hz, if you implement your filter with a sampling rate of 80 Hz versus 9.5 KHz, the number of taps needed to achieve the same filter decreases by a factor of nearly 120!! You may need to use the higher sampling rate elsewhere in your implementation, which is fine, but for the filter, and subsequent processing after your filter, you should highly consider running at a lower rate (multi-rate signal processing).

To properly decimate, you must filter the input signal first and then throw away samples (downsample). The filtering is important as there are bands in your overall frequency range that will fold right into your low pass signal of interest when you throw away samples, so very important to reject noise and interference at these locations. A very simple and effective way to do this, it to filter with a moving average, as it is just a summation with no multipliers and provides it's maximum rejection at all the frequency locations that would fold in. For example, if you want to decimate by 2 (reduce the sampling rate by 2), then you can sum adjacent samples in a moving average, and then throw away every other answer in the result; which is identical to doing a block by block average with 2 sample blocks. To decimate by 10, you can do a moving average over 10 samples, and then throw away all but every 10th sample in the result; which is identical to doing a block by block average with 10 sample blocks. Key point: A moving average filter (unity gain coeff FIR) has a $\textrm{sinc}$ function frequency response in frequency, with the first null at $1/T$ Hz, and subsequent nulls in integer multiples of $1/T$, where $T$ is the length of your filter in time. These null locations ARE the locations where the noise and interference folds into your "DC" area when you downsample. A CIC filter properly rejects these nulls and provides the necessary moving average for decimation. The passband droop which results can be easily compensated for with an inverse $\textrm{sinc}$ pre-emphasis filter. This inverse-$\textrm{sinc}$ filter can be easily done with a 3 tap FIR.

CIC Filter Compensation

I found this question and added an answer there detailing the simple 3-tap FIR filter I use for CIC compensation

how to make CIC compensation filter