lowpass filter - Fast Integer 8 Hz 2nd Order LP for Microcontroller

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?