matlab - NAudio real-time equalizer

I have designed 9 IIR bandpass filters (chebyshev type 1, order 4) using fdatool in Matlab. Then I use a and b filter coefficients to apply it on differential equation. So, my question is how to apply (real-time) gain of each bandpass filter using sliders? ( how to add gain or loss to certain filter?)

Here is my code in C#:

public void initializeFilters(){
    filters.Clear();
    //9 Band pass filters
    for (int i = 0; i < 9; i++)
        filters.Add(new Filter());

    //F1=40hz Fc=50Hz F2=60Hz
    float k = (float)Math.Pow(1.0, -5);
    filters[0].B = new float[] { k * 0.1992f, 0, k * -0.3983f, 0, k * 0.1992f }; 

    filters[0].A = new float[] { 1.0000f, -3.9968f, 5.9904f, -3.9905f, 0.9969f };

    //60Hz 80Hz 100Hz
    k = (float)Math.Pow(1.0, -4);
    filters[1].B = new float[] { k * 0.0795f, 0, k * -0.1591f, 0, k * 0.0795f };
    filters[1].A = new float[] { 1.0000f, -3.9935f, 5.9807f, -3.9810f, 0.9938f };

    //100Hz 130Hz 160Hz
    filters[2].B = new float[] { k * 0.1787f, 0, k * -0.3574f, 0, k * 0.1787f };
    filters[2].A = new float[] { 1.0000f, -3.9899f, 5.9705f, -3.9713f, 0.9907f };


    //160Hz 350Hz 600Hz
    filters[3].B = new float[] { 0.0009f, 0, -0.0019f, 0, 0.0009f };
    filters[3].A = new float[] { 1.0000f, -3.9255f, 5.7847f, -3.7927f, 0.9335f };

    //600Hz 1300Hz 2000Hz
    filters[4].B = new float[] { 0.0088f, 0, -0.0176f, 0, 0.0088f };
    filters[4].A = new float[] { 1.0000f, -3.7187f, 5.2467f, -3.3315f, 0.8040f };

    //2000Hz 4000Hz 6000Hz

    filters[5].B = new float[] { 0.0595f, 0, -0.1191f, 0, 0.0595f };
    filters[5].A = new float[] { 1.0000f, -2.8765f, 3.4528f, -2.0797f, 0.5459f };

    //6000Hz 8000Hz 10000Hz
    filters[6].B = new float[] { 0.0595f, 0, -0.1191f, 0, 0.0595f };
    filters[6].A = new float[] { 1.0000f, -1.4273f, 1.8140f, -1.0319f, 0.5459f };

    //10000Hz 12000Hz 14000Hz
    filters[7].B = new float[] { 0.0595f, 0, -0.1191f, 0, 0.0595f };
    filters[7].A = new float[] { 1.0000f, 0.4731f, 1.3375f, 0.3420f, 0.5459f };


    //14000Hz 16000Hz 18000Hz
    filters[8].B = new float[] { 0.0595f, 0, -0.1191f, 0, 0.0595f };
    filters[8].A = new float[] { 1.0000f, 2.2239f, 2.5782f, 1.6078f, 0.5459f };
    }
    public override int Read(byte[] buffer, int offset, int count)
    {
    int read = sourceStream.Read(buffer, offset, count);
    float[] f_read = new float[read / 4], y = new float[read / 4];
    Buffer.BlockCopy(buffer, offset, f_read, 0, count);


    for (int n = 0; n < read / 4; n++)
        y[n] = 0;

    for (int i = 0; i < filters.Count; i++)
    {
        for (int n = 0; n < read / 4; n++)
        {
            for (int k = 0; k < filters[i].B.Length; k++)
            {

                if (n - k >= 0)
                    y[n] = y[n] + filters[i].B[k] * f_read[n - k];
            }
            for (int k = 1; k < filters[i].A.Length; k++)
            {
                if (n - k >= 0)
                    y[n] = y[n] - filters[i].A[k] * y[n - k];
            }
        }
    }


    for (int n = 0; n < read / 4; n++)
        y[n]= Math.Min(1, Math.Max(-1, y[n]));

    Buffer.BlockCopy(y, 0, buffer, offset, read / 4);
    return read;
    }

Here is GUI of the sliders: (0% - certain filter does not effect the signal, 100% - fully effect the signal) gui sliders

EDIT 1 I have redesigned my bandpass filters, using IIR chebyshev type I, order 2:

    Fc Flow Fhigh [in Hz]
50 35 71  b = [0.0050,0,-0.0050] a= [1.0000,-1.9899,0.9900]
80 57 113 b=[0.0078,0,-0.0078] a=[1.0000,-1.9843,0.9844]
130 92 184 b=[0.0127,0,-0.0127] a=[1.0000,-1.9742,0.9746]
350 247 495 b=[0.0336,0,-0.0336] a=[1.0000,-1.9305,0.9329]
1300 919 1838 b=[0.1141,0,-0.1141] a=[1.0000,-1.7414,0.7717]
4000 2828 5657 b=[0.2865,0,-0.2865] a=[1.0000,-1.1984,0.4270]
8000 5657 11314 b=[0.4559,0,-0.4559] a=[1.0000,-0.4188,0.0882]
12000 8485 16971 b=[ 0.5758,0,-0.5758] a=[1.0000,0.2477,-0.1517]
16000 11314 22040 b=[0.6532,0,-0.6532] a=[1.0000,0.6927,-0.3063]

