Thread Subject:

Matlba butterworth IIR filter implementation - sample by sample

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."

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."

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...

--

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."

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."

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? :)

--

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."

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

  

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...

--