Forums

Re: laplace tranform convert to code

Started by Tim Wescott July 24, 2006
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