Re: laplace tranform convert to code

tim w wrote:

> Tim Wescott <tim@seemywebsite.com> wrote in
> news:iv-dncPJu77Y1VzZnZ2dnUVZ_oadnZ2d@web-ster.com: 
> 
> 
>>tim w wrote:
>>
>>>Hello all,
>>>
>>>I need some guidance in programming a laplace transfer function into 
>>>computer language -- pseudocode for now.
>>>
>>>The transfer function is a second order function:
>>>
>>>To^2*s^2 + zeta1*To*s + 1
>>>-------------------------
>>>To^2*s^2 + zeta2*To*s + 1
>>>
>>>From what i've read the above transfer function is a bandstop filter,
>>>Zeta1 and Zeta2 are adjustable paramaters. A demand signal is
>>>conditioned by the function.
>>>
>>>What are the steps I need to take to convert to software? Should I
>>>convert to z? Or??? 
>>>
>>>Any web site or recomended book??
>>>
>>>by the way, i've never programmed from transfer functions so I am new
>>>at this.
>>>
>>>thanks
>>
>>This is a perfect question for the comp.dsp group; I am taking the 
>>liberty of cross-posting it there.
>>
>>Yes, if you want this to actually execute in the digital domain you'll
>>need to convert to the z domain one way or another, then write
>>software to the resulting transfer function.  If you must start from
>>the s domain and go to the z domain then your best bet is to use
>>Tustin's approximation ("bilinear transform") after prewarping the
>>poles.  When I can I prefer to start with the requirements for the
>>filter and just design the whole darn thing in the z domain -- this is
>>particularly important (albeit frustrating) for a notch filter where
>>you want to make sure the gain above the notch is the same as the gain
>>below. 
>>
>>"Understanding Digital Signal Processing" by Rick Lyons will get you a
>>long way down the road.  It includes approximating Laplace-domain 
>>filters in the z domain, but skimming through the table of contents
>>and flipping through the book I don't see anything that looks like a
>>promise of code (Rick?).
>>
>>My book, "Applied Control Theory for Embedded Systems" gives you the 
>>tools you need to go from a z-domain transfer function to code, but 
>>it'll be pretty light in getting from the s-domain to the z-domain --
>>I take refuge in claims of finite page counts and finite time.
>>
>>I hope this helps.
>>
> 
> 
> What a coincidence Tim.
> 
>  
> 
> Thanks fro crossposting on the other group, comp.dsp.group. I wandered
> over and found some interesting stuff. I came accross this web site
> which mentions that one can convert a laplace transfer function into the
> z plane by taking the bilinear transform. So I did and this is what I
> have, I just need to see how to convert or group together. Here it goes
> but before looking, take a look at the following web site that I used as
> a guide for the transfer function I want to convert. 
> 
>  
> 
> http://www.earlevel.com/Digital%20Audio/Bilinear.html
>  
> 
> So here is my transfer function.
> 
>  
> To^2*s^2 + zeta1*To*s + 1
> -------------------------
> To^2*s^2 + zeta2*To*s + 1
> 
> 
> I substituded the s into ((1/K * (z-1)/(z+1)). I also simplified and
> collected the z terms and this is what I ended up with: 
> 
> 
> (K^2+To^2+2*zeta1*To*K)*z^2 + (2*K^2-2*To^2)*z - 2*zeta1*To*K+K^2 +To^2
> -------------------------------------------------------------------------
> (K^2+To^2+2*zeta2*To*K)*z^2 + (2*K^2-2*To^2)*z - 2*zeta2*To*K+K^2 +To^2 
> 
> 
> I then multiplied by z^-2 and came up with the following:
>  
> 
> (K^2+To^2+2*zeta1*To*K)+(2*K^2-2*To^2)*z^-1+(-2*zeta1*To*K+K^2+To^2)*z^-2 
> -------------------------------------------------------------------------
> (K^2+To^2+2*zeta2*To*K)+(2*K^2-2*To^2)*z^-1+(-2*zeta2*To*K+K^2+To^2)*z^-2 
> 
> 
> So I collect the terms for a0, a1, a2, b0, b1, and b2. Now, in some
> other websites, they have the b0 to b2 terms in the numerator. In Earl's
> page, the have them on the numerator so that is what I am using. One
> last thing I did and that was normalize by dividing each by b0. 
> 
> 
> a0 = (K^2 + To^2 + 2*zeta1*To*K) / (K^2 + To^2 + 2*zeta2*To*K)    
> 
> a1 = (2*K^2 - 2*To^2) / (K^2 + To^2 + 2*zeta2*To*K)
> 
> a2 = (-2*zeta1*To*K + K^2 + To^2) / (K^2 + To^2 + 2*zeta2*To*K)
> 
> 
> b0 = 1
> b1 = a1
> b2 = a2
> 
> K = tan (pie (fc/fs))
> fs = sample rate
> 
> 
> Questions:
> 
> - Did I did okay in following the example in Earl's page?

It appears so, but you'll have to excuse me for not going through the 
math in detail -- any time you do something like this you should test 
it, which I suggest you do at least by doing a frequency response plot 
of the z-domain transfer function.

There's an article on my web site, 
http://www.wescottdesign.com/articles/zTransform/z-transforms.html, that 
covers this.  The information is also (I think) in Rick's book, and 
certainly in mine.

> - Does it seem like I took the right steps?

It appears so.

> - What do I now do with the above equations? May somebody post or
> provide a link for the difference equation I need to use to get the
> above in code??? 

There is an infinite number of ways to do this.  The "Direct Form 1" way is:

            Y(z)   b_2 z^2 + b_1 z + b_0
For H(z) = ---- = ---------------------
            U(z)     z^2 + a_1 z + a_0

you get

y(n) = b_2 u(n) + b_1 u(n-1) + b_2 u(n-2) - a_1 y(n-1) - a_0 y(n-2).

Rick's book goes into this in great detail, and you can probably do a 
web search on "direct form" and "biquad" to get some discussions about 
the most popular architectures.

> 
> -Also, the variable "fc" is the corner frequency. What is that on how do
> I set it??? 
> 
In a band stop filter fc is the center frequency of the stop band.  In 
your transfer function everything is scaled in time, with "To".  For fc 
in Hz, fc = 1/(2 * pi * To).
> 
> Next step would be test but I need to but the a's and b's together
> before continuing. 
> 
> Any help/guidance is appreciated.


-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Posting from Google?  See http://cfaj.freeshell.org/google/

"Applied Control Theory for Embedded Systems" came out in April.
See details at http://www.wescottdesign.com/actfes/actfes.html