Back to the Averager (and a bit of digital filter groundwork for lunch)

Well, last night I wanted to noodle around with some more sounds and get the basic KS on a little further so ended up deep in digital filter land, more of which later. I ended up the evening having implemented all the digital filter classes in the STK as they proved to be a good set of building blocks and an easy way to bring the concepts together in an hour or so. I'm going to blog on that next, but before going there I wanted to tie-up the lose end I mentioned with the Averager and in doing so take the opportunity to dip into digital filter land. So, this is kind of an explanation and pre-step to the next journey.

Intuitively it's obvious that a basic two point moving averager is a low-pass filter, and we did the hand-waving bit by putting the phase delay to be 0.5 without taking that any further. Before going into digital filters proper, it's worth understanding these two points (ha ha) as the thinking is good to explore the more complicated stuff later on.

As an online introduction reference for digital filters, keeping with the STK and audio theme, the best reference is Julius O Smith's Introduction to Digital Filters with Audio Applications.

So, stepping straight in, we're going to go for the general case here to setup some tools and then when I blog about the STK filter implementation we'll use those tools and build upwards. Constraining the problem, we're only going to look at Linear Time Invariant (LTI) filters. The general form we'll be looking at are Infinite Impulse Response (IIR) filters.

IIR filter
In block diagram terms they look like the diagram below. In simple terms there are a set of delayed inputs, each scaled by a value b added together with a set of delayed outputs (feedback) each scaled by a value a.

What this means is many of the different forms of delay and filter (a delay is just a kind of filter) we have already seen and will go on to look at can all be generalised to the same model. A simplification of this model where there is no feedback present and all a values are set to zero is called a Finite Impulse Response (FIR) filter. For example the simple non-interplotating delay line could be considered as a FIR filter where the delay items b0,b1,.... b(N-1) are all zero and bN is given the value 1. This would just delay the input by the number of samples N.

Back to our two point moving averager, this is actually a FIR filter with b0 = 0.5 and b1 = 0.5. Explaining this, b0 is the current input and we take the last input (delayed by one sample) add the two together and divide by two (averaging) or as we can see, mix half of each together.

So, in our Averager class implementation as this is a subclass of Filter we can set the two parameters b[0] and b[1] = 0.5 which will allow us to use some of the other tools we're going to look at. It will also mean that the overriden phaseDelay function can be removed and the general one in Filter can be used instead as we'll see a bit later on.

        public Averager()
        {
            b = new List<double>() { 0.5, 0.5 };
        }

Before going further, we should have a look at how the block diagram seems in notation which will help us navigate some of the maths going on and the code algorithms after that. The difference equation and transfer function for an IIR are shown below:

As there are a lot of symbols here to start with, I've taken the small simplification of making a0=1 which is the usual case as it looks slightly neater and will not get in the way where we're going. Please consult one of the references to dig deeper if needed.

What the equations are saying is that the output is a sum of the scaled (by b-values) of delayed inputs x(n-p) along with scaled (by a-values) of the feedback outputs y[n-q].

So if you set all the a's to zero half of the equation disappears like the block diagram. 

Frequency Response

Now, what we're interested in with filters is getting some kind of frequency response. What this means is that for certain frequencies we would like the input signal to have the amplitude of those frequencies scaled by a certain amount. I'm assuming here that this grounding of that is in place, so will not distract us by going into details in this post. Now, this is the usual departure point to go into lots of maths, but I'll leave that for the textbooks. What we want is some way numerically (e.g. in code) of being able to plug a set of these a and b parameters into a model and look at how the filter affects certain frequencies to see what effect we're getting.

Robin Schmit has a nice little document on his site with some of the equations we need. The one below calculates the magnitude response of a filter at a certain frequency:

This is basically taking apart bits of the equations and considering the phasors at a particular frequency and seeing how they are scaled. I'm not going to go any deeper with this as the equation is enough to get us back into code.

The other companion equation is the phase response.

This is like the magnitude response, but shows how the phase of a certain input frequency signal is affected by the filter.

Now, bringining this together. Remember that equation phaseDelay that we glossed over in the Filter class? Hmmm, well, looks pretty similar to the phase response, right? Nope, well take a look, it's summing up the reala and imaga (which is ca and sa in the above equations) and then doing a little atan2 calculation.

Well, we're nearly there, the eagle eyed will have noticed that this equation is the phase response and the code is calculating the phaseDelay. Well, that's simply just a shift in terms to see how much delay the phase shift means, the concepts are interchangable. 

So, we'd better get back into some code and put in this magnitude response function, which can be added to Filter:

        public virtual double magnitude(double frequency)
        {
            if (frequency <= 0.0 || frequency > 0.5 * PcmStreamSource.SampleRate)
            {
                throw new Exception();
            }

            double omegaT = 2 * Math.PI * frequency / PcmStreamSource.SampleRate;
            double reala = 1.0;
            double imaga = 0.0;
            double realb = 0.0;
            double imagb = 0.0;
 
            for (int i = 0; i < a.Count; i++)
            {
                reala += a[i] * Math.Cos(i * omegaT);
                imaga -= a[i] * Math.Sin(i * omegaT);
            }

            for (int i = 0; i < b.Count; i++)
            {
                realb += b[i] * Math.Cos(i * omegaT);
                imagb -= b[i] * Math.Sin(i * omegaT);
            }

            reala *= gain;
            realb *= gain;
            imaga *= gain;
            imagb *= gain;

            double mag = Math.Sqrt((Math.Pow(realb,2) + Math.Pow(imagb,2)) / 
                                   (Math.Pow(reala,2) + Math.Pow(imaga,2)));

            return mag;
        }

With this in place we're ready to rock and roll and take a look at the magnitude and phase response for the averager and any of these type of filters that can be made from this in the future.

It's a bit kludgy, but here's a simple algorithm to show using the functions with the averager and putting the values into a long text string. Why? Well it's an easy way to get the values out and into some graphing tools.

         private void AnalyseResponse()
        {
            Averager filter = new Averager();

            string phase = "";
            string amp = "";

            int SAMPLES = 100;
            for (int i = 1; i < SAMPLES; i++)
            {
                double ph = filter.phaseDelay(PcmStreamSource.SampleRate / 2.0 * ((double)i / (double)SAMPLES));
                double a = filter.amplitude(PcmStreamSource.SampleRate / 2.0 * ((double)i / (double)SAMPLES));

                phase += String.Format("{0:0.00},", ph);
                amp += String.Format("{0:0.00},", a);
            }
        }

Having promised to explain the averager response, I'm not quite going to get there in this post. So, finishing up for lunch, I'll quickly post how to show this response in another blog note.