Matlba butterworth IIR filter implementation - sample by sample filtering basis

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 Sun, 17 Jun 2012 16:17:45 -0400, robert bristow-johnson wrote:

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

-->scf(0); clf; bode(H, 1.0001e-5, 0.4999, 0.01)

But isn't it easier to just use the native transfer function form that 
Scilab knows and understands?

-->H = syslin('d', (0.0001 / 3) * (%z^2 + %z + 1) / (%z^2 - 1.9858 * %z + 
0.9859))
 H  =
 
                                       2  
    0.0000333 + 0.0000333z + 0.0000333z   
    -----------------------------------   
                               2          
           0.9859 - 1.9858z + z           
 
-->scf(0); clf; bode(H, 1.0001e-5, 0.4999, 0.01)

The resulting transfer functions -- and yes, H is a _transfer function_ 
in Scilab, as a _built in native type_ -- can be run through a filter 
routine, too:

-->scf(1); clf; h = flts(ones(1, 1000), H); plot2d(h)

(Note that I'm plotting these -- that's to make sure that I'm saying 
sensible things so no one says "but that's an unstable filter, you 
doof!".  I trust you can see the underlying functions, though).

Note, too, that you don't have to stop at transfer functions:  you can 
define a linear system in state-space form:

Hss = syslin('d', [1.9858 -0.9953; 0.9906 0], [-0.01; 0.00005], [-0.01 
0]);

(I kinda messed this up with rounding, but you get the general idea).

With the state-space form, you can make a bunch of cascaded 2nd-order 
filters that executes (or Bode plots) all in one thump, and from the 
outside you manipulate the system just the same as if it were in transfer 
function form -- in fact, the answer to H * Hss is the filter, in state-
space form -- that you get from cascading H and Hss.

(If you haven't guessed, I am _not_ irate with the folks at Scilab!)

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

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

--

Hi,

regarding the original question:

>> But I want to process sample by sample 

check the "doc filter" page in Matlab:
http://www.mathworks.se/help/techdoc/ref/filter.html

What you need is
[y,zf] = filter(b,a,X) returns the final conditions

and

[y,zf] = filter(b,a,X,zi) accepts initial conditions,

Then run it on a single-sample input for each round, save and restore the
initial conditions.

A FIR filter in a PLL might be problematic unless you design it explicitly
as nonlinear phase (asymmetric impulse response), which is a FIR filter
pretending very hard to be an IIR filter, which usually doesn't make too
much sense. 
Try the run-of-the-mill "[b, a] = butter(3, 0.5)" filter (for example),
plug the [b, a] into the filter() function and make it work. 

On Sun, 17 Jun 2012 18:11:58 -0400, robert bristow-johnson wrote:

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

Unless the polynomial is in z^-1.

-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com