and here is my new code:

        private void initializeFilters()
    {
        filters = new List();
        filters.Add(new Filter(new float[] { 0.0050f, 0f, -0.0050f }, new float[] { 1.0000f, -1.9899f, 0.9900f }));
        filters.Add(new Filter(new float[] { 0.0078f, 0f, -0.0078f }, new float[] { 1.0000f, -1.9843f, 0.9844f }));
        filters.Add(new Filter(new float[] { 0.0127f, 0f, -0.0127f }, new float[] { 1.0000f, -1.9742f, 0.9746f }));
        filters.Add(new Filter(new float[] { 0.0336f, 0f, -0.0336f }, new float[] { 1.0000f, -1.9305f, 0.9329f }));
        filters.Add(new Filter(new float[] { 0.1141f, 0f, -0.1141f }, new float[] { 1.0000f, -1.7414f, 0.7717f }));

        filters.Add(new Filter(new float[] { 0.2865f, 0f, -0.2865f }, new float[] { 1.0000f, -1.1984f, 0.4270f }));
        filters.Add(new Filter(new float[] { 0.4559f, 0f, -0.4559f }, new float[] { 1.0000f, -0.4188f, 0.0882f }));
        filters.Add(new Filter(new float[] { 0.5758f, 0f, -0.5758f }, new float[] { 1.0000f, 0.2477f, -0.1517f }));
        filters.Add(new Filter(new float[] { 0.6532f, 0f, -0.6532f }, new float[] { 1.0000f, 0.6927f, -0.3063f }));
    }
    public void process(float[] x, ref float[] y, float value, int i)
    {
        float f = ((100 * value) / 20)/100; //mapping values to 0.0 - 1.0
        float[] tmp = new float[y.Length];
        for (int n = 0; n < x.Length; n++)

        {
            for (int k = 0; k < filters[i].B.Length; k++)
            {
                if (n - k >= 0)
                {
                    tmp[n] += filters[i].B[k] * x[n - k];
                    y[n] += f * tmp[n];
                }
            }
            for (int k = 1; k < filters[i].A.Length; k++)

            {
                if (n - k >= 0)
                {
                    tmp[n] -= filters[i].A[k] * tmp[n - k];
                    y[n] -= f * tmp[n];
                }
            }
        }
    }
    public override int Read(byte[] buffer, int offset, int count)

    {
        int read = sourceStream.Read(buffer, offset, count);
        float[] f_read = new float[read / 4], y = new float[read / 4];
        Buffer.BlockCopy(buffer, offset, f_read, 0, count);

        int st = 1;
        for (int i = 0; i < Form1.trackBarValue.Length; i++)
        {
            if (Form1.trackBarValue[i] > 0)
            {

                process(f_read, ref y, Form1.trackBarValue[i], i);
                st++;
            }
        }
        if (st == 1)
           process(f_read, ref y, 20, 0); //20 is just a maximum slider value

        for (int n = 0; n < read / 4; n++)
            y[n] = Math.Min(1, Math.Max(-1, y[n] / st));


        Buffer.BlockCopy(y, 0, buffer, offset, read);
        return read;
    }

I'm hearing a little noise in signal(a little bangs), can anybody tell me what am I doing wrong?

Answer

Bunch of things:

1. Your band definition is really odd. Some of the lower bands are really small while some of the mid bands are huge. Ideally you want octaves where the band edges are 1/sqrt(2) and sqrt(2) around the center frequency


2. You have designed bandpass filters in parallel. These are not "reconstructive" so you don't end up with a flat frequency response if you add the outputs of the band passes together.



3. You need a "parametric" or "peaking EQ" filters in cascade. These filters are defined by center frequency, gain (or cut) and Q. Q controls the bandwidth and Q=2 should give you something similar to ocatve bandwidth. See http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt .


4. You will have to update the filter coefficients every time the slider is moved. This may results in pops and clicks and there are various ways to deal with this. Post another questions if this is a problem

RE: EDIT 1

1. You are still using parallel bandpasses, so you will never get a perfectly falt frequency response if the all the band gains are "neutral"


2. The slider feels off. Typical graphical equalizers would do -15dB .. +15dB in logarithmic steps. Your slider is linear and it only provides cut


3. You have a "hard clip" to prevent overflow. That will sound pretty bad if overflow actually happens. There are better ways of dealing with that.


4. I think the main problem is your state management (or lack thereof). There is a conditional: $if (n - k >= 0)$. That's wrong, you always need to do the summation over the whole loop. In order to do this, you need to keep the last two input and output samples of the previous data block around. So x[n-1] for block K becomes x[-1] for block K+1. You need to keep track of x[-1], x[-2], y[-1] and y[-2].