Hi all, I want to implement an butterworth low pass filter to filter the DC component of a signal inside a PLL - phase locked loop. I am processing sample-by-sample. Therefore I have to filter the signal sample by sample (Not the whole signal or block by block). I have implemented a filter to filter the whole signal. But Im struggling to modify it make it work sample by sample basis. I highly appreciate if someone can advice me regarding this issue. I have included the code Im using to filter. [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter coefficient generation (Cutoff frequency is 2*fc) x=filter(b,a,y); % y is the input signal. % Here these two lines do the perfect filtering operation for the entire signal, But I want to process sample by sample (In-side a phase locked loop) Can someone please advice how to do that using available matlab function or different approach? Thank you,
Matlba butterworth IIR filter implementation - sample by sample filtering basis
Started by ●June 16, 2012
Reply by ●June 16, 20122012-06-16
On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote:> Hi all, > > I want to implement an butterworth low pass filter to filter the DC > component of a signal inside a PLL - phase locked loop. I am processing > sample-by-sample. Therefore I have to filter the signal sample by sample > (Not the whole signal or block by block). I have implemented a filter to > filter the whole signal. But Im struggling to modify it make it work > sample by sample basis. > > I highly appreciate if someone can advice me regarding this issue. I > have included the code Im using to filter. > > > [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter coefficient > generation (Cutoff frequency is 2*fc) > > x=filter(b,a,y); % y is the input signal. > > % Here these two lines do the perfect filtering operation for the entire > signal, But I want to process sample by sample (In-side a phase locked > loop) > > Can someone please advice how to do that using available matlab function > or different approach?1: I'm not sure what order Matlab uses for it's a and b vectors. I'm fairly sure that a(1) corresponds to the z^0 coefficient, a(2) corresponds to z^1, etc. -- but you need to check this. 2: You absolutely, positively, do not want to implement a 10th order IIR filter as a single direct-form filter. That's Very Bad -- you'll get no end of problems with coefficient rounding. In fact, I'm surprised you got sensible results in batch mode, and I strongly suggest that you check the Bode plot of the filter that you build with those A and B vectors. 2a: Split the filter up into five 2nd-order sections and cascade them. 3: I'm not even sure if your butter and filter functions are working in the Laplace domain or the z domain -- check. If they're in the Laplace domain (i.e., continuous time filtering), then check back here for more guidance. 4: If you are implementing this 10th order filter inside your control loop, just stop. You're doing something wrong. If you make the filter cutoff frequency low enough to be useful you will get so much phase shift that you'll kill all stability. If you make the filter cutoff high enough to maintain stability, your filter will be useless. The vectors a and b express a difference equation (if butter is designing a sampled-time filter): x(n) = y(n) * b(11) + y(n-1) + b(10) + ... + y(n-10) * b(1) + -x(n-1) * a(10) - x(n-2) * a(9) - ... - x(n-10) * a(1) This is a direct-form (type mumble mumble -- I can never keep them straight, but it runs from 1 to 4, or 1 and 2 and 1 modified and 2 modified -- it depends on the author). But you don't want to realize the above filter: you want to factor your a and b polynomials, assemble each one back into five 2nd-order polynomials, choose the b polynomials whose roots are closest to the a polynomials, and implement the filter whose transfer function is B/A = B1/A1 * B2/A2 * ... * B5/A5. In the time-domain, that means that you cascade filters 1 through 5 (i.e., run filter 1's output into filter 2, thence into filter 3, etc.). For each filter, just implement the above difference equation, only with delays that never exceed 2. -- Tim Wescott Control system and signal processing consulting www.wescottdesign.com
Reply by ●June 16, 20122012-06-16
On Sunday, June 17, 2012 12:35:49 PM UTC+10, Tim Wescott wrote:> On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote: > > > Hi all, > > > > I want to implement an butterworth low pass filter to filter the DC > > component of a signal inside a PLL - phase locked loop. I am processing > > sample-by-sample. Therefore I have to filter the signal sample by sample > > (Not the whole signal or block by block). I have implemented a filter to > > filter the whole signal. But Im struggling to modify it make it work > > sample by sample basis. > > > > I highly appreciate if someone can advice me regarding this issue. I > > have included the code Im using to filter. > > > > > > [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter coefficient > > generation (Cutoff frequency is 2*fc) > > > > x=filter(b,a,y); % y is the input signal. > > > > % Here these two lines do the perfect filtering operation for the entire > > signal, But I want to process sample by sample (In-side a phase locked > > loop) > > > > Can someone please advice how to do that using available matlab function > > or different approach? > > 1: I'm not sure what order Matlab uses for it's a and b vectors. I'm > fairly sure that a(1) corresponds to the z^0 coefficient, a(2) > corresponds to z^1, etc. -- but you need to check this. > > 2: You absolutely, positively, do not want to implement a 10th order IIR > filter as a single direct-form filter. That's Very Bad -- you'll get no > end of problems with coefficient rounding. In fact, I'm surprised you > got sensible results in batch mode, and I strongly suggest that you check > the Bode plot of the filter that you build with those A and B vectors. > > 2a: Split the filter up into five 2nd-order sections and cascade them. > > 3: I'm not even sure if your butter and filter functions are working in > the Laplace domain or the z domain -- check. If they're in the Laplace > domain (i.e., continuous time filtering), then check back here for more > guidance. > > 4: If you are implementing this 10th order filter inside your control > loop, just stop. You're doing something wrong. If you make the filter > cutoff frequency low enough to be useful you will get so much phase shift > that you'll kill all stability. If you make the filter cutoff high > enough to maintain stability, your filter will be useless. > > The vectors a and b express a difference equation (if butter is designing > a sampled-time filter): > > x(n) = y(n) * b(11) + y(n-1) + b(10) + ... + y(n-10) * b(1) + > -x(n-1) * a(10) - x(n-2) * a(9) - ... - x(n-10) * a(1) > > This is a direct-form (type mumble mumble -- I can never keep them > straight, but it runs from 1 to 4, or 1 and 2 and 1 modified and 2 > modified -- it depends on the author). > > But you don't want to realize the above filter: you want to factor your a > and b polynomials, assemble each one back into five 2nd-order > polynomials, choose the b polynomials whose roots are closest to the a > polynomials, and implement the filter whose transfer function is > > B/A = B1/A1 * B2/A2 * ... * B5/A5. > > In the time-domain, that means that you cascade filters 1 through 5 > (i.e., run filter 1's output into filter 2, thence into filter 3, etc.). > > For each filter, just implement the above difference equation, only with > delays that never exceed 2. > > -- > Tim Wescott > Control system and signal processing consulting > www.wescottdesign.comHi Tim, Thanks very much for your reply. I didn't look at the phase shift of the filter. But when I filter a signal and viewed the frequency spectrum it has filtered the higher frequency component. I choose the order as 10 because when I use low order like 2 it didnt remove the high frequency component completely. I didn't look at the body plots because the frequency response of the filter (freqz(b,a)) seems to be a non oscillating graph. I am not familiar with digital signal processing implementations. Therefore inorder to implement a required filter (to remove high frequency component 2*fc and take DC component from a sinusoidal multiplied signal) what should be the process. I highly appreciate your advice. Thank you.
Reply by ●June 17, 20122012-06-17
On Sat, 16 Jun 2012 19:52:29 -0700, Dinuka wrote:> On Sunday, June 17, 2012 12:35:49 PM UTC+10, Tim Wescott wrote: >> On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote: >> >> > Hi all, >> > >> > I want to implement an butterworth low pass filter to filter the DC >> > component of a signal inside a PLL - phase locked loop. I am >> > processing sample-by-sample. Therefore I have to filter the signal >> > sample by sample (Not the whole signal or block by block). I have >> > implemented a filter to filter the whole signal. But Im struggling to >> > modify it make it work sample by sample basis. >> > >> > I highly appreciate if someone can advice me regarding this issue. I >> > have included the code Im using to filter. >> > >> > >> > [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter >> > coefficient generation (Cutoff frequency is 2*fc) >> > >> > x=filter(b,a,y); % y is the input signal. >> > >> > % Here these two lines do the perfect filtering operation for the >> > entire signal, But I want to process sample by sample (In-side a >> > phase locked loop) >> > >> > Can someone please advice how to do that using available matlab >> > function or different approach? >> >> 1: I'm not sure what order Matlab uses for it's a and b vectors. I'm >> fairly sure that a(1) corresponds to the z^0 coefficient, a(2) >> corresponds to z^1, etc. -- but you need to check this. >> >> 2: You absolutely, positively, do not want to implement a 10th order >> IIR filter as a single direct-form filter. That's Very Bad -- you'll >> get no end of problems with coefficient rounding. In fact, I'm >> surprised you got sensible results in batch mode, and I strongly >> suggest that you check the Bode plot of the filter that you build with >> those A and B vectors. >> >> 2a: Split the filter up into five 2nd-order sections and cascade them. >> >> 3: I'm not even sure if your butter and filter functions are working in >> the Laplace domain or the z domain -- check. If they're in the Laplace >> domain (i.e., continuous time filtering), then check back here for more >> guidance. >> >> 4: If you are implementing this 10th order filter inside your control >> loop, just stop. You're doing something wrong. If you make the filter >> cutoff frequency low enough to be useful you will get so much phase >> shift that you'll kill all stability. If you make the filter cutoff >> high enough to maintain stability, your filter will be useless. >> >> The vectors a and b express a difference equation (if butter is >> designing a sampled-time filter): >> >> x(n) = y(n) * b(11) + y(n-1) + b(10) + ... + y(n-10) * b(1) + >> -x(n-1) * a(10) - x(n-2) * a(9) - ... - x(n-10) * a(1) >> >> This is a direct-form (type mumble mumble -- I can never keep them >> straight, but it runs from 1 to 4, or 1 and 2 and 1 modified and 2 >> modified -- it depends on the author). >> >> But you don't want to realize the above filter: you want to factor your >> a and b polynomials, assemble each one back into five 2nd-order >> polynomials, choose the b polynomials whose roots are closest to the a >> polynomials, and implement the filter whose transfer function is >> >> B/A = B1/A1 * B2/A2 * ... * B5/A5. >> >> In the time-domain, that means that you cascade filters 1 through 5 >> (i.e., run filter 1's output into filter 2, thence into filter 3, >> etc.). >> >> For each filter, just implement the above difference equation, only >> with delays that never exceed 2. >> >> -- >> Tim Wescott Control system and signal processing consulting >> www.wescottdesign.com > > Hi Tim, > > Thanks very much for your reply. > > I didn't look at the phase shift of the filter. But when I filter a > signal and viewed the frequency spectrum it has filtered the higher > frequency component. I choose the order as 10 because when I use low > order like 2 it didnt remove the high frequency component completely. > > I didn't look at the body plots because the frequency response of the > filter (freqz(b,a)) seems to be a non oscillating graph. > > I am not familiar with digital signal processing implementations. > Therefore inorder to implement a required filter (to remove high > frequency component 2*fc and take DC component from a sinusoidal > multiplied signal) what should be the process. > > I highly appreciate your advice. > > Thank you.Well, are you making the filter part of your loop filter, or are you filtering the phase error and using it elsewhere? That's what's important to the loop stability. A notch filter at 2*fc would do what you are asking for -- but is what you are asking for what you need? What are you trying to _do_? -- Tim Wescott Control system and signal processing consulting www.wescottdesign.com
Reply by ●June 17, 20122012-06-17
On Sunday, June 17, 2012 1:57:54 PM UTC+10, Tim Wescott wrote:> On Sat, 16 Jun 2012 19:52:29 -0700, Dinuka wrote: > > > On Sunday, June 17, 2012 12:35:49 PM UTC+10, Tim Wescott wrote: > >> On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote: > >> > >> > Hi all, > >> > > >> > I want to implement an butterworth low pass filter to filter the DC > >> > component of a signal inside a PLL - phase locked loop. I am > >> > processing sample-by-sample. Therefore I have to filter the signal > >> > sample by sample (Not the whole signal or block by block). I have > >> > implemented a filter to filter the whole signal. But Im struggling to > >> > modify it make it work sample by sample basis. > >> > > >> > I highly appreciate if someone can advice me regarding this issue. I > >> > have included the code Im using to filter. > >> > > >> > > >> > [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter > >> > coefficient generation (Cutoff frequency is 2*fc) > >> > > >> > x=filter(b,a,y); % y is the input signal. > >> > > >> > % Here these two lines do the perfect filtering operation for the > >> > entire signal, But I want to process sample by sample (In-side a > >> > phase locked loop) > >> > > >> > Can someone please advice how to do that using available matlab > >> > function or different approach? > >> > >> 1: I'm not sure what order Matlab uses for it's a and b vectors. I'm > >> fairly sure that a(1) corresponds to the z^0 coefficient, a(2) > >> corresponds to z^1, etc. -- but you need to check this. > >> > >> 2: You absolutely, positively, do not want to implement a 10th order > >> IIR filter as a single direct-form filter. That's Very Bad -- you'll > >> get no end of problems with coefficient rounding. In fact, I'm > >> surprised you got sensible results in batch mode, and I strongly > >> suggest that you check the Bode plot of the filter that you build with > >> those A and B vectors. > >> > >> 2a: Split the filter up into five 2nd-order sections and cascade them. > >> > >> 3: I'm not even sure if your butter and filter functions are working in > >> the Laplace domain or the z domain -- check. If they're in the Laplace > >> domain (i.e., continuous time filtering), then check back here for more > >> guidance. > >> > >> 4: If you are implementing this 10th order filter inside your control > >> loop, just stop. You're doing something wrong. If you make the filter > >> cutoff frequency low enough to be useful you will get so much phase > >> shift that you'll kill all stability. If you make the filter cutoff > >> high enough to maintain stability, your filter will be useless. > >> > >> The vectors a and b express a difference equation (if butter is > >> designing a sampled-time filter): > >> > >> x(n) = y(n) * b(11) + y(n-1) + b(10) + ... + y(n-10) * b(1) + > >> -x(n-1) * a(10) - x(n-2) * a(9) - ... - x(n-10) * a(1) > >> > >> This is a direct-form (type mumble mumble -- I can never keep them > >> straight, but it runs from 1 to 4, or 1 and 2 and 1 modified and 2 > >> modified -- it depends on the author). > >> > >> But you don't want to realize the above filter: you want to factor your > >> a and b polynomials, assemble each one back into five 2nd-order > >> polynomials, choose the b polynomials whose roots are closest to the a > >> polynomials, and implement the filter whose transfer function is > >> > >> B/A = B1/A1 * B2/A2 * ... * B5/A5. > >> > >> In the time-domain, that means that you cascade filters 1 through 5 > >> (i.e., run filter 1's output into filter 2, thence into filter 3, > >> etc.). > >> > >> For each filter, just implement the above difference equation, only > >> with delays that never exceed 2. > >> > >> -- > >> Tim Wescott Control system and signal processing consulting > >> www.wescottdesign.com > > > > Hi Tim, > > > > Thanks very much for your reply. > > > > I didn't look at the phase shift of the filter. But when I filter a > > signal and viewed the frequency spectrum it has filtered the higher > > frequency component. I choose the order as 10 because when I use low > > order like 2 it didnt remove the high frequency component completely. > > > > I didn't look at the body plots because the frequency response of the > > filter (freqz(b,a)) seems to be a non oscillating graph. > > > > I am not familiar with digital signal processing implementations. > > Therefore inorder to implement a required filter (to remove high > > frequency component 2*fc and take DC component from a sinusoidal > > multiplied signal) what should be the process. > > > > I highly appreciate your advice. > > > > Thank you. > > Well, are you making the filter part of your loop filter, or are you > filtering the phase error and using it elsewhere? That's what's > important to the loop stability. > > A notch filter at 2*fc would do what you are asking for -- but is what > you are asking for what you need? > > What are you trying to _do_? > > -- > Tim Wescott > Control system and signal processing consulting > www.wescottdesign.comHi, I am trying to filter the phase error in my loop filter. Actually I have used a FIR filter with 200 coefficients before to remove the unnecessary frequency components. And the PLL worked when there is no noise. But when I add noise (at SNR of 10 dB) the PLL fails to lock. So I assume that the delay in the Loop filter(Low pass filter) might cause the problem. So I decided to use a IIR filter with low length so that the delay would be low. Im sorry for my ignorance, but I am struggling to make a suitable low pass filter. Tim, can I have your email address? Thanks very much.
Reply by ●June 17, 20122012-06-17
On Sunday, June 17, 2012 1:57:54 PM UTC+10, Tim Wescott wrote:> On Sat, 16 Jun 2012 19:52:29 -0700, Dinuka wrote: > > > On Sunday, June 17, 2012 12:35:49 PM UTC+10, Tim Wescott wrote: > >> On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote: > >> > >> > Hi all, > >> > > >> > I want to implement an butterworth low pass filter to filter the DC > >> > component of a signal inside a PLL - phase locked loop. I am > >> > processing sample-by-sample. Therefore I have to filter the signal > >> > sample by sample (Not the whole signal or block by block). I have > >> > implemented a filter to filter the whole signal. But Im struggling to > >> > modify it make it work sample by sample basis. > >> > > >> > I highly appreciate if someone can advice me regarding this issue. I > >> > have included the code Im using to filter. > >> > > >> > > >> > [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter > >> > coefficient generation (Cutoff frequency is 2*fc) > >> > > >> > x=filter(b,a,y); % y is the input signal. > >> > > >> > % Here these two lines do the perfect filtering operation for the > >> > entire signal, But I want to process sample by sample (In-side a > >> > phase locked loop) > >> > > >> > Can someone please advice how to do that using available matlab > >> > function or different approach? > >> > >> 1: I'm not sure what order Matlab uses for it's a and b vectors. I'm > >> fairly sure that a(1) corresponds to the z^0 coefficient, a(2) > >> corresponds to z^1, etc. -- but you need to check this. > >> > >> 2: You absolutely, positively, do not want to implement a 10th order > >> IIR filter as a single direct-form filter. That's Very Bad -- you'll > >> get no end of problems with coefficient rounding. In fact, I'm > >> surprised you got sensible results in batch mode, and I strongly > >> suggest that you check the Bode plot of the filter that you build with > >> those A and B vectors. > >> > >> 2a: Split the filter up into five 2nd-order sections and cascade them. > >> > >> 3: I'm not even sure if your butter and filter functions are working in > >> the Laplace domain or the z domain -- check. If they're in the Laplace > >> domain (i.e., continuous time filtering), then check back here for more > >> guidance. > >> > >> 4: If you are implementing this 10th order filter inside your control > >> loop, just stop. You're doing something wrong. If you make the filter > >> cutoff frequency low enough to be useful you will get so much phase > >> shift that you'll kill all stability. If you make the filter cutoff > >> high enough to maintain stability, your filter will be useless. > >> > >> The vectors a and b express a difference equation (if butter is > >> designing a sampled-time filter): > >> > >> x(n) = y(n) * b(11) + y(n-1) + b(10) + ... + y(n-10) * b(1) + > >> -x(n-1) * a(10) - x(n-2) * a(9) - ... - x(n-10) * a(1) > >> > >> This is a direct-form (type mumble mumble -- I can never keep them > >> straight, but it runs from 1 to 4, or 1 and 2 and 1 modified and 2 > >> modified -- it depends on the author). > >> > >> But you don't want to realize the above filter: you want to factor your > >> a and b polynomials, assemble each one back into five 2nd-order > >> polynomials, choose the b polynomials whose roots are closest to the a > >> polynomials, and implement the filter whose transfer function is > >> > >> B/A = B1/A1 * B2/A2 * ... * B5/A5. > >> > >> In the time-domain, that means that you cascade filters 1 through 5 > >> (i.e., run filter 1's output into filter 2, thence into filter 3, > >> etc.). > >> > >> For each filter, just implement the above difference equation, only > >> with delays that never exceed 2. > >> > >> -- > >> Tim Wescott Control system and signal processing consulting > >> www.wescottdesign.com > > > > Hi Tim, > > > > Thanks very much for your reply. > > > > I didn't look at the phase shift of the filter. But when I filter a > > signal and viewed the frequency spectrum it has filtered the higher > > frequency component. I choose the order as 10 because when I use low > > order like 2 it didnt remove the high frequency component completely. > > > > I didn't look at the body plots because the frequency response of the > > filter (freqz(b,a)) seems to be a non oscillating graph. > > > > I am not familiar with digital signal processing implementations. > > Therefore inorder to implement a required filter (to remove high > > frequency component 2*fc and take DC component from a sinusoidal > > multiplied signal) what should be the process. > > > > I highly appreciate your advice. > > > > Thank you. > > Well, are you making the filter part of your loop filter, or are you > filtering the phase error and using it elsewhere? That's what's > important to the loop stability. > > A notch filter at 2*fc would do what you are asking for -- but is what > you are asking for what you need? > > What are you trying to _do_? > > -- > Tim Wescott > Control system and signal processing consulting > www.wescottdesign.comHi Tim, I am trying to implement a low pass filter in my loop filter. I want to remove higher carrier and filter the phase error to feed to VCO. Initially I used a FIR filter with 200 coefficients but it failed at noisy environment (SNR of 10dB). Tim, can I have your email address? Thank you.
Reply by ●June 17, 20122012-06-17
On 6/16/12 10:35 PM, Tim Wescott wrote:> On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote: >> >> I want to implement an butterworth low pass filter to filter the DC >> component of a signal inside a PLL - phase locked loop. I am processing >> sample-by-sample. Therefore I have to filter the signal sample by sample >> (Not the whole signal or block by block). I have implemented a filter to >> filter the whole signal. But Im struggling to modify it make it work >> sample by sample basis. >> >> I highly appreciate if someone can advice me regarding this issue. I >> have included the code Im using to filter. >> >> >> [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter coefficient >> generation (Cutoff frequency is 2*fc) >> >> x=filter(b,a,y); % y is the input signal. >> >> % Here these two lines do the perfect filtering operation for the entire >> signal, But I want to process sample by sample (In-side a phase locked >> loop) >> >> Can someone please advice how to do that using available matlab function >> or different approach? > > 1: I'm not sure what order Matlab uses for it's a and b vectors.they do it fucking backwards, Tim. just like they do their polynomials ass-backwards. and it's related to the fact that DC is at X(1) in matlab. how is it that Cleve or whoever decided that the first (since they don't have zeroth) coefficient should go to x^N. but if they had the ability to have 0 (or even negative) indices, they could have associated a(0) with x^0 . what an oversight. a shame it got carved into stone.> I'm > fairly sure that a(1) corresponds to the z^0 coefficient, a(2) > corresponds to z^1, etc. -- but you need to check this.no, Tim. it's fucked up. you have to do this backwards, and you have to subtract 1 from the index after using max() in an FFT, if you plan on doing any math to the frequency of that bin.> > 2: You absolutely, positively, do not want to implement a 10th order IIR > filter as a single direct-form filter. That's Very Bad -- you'll get no > end of problems with coefficient rounding.it might even go boom. (roots of coef set might wander spuriously into a bad neighborhood.)> In fact, I'm surprised you > got sensible results in batch mode, and I strongly suggest that you check > the Bode plot of the filter that you build with those A and B vectors. > > 2a: Split the filter up into five 2nd-order sections and cascade them. > > 3: I'm not even sure if your butter and filter functions are working in > the Laplace domain or the z domain -- check. If they're in the Laplace > domain (i.e., continuous time filtering), then check back here for more > guidance.it normally does it in the digital domain. and i think it gives results that fit right into the filter() function. while it's consistent with itself (to some extent ...), MATLAB still defines it backwards. (because if you're doing it backwards, you'll have trouble remaining self-consistent in every situation and you'll be needing the fliplr() function because they defined it wrong.)> 4: If you are implementing this 10th order filter inside your control > loop, just stop. You're doing something wrong.i dunno, Tim, a 10th order system in a closed loop? what could go wrong? ship it. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Reply by ●June 17, 20122012-06-17
On Sun, 17 Jun 2012 01:29:55 -0400, robert bristow-johnson wrote:> On 6/16/12 10:35 PM, Tim Wescott wrote: >> On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote: >>> >>> I want to implement an butterworth low pass filter to filter the DC >>> component of a signal inside a PLL - phase locked loop. I am >>> processing sample-by-sample. Therefore I have to filter the signal >>> sample by sample (Not the whole signal or block by block). I have >>> implemented a filter to filter the whole signal. But Im struggling to >>> modify it make it work sample by sample basis. >>> >>> I highly appreciate if someone can advice me regarding this issue. I >>> have included the code Im using to filter. >>> >>> >>> [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter >>> coefficient generation (Cutoff frequency is 2*fc) >>> >>> x=filter(b,a,y); % y is the input signal. >>> >>> % Here these two lines do the perfect filtering operation for the >>> entire signal, But I want to process sample by sample (In-side a phase >>> locked loop) >>> >>> Can someone please advice how to do that using available matlab >>> function or different approach? >> >> 1: I'm not sure what order Matlab uses for it's a and b vectors. > > they do it fucking backwards, Tim. just like they do their polynomials > > ass-backwards. > > and it's related to the fact that DC is at X(1) in matlab. > > how is it that Cleve or whoever decided that the first (since they don't > have zeroth) coefficient should go to x^N. but if they had the ability > to have 0 (or even negative) indices, they could have associated a(0) > with x^0 . > > what an oversight. a shame it got carved into stone. > >> I'm >> fairly sure that a(1) corresponds to the z^0 coefficient, a(2) >> corresponds to z^1, etc. -- but you need to check this. > > no, Tim. it's fucked up. you have to do this backwards, and you have > to subtract 1 from the index after using max() in an FFT, if you plan on > doing any math to the frequency of that bin. > > >> 2: You absolutely, positively, do not want to implement a 10th order >> IIR filter as a single direct-form filter. That's Very Bad -- you'll >> get no end of problems with coefficient rounding. > > it might even go boom. (roots of coef set might wander spuriously into > a bad neighborhood.) > >> In fact, I'm surprised you >> got sensible results in batch mode, and I strongly suggest that you >> check the Bode plot of the filter that you build with those A and B >> vectors. >> >> 2a: Split the filter up into five 2nd-order sections and cascade them. >> >> 3: I'm not even sure if your butter and filter functions are working in >> the Laplace domain or the z domain -- check. If they're in the Laplace >> domain (i.e., continuous time filtering), then check back here for more >> guidance. > > it normally does it in the digital domain. and i think it gives results > that fit right into the filter() function. while it's consistent with > itself (to some extent ...), MATLAB still defines it backwards. (because > if you're doing it backwards, you'll have trouble remaining > self-consistent in every situation and you'll be needing the fliplr() > function because they defined it wrong.) > >> 4: If you are implementing this 10th order filter inside your control >> loop, just stop. You're doing something wrong. > > i dunno, Tim, a 10th order system in a closed loop? what could go > wrong? ship it.Boy, I get this dim notion that you might be a little irate at the folks at the MathWorks, Robert. The way to define a transfer function in Scilab is H = (b2 * %z^2 + b1 * %z + b0)/(%z^2 + a1*%z + a2), so "backwards" doesn't even enter into it. And if you want Laplace domain, you say H = (b1 * %s + b0) / (%s + a0) etc. At any rate, _part_ of the reason I use Scilab is because I'm a cheap bastard -- but part is because it just works better for the stuff I do all the time. -- 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 ●June 17, 20122012-06-17
On Sat, 16 Jun 2012 21:12:04 -0700, Dinuka wrote:> On Sunday, June 17, 2012 1:57:54 PM UTC+10, Tim Wescott wrote: >> On Sat, 16 Jun 2012 19:52:29 -0700, Dinuka wrote: >> >> > On Sunday, June 17, 2012 12:35:49 PM UTC+10, Tim Wescott wrote: >> >> On Sat, 16 Jun 2012 17:36:45 -0700, Dinuka wrote: >> >> >> >> > Hi all, >> >> > >> >> > I want to implement an butterworth low pass filter to filter the >> >> > DC component of a signal inside a PLL - phase locked loop. I am >> >> > processing sample-by-sample. Therefore I have to filter the signal >> >> > sample by sample (Not the whole signal or block by block). I have >> >> > implemented a filter to filter the whole signal. But Im struggling >> >> > to modify it make it work sample by sample basis. >> >> > >> >> > I highly appreciate if someone can advice me regarding this issue. >> >> > I have included the code Im using to filter. >> >> > >> >> > >> >> > [b,a] = butter(10, 4*fc/Fs); % 10th Order butterworth filter >> >> > coefficient generation (Cutoff frequency is 2*fc) >> >> > >> >> > x=filter(b,a,y); % y is the input signal. >> >> > >> >> > % Here these two lines do the perfect filtering operation for the >> >> > entire signal, But I want to process sample by sample (In-side a >> >> > phase locked loop) >> >> > >> >> > Can someone please advice how to do that using available matlab >> >> > function or different approach? >> >> >> >> 1: I'm not sure what order Matlab uses for it's a and b vectors. >> >> I'm fairly sure that a(1) corresponds to the z^0 coefficient, a(2) >> >> corresponds to z^1, etc. -- but you need to check this. >> >> >> >> 2: You absolutely, positively, do not want to implement a 10th order >> >> IIR filter as a single direct-form filter. That's Very Bad -- >> >> you'll get no end of problems with coefficient rounding. In fact, >> >> I'm surprised you got sensible results in batch mode, and I strongly >> >> suggest that you check the Bode plot of the filter that you build >> >> with those A and B vectors. >> >> >> >> 2a: Split the filter up into five 2nd-order sections and cascade >> >> them. >> >> >> >> 3: I'm not even sure if your butter and filter functions are working >> >> in the Laplace domain or the z domain -- check. If they're in the >> >> Laplace domain (i.e., continuous time filtering), then check back >> >> here for more guidance. >> >> >> >> 4: If you are implementing this 10th order filter inside your >> >> control loop, just stop. You're doing something wrong. If you make >> >> the filter cutoff frequency low enough to be useful you will get so >> >> much phase shift that you'll kill all stability. If you make the >> >> filter cutoff high enough to maintain stability, your filter will be >> >> useless. >> >> >> >> The vectors a and b express a difference equation (if butter is >> >> designing a sampled-time filter): >> >> >> >> x(n) = y(n) * b(11) + y(n-1) + b(10) + ... + y(n-10) * b(1) + >> >> -x(n-1) * a(10) - x(n-2) * a(9) - ... - x(n-10) * a(1) >> >> >> >> This is a direct-form (type mumble mumble -- I can never keep them >> >> straight, but it runs from 1 to 4, or 1 and 2 and 1 modified and 2 >> >> modified -- it depends on the author). >> >> >> >> But you don't want to realize the above filter: you want to factor >> >> your a and b polynomials, assemble each one back into five 2nd-order >> >> polynomials, choose the b polynomials whose roots are closest to the >> >> a polynomials, and implement the filter whose transfer function is >> >> >> >> B/A = B1/A1 * B2/A2 * ... * B5/A5. >> >> >> >> In the time-domain, that means that you cascade filters 1 through 5 >> >> (i.e., run filter 1's output into filter 2, thence into filter 3, >> >> etc.). >> >> >> >> For each filter, just implement the above difference equation, only >> >> with delays that never exceed 2. >> >> >> >> -- >> >> Tim Wescott Control system and signal processing consulting >> >> www.wescottdesign.com >> > >> > Hi Tim, >> > >> > Thanks very much for your reply. >> > >> > I didn't look at the phase shift of the filter. But when I filter a >> > signal and viewed the frequency spectrum it has filtered the higher >> > frequency component. I choose the order as 10 because when I use low >> > order like 2 it didnt remove the high frequency component completely. >> > >> > I didn't look at the body plots because the frequency response of the >> > filter (freqz(b,a)) seems to be a non oscillating graph. >> > >> > I am not familiar with digital signal processing implementations. >> > Therefore inorder to implement a required filter (to remove high >> > frequency component 2*fc and take DC component from a sinusoidal >> > multiplied signal) what should be the process. >> > >> > I highly appreciate your advice. >> > >> > Thank you. >> >> Well, are you making the filter part of your loop filter, or are you >> filtering the phase error and using it elsewhere? That's what's >> important to the loop stability. >> >> A notch filter at 2*fc would do what you are asking for -- but is what >> you are asking for what you need? >> >> What are you trying to _do_? >> >> -- >> Tim Wescott >> Control system and signal processing consulting www.wescottdesign.com > > Hi, > > I am trying to filter the phase error in my loop filter. Actually I have > used a FIR filter with 200 coefficients before to remove the unnecessary > frequency components. And the PLL worked when there is no noise. But > when I add noise (at SNR of 10 dB) the PLL fails to lock. So I assume > that the delay in the Loop filter(Low pass filter) might cause the > problem. So I decided to use a IIR filter with low length so that the > delay would be low. > > Im sorry for my ignorance, but I am struggling to make a suitable low > pass filter.I think that if your approach requires a 10th order filter inside your loop, then -- unless you've got a bizarre environment or bizarre specifications -- you have a problem with your approach. Generally if you've got noise problems, then the first thing you want to do is narrow your loop bandwidth by reducing your loop filter gains; if that solution gives you a bandwidth that is too small before you get the noise reduction you need, then you have to consider fancier methods. I may have mentioned this before, but you may want to get a book on PLL techniques. The one that I know and love is a bit obscure: "Phase Locked Loop Circuit Design" by Dan Wolavar. The other classic one is by Floyd Gardener, with a title similar to "Phase Locked Loop Design". Both of these books don't address digital PLL techniques at all -- but the underlying theory is the same, regardless of whether you're implementing your loop filter in op-amps or lines of code.> Tim, can I have your email address?tim at wescott design dot com (remove spaces and make the obvious changes to "at" and "dot"). But please keep in mind that with very few exceptions I do not give away free advise over email or by phone -- that's what USENET is for. -- 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 ●June 17, 20122012-06-17
On 6/17/12 3:50 PM, Tim Wescott wrote:> On Sun, 17 Jun 2012 01:29:55 -0400, robert bristow-johnson wrote: > >> On 6/16/12 10:35 PM, Tim Wescott wrote:...>>> 1: I'm not sure what order Matlab uses for it's a and b vectors. >> >> they do it fucking backwards, Tim. just like they do their polynomials >> >> ass-backwards. >> >> and it's related to the fact that DC is at X(1) in matlab. >> > > Boy, I get this dim notion that you might be a little irate at the folks > at the MathWorks, Robert.ya think. back in 1995 or so, i had a phone conversation with Cleve Moler and subsequently some email exchanges, and a couple of exchanges on comp.soft-sys.matlab. while he listened and the exchange is polite, i came away from it that TMW is quite self-satisfied regarding their decisions regarding indexing.> > The way to define a transfer function in Scilab is > > H = (b2 * %z^2 + b1 * %z + b0)/(%z^2 + a1*%z + a2), > > so "backwards" doesn't even enter into it.we, that's a biquad. but if you were to define a general Nth-order IIR function with two vectors, one for the numerator coefs and one for the denominator (that you pass to a routine, which i see you're not doing with Scilab), what is the logical way to do it? and there *is* a frontwards and backwards to it. what is the logical order to specifying the coefficients to a general Nth-order polynomial? about the other old issue, what is the logical index for the DC bin coming back outa the FFT? and what is the logical manner to represent a non-causal FIR filter? the original authors of the very first MATLAB toolbox (the signal processing toolbox) were Thomas Krauss and Loren Shure. i had a brief email exchange with Mr. Krauss in the 90s (after he left TMW) and he told me that before the toolbox came to market, he told the TMW leadership that they should do something about the indexing issue (because DC doesn't have a frequency of 1 in *any* frequency unit) and this advise fell on deaf ears. i met Ms. Shure at a workshop on wavelets and filter banks (i went to this with Randy Yates) and she just sang the company line (that extending the indices to non-positive integers would "break" MATLAB, which just is not true). in our email conversations and on the newsgroup (in the 90s) Cleve also said that it wasn't backward compatible and then i took a great deal of time to spell out exactly how the MATLAB "variable" would be extended to accommodate non-positive indices and how it would be made backward compatible and how various common operations (like multiplication of arrays) would work with it. deaf. TMW is big, monolithic, and quite self-satisfied that they just cannot imagine that they're doing something fundamentally flawed. and, in this case, fundamental flaws lead to bad conventions and bad conventions lead to unnecessary and hard-to-find errors.> And if you want Laplace > domain, you say > > H = (b1 * %s + b0) / (%s + a0) > > etc. > > At any rate, _part_ of the reason I use Scilab is because I'm a cheap > bastard -- but part is because it just works better for the stuff I do > all the time.i can believe that. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
