Programming a low-pass filter

Question

I have programmed a Sega Master System emulator in Java (although this question is not Java specific) and have finished everything except for the SN76489 sound chip. How this chip makes sound is easy enough - the problem I have is converting it to a form that is playable through a PC/laptop/whatever the JVM is running on.

I have identified the steps as follows;

As the SN76489 runs at a sample rate of roughly 221khz - this means the waves it outputs have a frequency of up to 110khz (although in practice I doubt anything ever goes this high). I need to therefore implement a low-pass filter before I downsample it.

I then want to downsample it to 44.1khz, so I can output it through an audio line (in my case a Java Source Data Line).

For this I need to set the low-pass filter at 22.05khz, but the problem is I have no idea (mathematically speaking) how a low-pass filter actually works. I need this in order to be able to write one.

At present, my sound chip creates a 0.2 second buffer, and stores samples in this, at 221khz as mentioned above. I can downsample, as I understand this - but if I downsample without a low-pass filter being applied first, I understand that I may potentially get aliasing glitches in the resulting sound stream.

Could anyone recommend the simplest mathematically minded algorithm for doing this - I realise due to the variables involved that a low-pass is never 'exact', but I just need a sound explanation that is dummed down enough for my brain (which has not really dealt with wave processing before now) to understand.

If it helps, the specifics are: The SN76489 outputs three square waves and one noise channnel simultaneously. These are summed together, and output to the mixer/amplifier - this stage in the chain is where I want to run the low-pass filter before I downsample and then amplify the wave. Any help people can give me is much appreciated. I realise background reading is required but I want to know 'what' I need to read. Many thanks.

UPDATE: I came up with a simpler approach in the end - still not quite there though. The SN76489 works by generating each tone channel from a register value - a polarity of 1 is output, the value decremented, and so on - until the value is 0, then the value is reset and the polarity switched to - 1, and so on. This value is then multiplied by a volume to get the final amplitude for that sample, and summed with the other channels.

I now simply prevent a register value that will produce a square wave above the nyquist limit I require from being produced. This leaves me with a much better signal, but it still has some buzzing/popping in it - not sure why as the maximum possible frequency should be 18,473Hz. Could this popping/buzzing be because when the chip switches a channel from one frequency to the next, it doesn't allow the current wave form to completely finish? As an example, the chip outputs 1111, then 00 - instead of the full four zeros and switches to a new freqency - this could cause aliasing could it not?

Answer 1

EDIT: I have included a filter implementation below that answers your question as asked. However, implementing a high order filter with a signal at such a high sample rate is going to consume many, many millions of operations per second. It might be best if you do a spectrum analysis of the output of the chip first. If there is no sound energy above a few kHz then the anti-aliasing filter is a waste of processing resources. Even if there is energy up to moderately high frequencies, it may be worth decimating the signal first, then filtering before applying the second stage decimation. As a side note, you might also want to decimate down to a much lower rate than 44.1 kHz. You'd probably be fine with an 8- or 10-kHz sample rate for a master system emulator (we're not talking hi-fi here). But anyway, to answer your question about how to implement a low-pass filter with the sampling rate and cutoff that you specify . . .

Ok, firstly to design a lowpass filter. The matlab decimate function sounds good to my ears so we will copy that approach for this example. The documentation says the following

The decimated vector y is r times shorter in length than the input vector x. By default, decimate employs an eighth-order lowpass Chebyshev Type I filter with a cutoff frequency of 0.8*(Fs/2)/r. It filters the input sequence in both the forward and reverse directions to remove all phase distortion, effectively doubling the filter order.

Cheby filters are a good choice as they have a steeper rejection than Butterworth designs at the expense of a little passband ripple. We cannot do IIR filtering in both directions in a real time system, but this should be ok for your purposes. We can make the filter coefficients using the following Matlab code . . . .

sr = 221e3;
srDesired = 44.1e3;
order = 8;
passBandRipple = 1; %//dB

Wp = 0.8 * (srDesired/2) / (sr/2);

[b,a] = cheby1 (order, passBandRipple, Wp, 'low');

freqz(b,a,[],sr, 'half');

sos = tf2sos(b,a)

This gives us an 8th order IIR filter with a response like the following. This looks like what we want. The phase response is unimportant for this application. The cutoff is 0.8* 22.050 kHz as you want the signal near to the Nyquist limit to be well attenuated before decimation.

The tf2sos command at the end converts the filter that we have just created into second order sections that you can realise using a cascade of biquad filter sections. The output of this command is as follows . . .

SECTION A

b=1.98795003258633e-07, 3.97711540624783e-07, 1.98854354149782e-07,
a=1 -1.81843900641769, 0.835282840946310

SECTION B

b=1, 2.02501937393162, 1.02534004997240,
a=1, -1.77945624664044, 0.860871442492022

SECTION C

b=1, 1.99938921206706, 0.999702296645625,
a=1, -1.73415546937221, 0.907015729252152

SECTION D

b=1, 1.97498006006623, 0.975285456788754,
a=1, -1.72600279649390, 0.966884508765457

Now, you can use these filter coefficients for each biquad stage in the filter cascade. You can realise this filter using code like the that in the following example. It is C code, but you should be able to convert it to java pretty easily. Note that there is no a0 coefficient in the code below. The second order sections above are correctly normalised so that a0 always equals 1. Just leave it out.

//Make the following vars private members of your filter class
// b0, b1, b2, a1, a2 calculated above
// m1, m2 are the memory locations
// dn is the de-denormal coeff (=1.0e-20f) 

void processBiquad(const float* in, float* out, unsigned length)
{
    for(unsigned i = 0; i < length; ++i)
    {
        register float w = in[i] - a1*m1 - a2*m2 + dn;
        out[i] = b1*m1 + b2*m2 + b0*w;
        m2 = m1; m1 = w;
    }
    dn = -dn;
}

You should make a class for this filter and then instantiate 4 separate classes (1 for each filter), setting the a, and b values to those specified above. Next, hook the input of one stage to the ouput of the next to give you your cascade.

Answer 2

http://en.wikipedia.org/wiki/Butterworth_filter

C code generator: http://www-users.cs.york.ac.uk/~fisher/mkfilter/ (should be translatable to Java easily enough)