Re: laplace tranform convert to code

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?
- Does it seem like I took the right steps?
- 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??? 

-Also, the variable "fc" is the corner frequency. What is that on how do
I set it??? 


Next step would be test but I need to but the a's and b's together
before continuing. 

Any help/guidance is appreciated.