# Linear phase Notch filter / Band reject filter implementation in C++

I'd like to implement a digital notch filter, which has linear phase, in C++ or C#. Example : a notch filter that removes 400hz frequency on a .wav file (16 bit, 44.1 Khz, stereo). I'm a bit lost with many examples, I don't achieve to implement them. Maybe some of you has a working example?

Thank you in advance.

• See equation 19-8 on this page: dspguide.com/ch19/3.htm – Paul R Oct 12 '13 at 14:20
• @PaulR thanks, but I already have some good theory books... The problem is that I don't know achieve to implement them with a practical code in C++... Would you have a working example in C++? Thanks in advance. – Basj Oct 12 '13 at 16:17
• The page I linked to gives you the formula for the coefficients for a simple notch filter - all you have to do is calculate the coefficients from the formula and then apply that. If you don't know how to implement a filter once you have the coefficients then I'll post some code later. – Paul R Oct 12 '13 at 16:47
• @PaulR I understand the maths, the coefficient computation, but I don't know how to make a C++ code running on WAV files with such a filter... – Basj Oct 12 '13 at 18:28
• OK - I've put some example C/C++ code in an answer below which shows how to implement a simple recursive filter using coefficients a0, a1, a2, b1, b2 - this can be used for any 2nd order IIR filter, but if you use coefficients from the link above you'll get a notch filter. – Paul R Oct 12 '13 at 20:58

You can implement a simple recursive filter to process blocks of samples like this:

void filter(const int *x, int *y, int n)
{
static float x_2 = 0.0f;                    // delayed x, y samples
static float x_1 = 0.0f;
static float y_2 = 0.0f;
static float y_1 = 0.0f;

for (i = 0; i < n; ++i)
{
y[i] = a0 * x[i] + a1 * x_1 + a2 * x_2  // IIR difference equation
+ b1 * y_1 + b2 * y_2;
x_2 = x_1;                              // shift delayed x, y samples
x_1 = x[i];
y_2 = y_1;
y_1 = y[i];
}
}


Notch filter coefficients a0, a1, a2, b1, b2 can be calculated using eq 19-8 on this page: http://dspguide.com/ch19/3.htm

• Thanks a lot PaulR, I will try to paste this in my code! Is it possible to include the computation of a0, a1, a2, b1, b2 directly in the code, then the only parameter would be f (frequency) ? (with some math C++ functions) ? Last question : do you think this method can be used in order to do super-precise notch filter like this one (super high Q value) (example using a VST) dl.dropboxusercontent.com/u/83031018/forum/… – Basj Oct 12 '13 at 21:16
• The coefficients would typically be calculated once, when the centre frequency is defined or changed, and then they are just passed into the filtering function as parameters. I don't know what you mean by "super-precise" - a filter is just something that matches its specification - what aspect of the filter specification are you particularly concerned about ? – Paul R Oct 13 '13 at 7:04
• I need linear phase because of further application of my algorithm (too long to explain in this small comment ;) )... How to increase the order of the filter ? – Basj Oct 13 '13 at 9:47
• The filter above is IIR and so not linear phase - you'll need to design an FIR filter with a large number of terms if you want a high Q notch and linear phase. – Paul R Oct 13 '13 at 21:24
• Assuming this is for a fixed frequency you will need a filter design program to generate the coefficients. It's going to be tedious though - if you really do need high Q and linear phase then the filter may well have hundreds of coefficients and the design process could be slow and painful. I don't know how big a filter MATLAB can design, but you might try it, otherwise you may need a specific filter design program. At this point I would ask on DSP.SE as your problem has become quite challenging now. – Paul R Oct 14 '13 at 5:44

If you usually use Matlab for analysis the files, you might want to consider using Simulink and then just generate the code of the filter. Especially if this is something you run in real time (ie inputting audio and processing it on the fly).

Another approach is using the coefficients as was suggested above.Once you have those you can use the following formula:

output[i] = A[0] * in[2] + A[1] * in[1] + A[2] * in[0] - B[1] * out[1] - B[2] * out[0];


Where A/B are you coefficients, in[0,1,2] are 3 points of you original file and out [0,1,2] is the filtered one. Of-course for the first few points there are no values of out[0,1] and you should erase them after filtration.

By running this formula and shifting both in and out at each point you will result in your filtered signal.

See this thread for further explanation

• You are right! Fixed it... – PolinaK Oct 12 '13 at 21:18

I found a proper C-Code example on this page

http://www.dsprelated.com/showcode/173.php

It helped me a lot, so maybe some of you are also interested...

• What is "proper C-Code"? Or are you somehow suggesting that doing it in C++ is improper!? – Clifford Apr 27 '16 at 16:52