# Thread Subject: Matlba butterworth IIR filter implementation - sample by sample

 Subject: Matlba butterworth IIR filter implementation - sample by sample Date: 17 Jun, 2012 05:29:55 Message: 1 of 9 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."
 Subject: Matlba butterworth IIR filter implementation - sample by sample From: dpb Date: 17 Jun, 2012 20:31:32 Message: 3 of 9 On 6/17/2012 12:29 AM, robert bristow-johnson wrote: ... > 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 . ... This thread just appeared in mid-stream at cs-sm... Just a comment on the above diatribe-- It's the normal way in which mathematics writes a polynomial...consider the time-honored quadratic formula for y= ax^2 + bx + c polyval() and friends are consistent with that time-honored formulation. That it's convenient depending on the problem domain to reverse the order is another issue altogether; one can't infer Matlab is somehow [expletive expression deleted]. It is indeed a pity that Matlab didn't include the facility for choosing the origin from the git-go. However, since it evolved from even earlier work on LINPACK and was initially coded in the 70s before F77, the native support in a Fortran compiler of the time wasn't around (or if was, was a vendor extension) so it's definitely unfair to point fingers at this point on what was then. The point at which it _could_ have happened was the time at which it was rewritten in C... --
 Subject: Matlba butterworth IIR filter implementation - sample by sample Date: 17 Jun, 2012 22:11:58 Message: 4 of 9 On 6/17/12 4:52 PM, Tim Wescott wrote: > > Well, if you feel bound and determined to do it with vectors, here's the > example (with Scilab's response) > > -->H = syslin('d', poly((0.0001/3) * [1 1 1], 'z', 'coeff') / > poly ([0.9859 -1.9858 1], 'z', 'coeff')) > H = > > 2 > 0.0000333 + 0.0000333z + 0.0000333z > ----------------------------------- > 2 > 0.9859 - 1.9858z + z > well, i can't tell with the numerator, but looking at the denominator it is apparent that Scilab differs from MATLAB regarding convention of ordering the coefficients of a polynomial. MATLAb would say it's [1 -1.9858 0.9859]. which is, of course, ass-backwards. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
 Subject: Matlba butterworth IIR filter implementation - sample by sample Date: 17 Jun, 2012 22:43:45 Message: 5 of 9 On 6/17/12 4:31 PM, dpb wrote: > On 6/17/2012 12:29 AM, robert bristow-johnson wrote: > ... > >> 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 . > ... > > This thread just appeared in mid-stream at cs-sm... > > Just a comment on the above diatribe-- > > It's the normal way in which mathematics writes a polynomial...consider > the time-honored quadratic formula for > > y = ax^2 + bx + c > > polyval() and friends are consistent with that time-honored formulation. so this is the "time-honored" representation of a power series (i do ASCII math with mono-spaced plain text on USENET):            N+1     y = SUM{ a_(N+1-n) * x^(N-1+n) }            n-1 my that looks clean. certainly cleaner than:             N     y = SUM{ a_n * x^n }            n=0 and, for the DFT, i'm sure this is time honored:                     N     x(n) = 1/N SUM{ X(k) * e^(i*2*pi*(n-1)*(k-1)/N) }                    k=1 so time-honored that new DSP textbooks are now starting to change their definitions in print, just so that the book can be consistent with MATLAB. DC has a frequency of 1. thanks to TMW, we'll all be directed back to the proper "time-honored" convention of counting from one. indices can have not other meaning than that. and i'm sure Edsger Dijkstra is happy (wherever he is now) about that: http://www.cs.utexas.edu/users/EWD/transcriptions/EWD08xx/EWD831.html > That it's convenient depending on the problem domain to reverse the > order is another issue altogether; one can't infer Matlab is somehow > [expletive expression deleted]. the purpose of having a good convention, rather than a bad convention (whether some have honored it over time or not) is so that you don't have to remember illogical details, just because someone had insisted on such illogical details. > It is indeed a pity that Matlab didn't include the facility for choosing > the origin from the git-go. i'm glad to read someone outside of comp.dsp saying so. but it wasn't to late to fix this in the 90s, and it isn't even too late now, unless there is little future to MATLAB. it can be easily fixed, and be perfectly backward compatible, because new arrays that are declared will continue to have their upper left corner at (1,1). > However, since it evolved from even earlier > work on LINPACK and was initially coded in the 70s before F77, the > native support in a Fortran compiler of the time wasn't around (or if > was, was a vendor extension) so it's definitely unfair to point fingers > at this point on what was then. i disagree. i've written Fortran IV, and i hated the indexing, the Hollerith format, etc. but i wrote user-friendly code (by adjusting from bad conventions to good, in my code). it's what we have to do now, writing something in MATLAB that other people are intended to use. so *i* need to remember to subtract 1 from indices so that users don't need to and *i* need to take the coefficient vector, ordered naturally, and fliplr() the damn thing, so that other users don't need to. but the purpose of MATLAB is so that *don't* have to fiddle around like that, so that equations we type into the program resemble, as much as possible, the compact equations we write down on paper. MATLAB v4 (in the 90s) used to claim you could do that in their paper manual. of course it's not true. > > The point at which it _could_ have happened was the time at which it was > rewritten in C... as it evolved from LINPACK to a commercial product. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
 Subject: Matlba butterworth IIR filter implementation - sample by sample From: dpb Date: 17 Jun, 2012 23:39:05 Message: 6 of 9 On 6/17/2012 5:43 PM, robert bristow-johnson wrote: > On 6/17/12 4:31 PM, dpb wrote: ...[earlier portions of book elided for brevity and lack of interest in debating]... > > i disagree. i've written Fortran IV, and i hated the indexing, the > Hollerith format, etc. but i wrote user-friendly code (by adjusting from > bad conventions to good, in my code). it's what we have to do now, > writing something in MATLAB that other people are intended to use. so > *i* need to remember to subtract 1 from indices so that users don't need > to and *i* need to take the coefficient vector, ordered naturally, and > fliplr() the damn thing, so that other users don't need to. > > but the purpose of MATLAB is so that *don't* have to fiddle around like > that, so that equations we type into the program resemble, as much as > possible, the compact equations we write down on paper. MATLAB v4 (in > the 90s) used to claim you could do that in their paper manual. of > course it's not true. > ... Well, I've done the same thing as well. Unfortunately, at this point it seems to be tilting at windmills as far as Matlab goes, for good or ill. If it bothers you sufficiently I guess you have a couple of choices--don't use Matlab or apply the same concept and write a set of wrappers around the implementations to deal with the order and offset so you can refer to them in the way you would like. The "write equations as paper" is reasonably true for a given class of matrix algebra problems--it just so happens those aren't necessarily the same class of problems as found in much of DSP. Changing conventions for that purpose breaks the convention for others--who's to say it should be for your particular set of problems rather than somebody else's? :) --
 Subject: Matlba butterworth IIR filter implementation - sample by sample Date: 18 Jun, 2012 00:06:51 Message: 7 of 9 On 6/17/12 7:39 PM, dpb wrote: > ... > The "write equations as paper" is reasonably true for a given class of > matrix algebra problems-- a given class. there is more to applied mathematics than this given class. > it just so happens those aren't necessarily the > same class of problems as found in much of DSP. oh? take a look at IEEE transactions. nobody seems to publish there without some attempt (usually it's an effort to over-mathematicalize to impress reviewers and other readers) to express the idea in matrix form > Changing conventions for > that purpose breaks the convention for others--who's to say it should be > for your particular set of problems rather than somebody else's? what i have been advocating, from square one (or square "zero" for the rest of us) is not *changing* the convention. it is *broadening* the convention to allow users to change what the origin index is for each dimension. and since any array (or "matrix", if you prefer) is instantiated with the origin set to 1, it does *not* break it for anyone who doesn't make use of this extension. Cleve and others tried that argument on me better than a decade ago and it didn't fly then and it doesn't fly now. you probably have to look this up in google groups to see what the particulars were (it really isn't all that complicated, but refuting objections is not effortless). i wrote tens of thousands of words and i'm not going to do it again, since the minds at TMW are closed. maybe i'll try to find you a link, if it interests you. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
 Subject: Matlba butterworth IIR filter implementation - sample by sample From: Nasser M. Abbasi Date: 18 Jun, 2012 02:47:39 Message: 8 of 9 On 6/17/2012 7:06 PM, robert bristow-johnson wrote: > >i wrote tens of thousands of words and > i'm not going to do it again, since the minds at TMW are closed. maybe > i'll try to find you a link, if it interests you. > here is one I remember http://www.mathworks.com/matlabcentral/newsreader/view_thread/285566 --Nasser
 Subject: Matlba butterworth IIR filter implementation - sample by sample From: dpb Date: 18 Jun, 2012 03:05:03 Message: 9 of 9 On 6/17/2012 7:06 PM, robert bristow-johnson wrote: > On 6/17/12 7:39 PM, dpb wrote: >> > .... >> The "write equations as paper" is reasonably true for a given class of >> matrix algebra problems-- > > a given class. there is more to applied mathematics than this given class. > >> it just so happens those aren't necessarily the >> same class of problems as found in much of DSP. > > oh? take a look at IEEE transactions. nobody seems to publish there > without some attempt (usually it's an effort to over-mathematicalize to > impress reviewers and other readers) to express the idea in matrix form Well, there's a _very_ large literature outside of IEEE Transactions as well for whatever that's worth...I don't see how that has a bearing on the complaint you have for a particular set of DSP >> Changing conventions for >> that purpose breaks the convention for others--who's to say it should be >> for your particular set of problems rather than somebody else's? > > what i have been advocating, from square one (or square "zero" for the > rest of us) is not *changing* the convention. it is *broadening* the > convention to allow users to change what the origin index is for each > dimension. and since any array (or "matrix", if you prefer) is > instantiated with the origin set to 1, it does *not* break it for anyone > who doesn't make use of this extension. Cleve and others tried that > argument on me better than a decade ago and it didn't fly then and it > doesn't fly now. > > you probably have to look this up in google groups to see what the > particulars were (it really isn't all that complicated, but refuting > objections is not effortless). i wrote tens of thousands of words and > i'm not going to do it again, since the minds at TMW are closed. maybe > i'll try to find you a link, if it interests you. I'll grant there's possibly a way to extend syntax to include other-than-unity-origin vectors. I'd think the only way TMW will consider it will be in the normal sequence of submitting enhancement request through the official support channel. I wouldn't hold my breath on "real soon now", though... :) You'll have more results quicker from one of the other choices I outlined earlier I expect unless you do just enjoy the battle, though. I sympathize w/ the lack; just not expecting it to change means I've learned to adjust. It's something like the new users who are used to that other language that want to try to turn Matlab into it as well--Matlab is what it is so "when in Rome..." applies. I've gotten more involved here than should have already... --

Separated by commas
Ex.: root locus, bode

### What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.