DSPRelated.com
Forums

How to write a NotchFilter procedure

Started by gpdsoft June 5, 2006
I'm not an expert of signal processing so I wonder if there is a formula or
an algoritm to cutoff the 50Hz from my signal...

something like:
CutFreq:=50;
SampleFreq:=50000;
SampleNum:=5000;

NotchFilter(MySignal,cutFreq,SampleFreq,SampleNum);

I've used the "DewResearch" library but the notch filter on 50Hz cut a lot
of frequency (from 10Hz to 150Hz)... It works really bad!

The signal I have to filter is the 50Hz (or 60Hz) of the electrical
system.

So Please help me!

gpdsoft wrote:
> I'm not an expert of signal processing so I wonder if there is a formula or > an algoritm to cutoff the 50Hz from my signal... > > something like: > CutFreq:=50; > SampleFreq:=50000; > SampleNum:=5000; > > NotchFilter(MySignal,cutFreq,SampleFreq,SampleNum);
Do you know what a biquad filter is? If so, you can filter the data with a biquad which has transfer function H(z) = (1 - 2 cos(w) z^-1 + z^-2)/(1 - 2 a cos(w) z^-1 + a z^-1), where w = CutFreq / SampleFreq, and 0 < a < 1 is a parameter that determines the width of the notch. As a -> 1, the notch becomes narrower. Regards, Andor
Andor wrote:

> gpdsoft wrote: > >>I'm not an expert of signal processing so I wonder if there is a formula or >>an algoritm to cutoff the 50Hz from my signal... >> >>something like: >>CutFreq:=50; >>SampleFreq:=50000; >>SampleNum:=5000; >> >>NotchFilter(MySignal,cutFreq,SampleFreq,SampleNum); > > > Do you know what a biquad filter is? If so, you can filter the data > with a biquad which has transfer function > > H(z) = (1 - 2 cos(w) z^-1 + z^-2)/(1 - 2 a cos(w) z^-1 + a z^-1), > > where > > w = CutFreq / SampleFreq, > > and 0 < a < 1 is a parameter that determines the width of the notch. As > a -> 1, the notch becomes narrower. > > Regards, > Andor >
It's not that easy if you want the high frequency gain to match the low frequency gain -- the 'w' in the denominator needs to be skewed from the notch frequency. I haven't been able to find a tidy closed-form solution (I have found some very messy ones, however). Perhaps Peter's spreadsheet* has a nice solution, or at least a good explanation of the messy one. * It's a measure of my compulsive alliteration that I spent about a minute trying to think of a synonym for 'spreadsheet' or 'formula' that started with a 'P'. Sigh. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html
Anonymous wrote:

>>>One tip: Rather than try and design a notch at 50 Hz it should be easier > > to > >>>put the notch at dc and then mix the input signal down 50 Hz before the >>>notch and back up after the notch. >> >>Hmmmm.... >> >>I am not sure this would be a good idea. Using a notch filter, not a >>high >>pass filter, implies that there is useful data in the band DC -- 50 Hz. >>Mixing down by 50 Hz would cause aliasing, which you would not >>be able to recover from. >> >> > > > I'm assuming that the processing/mixer would be complex so there should be > no alias problem.
How do you get the quadrature signal needed for the complex mix if DC needs to be preserved? It takes a long time to delay DC by 90 degrees. Jerry -- Engineering is the art of making what you want from things you can get. &#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;
gpdsoft wrote:

> I'm not an expert of signal processing so I wonder if there is a formula or > an algoritm to cutoff the 50Hz from my signal... > > something like: > CutFreq:=50; > SampleFreq:=50000; > SampleNum:=5000; > > NotchFilter(MySignal,cutFreq,SampleFreq,SampleNum); > > I've used the "DewResearch" library but the notch filter on 50Hz cut a lot > of frequency (from 10Hz to 150Hz)... It works really bad! > > The signal I have to filter is the 50Hz (or 60Hz) of the electrical > system. > > So Please help me! >
Actually, we've given you at least two. 1. "theirs" -- Make an IIR notch filter. The algorithm to do this is: a. Web search on "IIR filter", find a site that looks sensible b. Learn how to write IIR filter code c. Web search on "Notch filter" for coefficients, or use Peter's d. Write your notch filter e. Use it. 2. "mine" -- Do a post-facto removal of 50 Hz from your data a. Estimate the 50Hz content of your data i. Multiply your data by cos(50Hz), sum the result, call it I ii. Multiply your data by sin(50Hz), sum the result, call it Q b. Calculate the 50Hz tone that's buried in your data: i. Tone = 2/N * (I cos(50Hz) + Q sin(50Hz)) c. Subtract it out: data - Tone = data_without_50Hz -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html
Tim Wescott wrote:

  ...

