Hi! I'm a novice in DSP programming so I'd like to ask for some help to solve my problem. I have a signal coming from an ADC, with 10 bits of resolution and the values are converted back to voltage between 0...2.5V. The sampling rate is 500Hz. The signal has DC offset and a baseline wander which is in the sub 1 Hz range. I'd like to remove this offset and also the slow wander. My idea is was to use a simple, 2nd order HP filter, where Fc=0.3...0.5Hz. I used fdatool to calculate the coefficients and for testing purposes I created a PC app just to check the behaviour of my filter. My filter equation is the following: y(n) = x(n)*b0 + x(n-1)*b1 + x(n-2)*b2 - y(n-1)*a1 - y(n-2)*a2 where b0=1, b1=-2, b2=1, a1=-1.9962 and a2=0.9962 However I'm not getting the expected results. It won't remove the offset but introduces an even bigger offset. For example if my signal is riding on 1.00V the filtered result would ride on 3.5V, which is strange. It also doesn't smooth out the <1Hz swing. Somewhere I found an IIR implementation where y(n) has been multiplied with a gain of <1. y(n) = x(n)*b0 + x(n-1)*b1 + x(n-2)*b2 - y(n-1)*a1 - y(n-2)*a2 y(n) = gain * y(n) It eventually does work with gain=0.9999 which removed the offset and smoothed out the swings but it had a quite nasty swings before settling and it looks a bit worrying to me. So, I'd like to know if it's the proper way of high pass filtering or I am making serious mistakes here...or any tip would be welcome. Regards, Greg
Basic IIR high pass filter help
Started by ●May 28, 2012
Reply by ●May 28, 20122012-05-28
"Greg73" <gsantha@n_o_s_p_a_m.protochem.hu> writes:> Hi! > > I'm a novice in DSP programming so I'd like to ask for some help to solve > my problem. I have a signal coming from an ADC, with 10 bits of resolution > and the values are converted back to voltage between 0...2.5V. The sampling > rate is 500Hz. The signal has DC offset and a baseline wander which is in > the sub 1 Hz range. I'd like to remove this offset and also the slow > wander. > > My idea is was to use a simple, 2nd order HP filter, where Fc=0.3...0.5Hz. > I used fdatool to calculate the coefficients and for testing purposes I > created a PC app just to check the behaviour of my filter. My filter > equation is the following: > > y(n) = x(n)*b0 + x(n-1)*b1 + x(n-2)*b2 - y(n-1)*a1 - y(n-2)*a2 > > where b0=1, b1=-2, b2=1, a1=-1.9962 and a2=0.9962Greg, How are you implementing this? In Matlab? In C? Are you using floating point or fixed point? Your equation and coefficients above look correct. I plotted the response using freqz() and it looks good. Suggestions 1. What system is this running on? Is it on some SDK? Or just on a PC with Matlab? Show us your actual implementation 2. Ensure the signs on a1 and a2 of your implementation are correct. As I said, they look correct in your filter equation, but double check your implementation. 3. If you're implementing the above equation in fixed-point, show us your coefficients and code. There can be very nasty non-linearities due to quantization when implementing such a filter in fixed-point. 4. Ensure you're converting the output to 10 bits appropriately. --Randy -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●May 28, 20122012-05-28
Hi Randy, and thank you for your help. Please see my answers below. I really appreciate your help.>How are you implementing this? In Matlab? In C? Are you using floating >point or fixed point?C# .NET and I'm using floating point.> >Your equation and coefficients above look correct. I plotted the >response using freqz() and it looks good. Suggestions > > 1. What system is this running on? Is it on some SDK? Or just on a PC > with Matlab? Show us your actual implementationBasically, this is an ECG data, converted back from a 10bit ADC to floating point values between 0V...2.5V as the ADC reference of the ADC was 2.5V. The converted ECG data is in the ecg_raw array which to be filtered. private void btn_hp_Click(object sender, EventArgs e) { FILTS filts = new FILTS(); //filter struct filts.a1 = Convert.ToDouble(tb_a1.Text); //read coeffs from GUI filts.a2 = Convert.ToDouble(tb_a2.Text); filts.b0 = Convert.ToDouble(tb_b0.Text); filts.b1 = Convert.ToDouble(tb_b1.Text); filts.b2 = Convert.ToDouble(tb_b2.Text); filts.gain = Convert.ToDouble(tb_gain.Text); chart1.Series["fRaw"].Points.Clear(); //Clear chart and result array ecg_raw_hp.Clear(); ecg_raw_hp.Add(0.0); //zero the first two array memmbers to avoid false indexing Exception ecg_raw_hp.Add(0.0); for (int i = 2; i < ecg_raw.Count; i++) //start from the second member { filts.input_2=ecg_raw[i-2]; filts.input_1=ecg_raw[i-1]; filts.input=ecg_raw[i]; filts.output_2=ecg_raw_hp[i-2]; filts.output_1=ecg_raw_hp[i-1]; ecg_raw_hp.Add(HP_IIR(filts)); //call IIR filter method (see below) chart1.Series["fRaw"].Points.AddXY((double)i, ecg_raw_hp[i]); //display result point on chart } clb_plot.SetItemChecked(2, true); //GUI specific not relevant } private double HP_IIR(FILTS fs) { double y; y = fs.b0 * fs.input + fs.b1 * fs.input_1 + fs.b2 * fs.input_2 - fs.a1 * fs.output_1 - fs.a2 * fs.output_2; return y*fs.gain; } Here you can check 2 screenshots of the issue. The blue signal is the converted floating point raw data, the red is the filtered. with no gain manipulation (gain = 1) http://kepfeltoltes.hu/120528/ecg_www.kepfeltoltes.hu_.png with gain = 0.9999 http://kepfeltoltes.hu/120528/ecg2_www.kepfeltoltes.hu_.png Best regards, Greg
Reply by ●May 28, 20122012-05-28
"Greg73" <gsantha@n_o_s_p_a_m.protochem.hu> writes:> Hi Randy, > > and thank you for your help. Please see my answers below. I really > appreciate your help. > >>How are you implementing this? In Matlab? In C? Are you using floating >>point or fixed point? > > C# .NET and I'm using floating point. >> >>Your equation and coefficients above look correct. I plotted the >>response using freqz() and it looks good. Suggestions >> >> 1. What system is this running on? Is it on some SDK? Or just on a PC >> with Matlab? Show us your actual implementation > > Basically, this is an ECG data, converted back from a 10bit ADC to floating > point values between 0V...2.5V as the ADC reference of the ADC was 2.5V. > The converted ECG data is in the ecg_raw array which to be filtered. > > private void btn_hp_Click(object sender, EventArgs e) > { > > FILTS filts = new FILTS(); > //filter struct > > filts.a1 = Convert.ToDouble(tb_a1.Text); //read coeffs > from GUI > filts.a2 = Convert.ToDouble(tb_a2.Text); > > filts.b0 = Convert.ToDouble(tb_b0.Text); > filts.b1 = Convert.ToDouble(tb_b1.Text); > filts.b2 = Convert.ToDouble(tb_b2.Text); > > filts.gain = Convert.ToDouble(tb_gain.Text); > > chart1.Series["fRaw"].Points.Clear(); //Clear chart and > result array > ecg_raw_hp.Clear(); > > ecg_raw_hp.Add(0.0); //zero the > first two array memmbers to avoid false indexing Exception > ecg_raw_hp.Add(0.0); > > > > for (int i = 2; i < ecg_raw.Count; i++) //start from the > second member > { > filts.input_2=ecg_raw[i-2]; > filts.input_1=ecg_raw[i-1]; > filts.input=ecg_raw[i]; > > filts.output_2=ecg_raw_hp[i-2]; > filts.output_1=ecg_raw_hp[i-1]; > > ecg_raw_hp.Add(HP_IIR(filts)); > //call IIR filter method (see below) > chart1.Series["fRaw"].Points.AddXY((double)i, > ecg_raw_hp[i]); //display result point on chart > } > clb_plot.SetItemChecked(2, true); > //GUI specific not relevant > } > > private double HP_IIR(FILTS fs) > { > double y; > > y = fs.b0 * fs.input + > fs.b1 * fs.input_1 + > fs.b2 * fs.input_2 - > fs.a1 * fs.output_1 - > fs.a2 * fs.output_2; > > > > return y*fs.gain; > } > > > Here you can check 2 screenshots of the issue. The blue signal is the > converted floating point raw data, the red is the filtered. > > with no gain manipulation (gain = 1) > > http://kepfeltoltes.hu/120528/ecg_www.kepfeltoltes.hu_.png > > with gain = 0.9999 > > http://kepfeltoltes.hu/120528/ecg2_www.kepfeltoltes.hu_.png > > Best regards, > > GregGreg, 1. Remove the fs.gain step and focus on getting your filter correct. This is not helping you find your problem. Your filter coefficients and equation are correct as you provided them initially, so there is something else wrong somewhere and adding a gain step obfuscates the issue. 2. Put a printf() inside HP_IIR and print out these values: fs.b0 fs.b1 fs.b2 fs.a1 fs.a2 fs.gain (or just check them inside the debugger). See if they match what you expect. 3. If you still can't get things to work, try changing your HP_IIR to this: private double HP_IIR(FILTS fs) { double y; y = fs.b0 * fs.input + fs.b1 * fs.input_1 + fs.b2 * fs.input_2; return y; } This is also a highpass filter, albeit one with a wider stop band. However, it should perfectly filter out DC if b0 = 1, b1 = -2, and b2 = 1. -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●May 28, 20122012-05-28
On Mon, 28 May 2012 11:05:45 -0400, Randy Yates wrote:> "Greg73" <gsantha@n_o_s_p_a_m.protochem.hu> writes: > >> Hi Randy, >> >> and thank you for your help. Please see my answers below. I really >> appreciate your help. >> >>>How are you implementing this? In Matlab? In C? Are you using floating >>>point or fixed point? >> >> C# .NET and I'm using floating point. >>> >>>Your equation and coefficients above look correct. I plotted the >>>response using freqz() and it looks good. Suggestions >>> >>> 1. What system is this running on? Is it on some SDK? Or just on a PC >>> with Matlab? Show us your actual implementation >> >> Basically, this is an ECG data, converted back from a 10bit ADC to >> floating point values between 0V...2.5V as the ADC reference of the ADC >> was 2.5V. The converted ECG data is in the ecg_raw array which to be >> filtered. >> >> private void btn_hp_Click(object sender, EventArgs e) >> { >> >> FILTS filts = new FILTS(); >> //filter struct >> >> filts.a1 = Convert.ToDouble(tb_a1.Text); //read >> coeffs >> from GUI >> filts.a2 = Convert.ToDouble(tb_a2.Text); >> >> filts.b0 = Convert.ToDouble(tb_b0.Text); filts.b1 = >> Convert.ToDouble(tb_b1.Text); filts.b2 = >> Convert.ToDouble(tb_b2.Text); >> >> filts.gain = Convert.ToDouble(tb_gain.Text); >> >> chart1.Series["fRaw"].Points.Clear(); //Clear chart >> and >> result array >> ecg_raw_hp.Clear(); >> >> ecg_raw_hp.Add(0.0); //zero >> the >> first two array memmbers to avoid false indexing Exception >> ecg_raw_hp.Add(0.0); >> >> >> >> for (int i = 2; i < ecg_raw.Count; i++) //start from >> the >> second member >> { >> filts.input_2=ecg_raw[i-2]; >> filts.input_1=ecg_raw[i-1]; >> filts.input=ecg_raw[i]; >> >> filts.output_2=ecg_raw_hp[i-2]; >> filts.output_1=ecg_raw_hp[i-1]; >> >> ecg_raw_hp.Add(HP_IIR(filts)); >> //call IIR filter method (see below) >> chart1.Series["fRaw"].Points.AddXY((double)i, >> ecg_raw_hp[i]); //display result point on chart >> } >> clb_plot.SetItemChecked(2, true); >> //GUI specific not relevant >> } >> >> private double HP_IIR(FILTS fs) >> { >> double y; >> >> y = fs.b0 * fs.input + >> fs.b1 * fs.input_1 + >> fs.b2 * fs.input_2 - >> fs.a1 * fs.output_1 - >> fs.a2 * fs.output_2; >> >> >> >> return y*fs.gain; >> } >> >> >> Here you can check 2 screenshots of the issue. The blue signal is the >> converted floating point raw data, the red is the filtered. >> >> with no gain manipulation (gain = 1) >> >> http://kepfeltoltes.hu/120528/ecg_www.kepfeltoltes.hu_.png >> >> with gain = 0.9999 >> >> http://kepfeltoltes.hu/120528/ecg2_www.kepfeltoltes.hu_.png >> >> Best regards, >> >> Greg > > Greg, > > 1. Remove the fs.gain step and focus on getting your filter correct. > This is not helping you find your problem. Your filter coefficients > and equation are correct as you provided them initially, so there is > something else wrong somewhere and adding a gain step obfuscates the > issue. > > 2. Put a printf() inside HP_IIR and print out these values: > > fs.b0 > fs.b1 > fs.b2 > fs.a1 > fs.a2 > fs.gain > > (or just check them inside the debugger). See if they match what you > expect. > > 3. If you still can't get things to work, try changing your HP_IIR to > this: > > private double HP_IIR(FILTS fs) > { > double y; > > y = fs.b0 * fs.input + > fs.b1 * fs.input_1 + > fs.b2 * fs.input_2; > > return y; > } > > This is also a highpass filter, albeit one with a wider stop band. > However, it should perfectly filter out DC if b0 = 1, b1 = -2, and b2 > = 1.I would do this if things were inexplicably going wonky, but only to verify that the rest of the system were working correctly. Clearly if that doesn't fix it then the problem isn't the filter -- if it does, I'd go back to the question of why in heck my filter isn't working correctly, perhaps by taking the filter out of the surrounding code and testing it in isolation. -- My liberal friends think I'm a conservative kook. My conservative friends think I'm a liberal kook. Why am I not happy that they have found common ground? Tim Wescott, Communications, Control, Circuits & Software http://www.wescottdesign.com
Reply by ●May 28, 20122012-05-28
Tim Wescott <tim@seemywebsite.com> writes:> On Mon, 28 May 2012 11:05:45 -0400, Randy Yates wrote: > >> "Greg73" <gsantha@n_o_s_p_a_m.protochem.hu> writes: >> >>> Hi Randy, >>> >>> and thank you for your help. Please see my answers below. I really >>> appreciate your help. >>> >>>>How are you implementing this? In Matlab? In C? Are you using floating >>>>point or fixed point? >>> >>> C# .NET and I'm using floating point. >>>> >>>>Your equation and coefficients above look correct. I plotted the >>>>response using freqz() and it looks good. Suggestions >>>> >>>> 1. What system is this running on? Is it on some SDK? Or just on a PC >>>> with Matlab? Show us your actual implementation >>> >>> Basically, this is an ECG data, converted back from a 10bit ADC to >>> floating point values between 0V...2.5V as the ADC reference of the ADC >>> was 2.5V. The converted ECG data is in the ecg_raw array which to be >>> filtered. >>> >>> private void btn_hp_Click(object sender, EventArgs e) >>> { >>> >>> FILTS filts = new FILTS(); >>> //filter struct >>> >>> filts.a1 = Convert.ToDouble(tb_a1.Text); //read >>> coeffs >>> from GUI >>> filts.a2 = Convert.ToDouble(tb_a2.Text); >>> >>> filts.b0 = Convert.ToDouble(tb_b0.Text); filts.b1 = >>> Convert.ToDouble(tb_b1.Text); filts.b2 = >>> Convert.ToDouble(tb_b2.Text); >>> >>> filts.gain = Convert.ToDouble(tb_gain.Text); >>> >>> chart1.Series["fRaw"].Points.Clear(); //Clear chart >>> and >>> result array >>> ecg_raw_hp.Clear(); >>> >>> ecg_raw_hp.Add(0.0); //zero >>> the >>> first two array memmbers to avoid false indexing Exception >>> ecg_raw_hp.Add(0.0); >>> >>> >>> >>> for (int i = 2; i < ecg_raw.Count; i++) //start from >>> the >>> second member >>> { >>> filts.input_2=ecg_raw[i-2]; >>> filts.input_1=ecg_raw[i-1]; >>> filts.input=ecg_raw[i]; >>> >>> filts.output_2=ecg_raw_hp[i-2]; >>> filts.output_1=ecg_raw_hp[i-1]; >>> >>> ecg_raw_hp.Add(HP_IIR(filts)); >>> //call IIR filter method (see below) >>> chart1.Series["fRaw"].Points.AddXY((double)i, >>> ecg_raw_hp[i]); //display result point on chart >>> } >>> clb_plot.SetItemChecked(2, true); >>> //GUI specific not relevant >>> } >>> >>> private double HP_IIR(FILTS fs) >>> { >>> double y; >>> >>> y = fs.b0 * fs.input + >>> fs.b1 * fs.input_1 + >>> fs.b2 * fs.input_2 - >>> fs.a1 * fs.output_1 - >>> fs.a2 * fs.output_2; >>> >>> >>> >>> return y*fs.gain; >>> } >>> >>> >>> Here you can check 2 screenshots of the issue. The blue signal is the >>> converted floating point raw data, the red is the filtered. >>> >>> with no gain manipulation (gain = 1) >>> >>> http://kepfeltoltes.hu/120528/ecg_www.kepfeltoltes.hu_.png >>> >>> with gain = 0.9999 >>> >>> http://kepfeltoltes.hu/120528/ecg2_www.kepfeltoltes.hu_.png >>> >>> Best regards, >>> >>> Greg >> >> Greg, >> >> 1. Remove the fs.gain step and focus on getting your filter correct. >> This is not helping you find your problem. Your filter coefficients >> and equation are correct as you provided them initially, so there is >> something else wrong somewhere and adding a gain step obfuscates the >> issue. >> >> 2. Put a printf() inside HP_IIR and print out these values: >> >> fs.b0 >> fs.b1 >> fs.b2 >> fs.a1 >> fs.a2 >> fs.gain >> >> (or just check them inside the debugger). See if they match what you >> expect. >> >> 3. If you still can't get things to work, try changing your HP_IIR to >> this: >> >> private double HP_IIR(FILTS fs) >> { >> double y; >> >> y = fs.b0 * fs.input + >> fs.b1 * fs.input_1 + >> fs.b2 * fs.input_2; >> >> return y; >> } >> >> This is also a highpass filter, albeit one with a wider stop band. >> However, it should perfectly filter out DC if b0 = 1, b1 = -2, and b2 >> = 1. > > I would do this if things were inexplicably going wonky, but only to > verify that the rest of the system were working correctly. Clearly if > that doesn't fix it then the problem isn't the filter -- if it does, I'd > go back to the question of why in heck my filter isn't working correctly, > perhaps by taking the filter out of the surrounding code and testing it > in isolation.I coded this in Matlab and was getting similarly weird results! Scratched my head awhile, then decided to check the roots of b and a. Bingo! a has a root at DC. So it was cancelling with the zero at DC in b (b factors to simple (z-1)(z-2)). But when implemented in direct form 1, it was apparently playing havoc with the numerical precision, even with Matlab doubles! -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●May 28, 20122012-05-28
On Mon, 28 May 2012 08:25:14 -0400, Randy Yates wrote:> "Greg73" <gsantha@n_o_s_p_a_m.protochem.hu> writes: > >> Hi! >> >> I'm a novice in DSP programming so I'd like to ask for some help to >> solve my problem. I have a signal coming from an ADC, with 10 bits of >> resolution and the values are converted back to voltage between >> 0...2.5V. The sampling rate is 500Hz. The signal has DC offset and a >> baseline wander which is in the sub 1 Hz range. I'd like to remove this >> offset and also the slow wander. >> >> My idea is was to use a simple, 2nd order HP filter, where >> Fc=0.3...0.5Hz. I used fdatool to calculate the coefficients and for >> testing purposes I created a PC app just to check the behaviour of my >> filter. My filter equation is the following: >> >> y(n) = x(n)*b0 + x(n-1)*b1 + x(n-2)*b2 - y(n-1)*a1 - y(n-2)*a2 >> >> where b0=1, b1=-2, b2=1, a1=-1.9962 and a2=0.9962 > > Greg, > > How are you implementing this? In Matlab? In C? Are you using floating > point or fixed point? > > Your equation and coefficients above look correct. I plotted the > response using freqz() and it looks good. Suggestions > > 1. What system is this running on? Is it on some SDK? Or just on a PC > with Matlab? Show us your actual implementation > > 2. Ensure the signs on a1 and a2 of your implementation are correct. > As I said, they look correct in your filter equation, but double check > your implementation. > > 3. If you're implementing the above equation in fixed-point, show us > your coefficients and code. There can be very nasty non-linearities > due to quantization when implementing such a filter in fixed-point. > > 4. Ensure you're converting the output to 10 bits appropriately. > > --RandyRandy, you missed something... If a1 = -1.9962 and a2 = 0.9962, then the filter has poles at z = 0.9962 and z = 1, and as such will be metastable. One of the zeros (since they're both at z = 1) will mask the metastable pole -- but that pole will still have a transient response, which is what's leaving the massive offset. Greg: Where did you get your coefficients, and what are you really trying to achieve? (I.e, are you aiming for a Butterworth response, Chebycheve, "anything as long as it solves my problem", etc.). Take the roots of the polynomial 1 + a1 + a2. If either root has an absolute value greater than one, then your filter is unstable. If the roots form a complex pair then your algorithm, with the chosen coefficients, is correct. If the roots are both real, then you are marginally better off implementing your filter as a pair of 1-st order highpass filters in cascade. You can get away with doing it the way you are, but it is bad style an folks like Randy and I will always be pointing it out to you and saying "Y'know, you could make that filter more numerically stable by...". I think that you were trying for a filter with a pair of poles at 0.9981, which gives you a1 = -1.9962 and a2 = 0.99620361. But you saw that minuscule little 0.00000361, and thought you could just round it off. You couldn't -- polynomial roots are quite sensitive to coefficients, and get ever more so in the lower-order coefficients of high-order polynomials (assuming that the polynomial is normalized). This sensitivity to roots is reflected in the behavior of IIR filters, and is why I would recommend that you split the thing up into a pair of 1st-order sections. At least keep all your digits, even if you require more precision to resolve your a2 exactly than you have digits in IEEE single-precision floating point (log_2(1e-8) is -26.6 or so, you have 25 bits to express that -- so sorry, can't get there from here...) -- My liberal friends think I'm a conservative kook. My conservative friends think I'm a liberal kook. Why am I not happy that they have found common ground? Tim Wescott, Communications, Control, Circuits & Software http://www.wescottdesign.com
Reply by ●May 28, 20122012-05-28
Randy Yates <yates@digitalsignallabs.com> writes:> [...] > (b factors to simple (z-1)(z-2)).Bah! "(b factors to simply (z-1)(z-1))". -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●May 28, 20122012-05-28
Tim Wescott <tim@seemywebsite.com> writes:> On Mon, 28 May 2012 08:25:14 -0400, Randy Yates wrote: > >> "Greg73" <gsantha@n_o_s_p_a_m.protochem.hu> writes: >> >>> Hi! >>> >>> I'm a novice in DSP programming so I'd like to ask for some help to >>> solve my problem. I have a signal coming from an ADC, with 10 bits of >>> resolution and the values are converted back to voltage between >>> 0...2.5V. The sampling rate is 500Hz. The signal has DC offset and a >>> baseline wander which is in the sub 1 Hz range. I'd like to remove this >>> offset and also the slow wander. >>> >>> My idea is was to use a simple, 2nd order HP filter, where >>> Fc=0.3...0.5Hz. I used fdatool to calculate the coefficients and for >>> testing purposes I created a PC app just to check the behaviour of my >>> filter. My filter equation is the following: >>> >>> y(n) = x(n)*b0 + x(n-1)*b1 + x(n-2)*b2 - y(n-1)*a1 - y(n-2)*a2 >>> >>> where b0=1, b1=-2, b2=1, a1=-1.9962 and a2=0.9962 >> >> Greg, >> >> How are you implementing this? In Matlab? In C? Are you using floating >> point or fixed point? >> >> Your equation and coefficients above look correct. I plotted the >> response using freqz() and it looks good. Suggestions >> >> 1. What system is this running on? Is it on some SDK? Or just on a PC >> with Matlab? Show us your actual implementation >> >> 2. Ensure the signs on a1 and a2 of your implementation are correct. >> As I said, they look correct in your filter equation, but double check >> your implementation. >> >> 3. If you're implementing the above equation in fixed-point, show us >> your coefficients and code. There can be very nasty non-linearities >> due to quantization when implementing such a filter in fixed-point. >> >> 4. Ensure you're converting the output to 10 bits appropriately. >> >> --Randy > > Randy, you missed something...So did you on your first pass... ;) -- Randy Yates Digital Signal Labs http://www.digitalsignallabs.com
Reply by ●May 28, 20122012-05-28
> >Greg: > >Where did you get your coefficients, and what are you really trying to >achieve? (I.e, are you aiming for a Butterworth response, Chebycheve, >"anything as long as it solves my problem", etc.). >I generated those coeffs using the fdatool in MATLAB. There are no strict specifications for the type of the filter or any response - I just need a HP filter to remove DC and baseline wander thus I can count the R peaks more accurately.>Take the roots of the polynomial > >1 + a1 + a2. > >If either root has an absolute value greater than one, then your filter >is unstable. > >If the roots form a complex pair then your algorithm, with the chosen >coefficients, is correct. > >If the roots are both real, then you are marginally better off >implementing your filter as a pair of 1-st order highpass filters in >cascade. You can get away with doing it the way you are, but it is bad >style an folks like Randy and I will always be pointing it out to you and>saying "Y'know, you could make that filter more numerically stableby...".>I'll try that. You mean I just feed F1's output into F2's input?>I think that you were trying for a filter with a pair of poles at 0.9981,>which gives you a1 = -1.9962 and a2 = 0.99620361. But you saw that >minuscule little 0.00000361, and thought you could just round it off. >You couldn't -- polynomial roots are quite sensitive to coefficients, and>get ever more so in the lower-order coefficients of high-order >polynomials (assuming that the polynomial is normalized). > >This sensitivity to roots is reflected in the behavior of IIR filters, >and is why I would recommend that you split the thing up into a pair of >1st-order sections. At least keep all your digits, even if you require >more precision to resolve your a2 exactly than you have digits in IEEE >single-precision floating point (log_2(1e-8) is -26.6 or so, you have 25 >bits to express that -- so sorry, can't get there from here...)I think that the rounding was what caused my problem. Initially, MATLAB gave me these coeffs: a1:-1.996240592289382 a2:0.99625345892156669 I was not aware that the digits beyond 1E-04 have such importance. So I just cut them off quite dumbly. Quickly I changed the coeffs and I'm getting the expected behaviour now: http://i47.tinypic.com/2mgsbkp.gif Thank you very Tim and Randy!