> * It's a measure of my compulsive alliteration that I spent about a > minute trying to think of a synonym for 'spreadsheet' or 'formula' that > started with a 'P'. Sigh.
Peter's presentation? Jerry -- Engineering is the art of making what you want from things you can get. &#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;
Tim Wescott wrote:
> Andor wrote: > > > gpdsoft wrote: > > > >>I'm not an expert of signal processing so I wonder if there is a formula or > >>an algoritm to cutoff the 50Hz from my signal... > >> > >>something like: > >>CutFreq:=50; > >>SampleFreq:=50000; > >>SampleNum:=5000; > >> > >>NotchFilter(MySignal,cutFreq,SampleFreq,SampleNum); > > > > > > Do you know what a biquad filter is? If so, you can filter the data > > with a biquad which has transfer function > > > > H(z) = (1 - 2 cos(w) z^-1 + z^-2)/(1 - 2 a cos(w) z^-1 + a z^-1), > > > > where > > > > w = CutFreq / SampleFreq, > > > > and 0 < a < 1 is a parameter that determines the width of the notch. As > > a -> 1, the notch becomes narrower. > > > > Regards, > > Andor > > > It's not that easy if you want the high frequency gain to match the low > frequency gain -- the 'w' in the denominator needs to be skewed from the > notch frequency. I haven't been able to find a tidy closed-form > solution (I have found some very messy ones, however).
There is something in r b-j's cookbook: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt If finite attenuation is enough (versus the true zero that the notch filter has at w), a peaking filter will also work as a notch filter.
> Perhaps Peter's > spreadsheet* has a nice solution, or at least a good explanation of the > messy one. > > * It's a measure of my compulsive alliteration that I spent about a > minute trying to think of a synonym for 'spreadsheet' or 'formula' that > started with a 'P'. Sigh.
Perhaps Peter's paper presents a practical plan? :-)
Tim Wescott skrev:
> Andor wrote: > > > gpdsoft wrote: > > > >>I'm not an expert of signal processing so I wonder if there is a formula or > >>an algoritm to cutoff the 50Hz from my signal... > >> > >>something like: > >>CutFreq:=50; > >>SampleFreq:=50000; > >>SampleNum:=5000; > >> > >>NotchFilter(MySignal,cutFreq,SampleFreq,SampleNum); > > > > > > Do you know what a biquad filter is? If so, you can filter the data > > with a biquad which has transfer function > > > > H(z) = (1 - 2 cos(w) z^-1 + z^-2)/(1 - 2 a cos(w) z^-1 + a z^-1), > > > > where > > > > w = CutFreq / SampleFreq, > > > > and 0 < a < 1 is a parameter that determines the width of the notch. As > > a -> 1, the notch becomes narrower. > > > > Regards, > > Andor > > > It's not that easy if you want the high frequency gain to match the low > frequency gain -- the 'w' in the denominator needs to be skewed from the > notch frequency. I haven't been able to find a tidy closed-form > solution (I have found some very messy ones, however). Perhaps Peter's > spreadsheet* has a nice solution, or at least a good explanation of the > messy one. > > * It's a measure of my compulsive alliteration that I spent about a > minute trying to think of a synonym for 'spreadsheet' or 'formula' that > started with a 'P'. Sigh.
Phormula? Rune
On Mon, 05 Jun 2006 12:51:37 -0500, "gpdsoft" <gpdsoftware@yahoo.it>
wrote:

>Hi, >I have do write a procedure in Delphi that take an array of samples and >filters it. > >The filter I've to do is Notch to cut out the 50Hz from my signal. > >My signal is sampled at 50Khz for a max of 250 millisec. > >Do you know where I can find the code or something similar to quickly >implement this function in my software? > >Thnx a lot
Hello gpdsoft, Implement a filter whose time-domain expression is: y(n) = x(n) + Kx(n-1) + x(n-2). Variable y(n) is the filter's output, and variable x(n) is the filter's input. The value for K is: K = -2cos(2*pi*fn/fs) where fn is 50 Hz, and fs is 50 kHz. Good Luck, [-Rick-]
If this is to cut out DC, isn't it the second harmonic that sounds
dissonant, and therefore you should be looking at a notch for 50Hz?


Rick Lyons wrote:
> On Mon, 05 Jun 2006 12:51:37 -0500, "gpdsoft" <gpdsoftware@yahoo.it> > wrote: > > >Hi, > >I have do write a procedure in Delphi that take an array of samples and > >filters it. > > > >The filter I've to do is Notch to cut out the 50Hz from my signal. > > > >My signal is sampled at 50Khz for a max of 250 millisec. > > > >Do you know where I can find the code or something similar to quickly > >implement this function in my software? > > > >Thnx a lot > > Hello gpdsoft, > > Implement a filter whose time-domain > expression is: > > y(n) = x(n) + Kx(n-1) + x(n-2). > > Variable y(n) is the filter's output, and variable > x(n) is the filter's input. > > The value for K is: > > K = -2cos(2*pi*fn/fs) > > where fn is 50 Hz, and fs is 50 kHz. > > Good Luck, > [-Rick-]