laplace tranform convert to code

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
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
tim w wrote:

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

Tim Wescott
Wescott Design Services
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

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.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
tim w wrote:

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.

It appears so.

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.

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

--

Tim Wescott
Wescott Design Services
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Tim Wescott wrote:

...

Other tim,
Remember to prewarp the critical frequencies. Rick's book goes into it, and so do many other sources. Robert Bristow-Johnson wrote a clear description (in comp.dsp, I think); Google should turn it up. If not, Tim or I can describe the process, the equation is simple.
Jerry
--
Engineering is the art of making what you want from things you can get.

  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Jerry Avins wrote:

Check his original URL: The author has a clever way of building the prewarping into the bilinear transform. Dunno if I'd want to make a habit of using it, but it's cute.
--

Tim Wescott
Wescott Design Services
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Tim Wescott wrote:
...

If w0 is half way between the edges of a bandpass filter in the s domain, it won't be after bilinear transformation. The simplest way I know is too specify the band edges separately, not band center and bandwidth. How does the prewarped transform handle that?o
Jerry
--
Engineering is the art of making what you want from things you can get.

  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Jerry Avins wrote:

Dunno -- I've never trusted it enough to try. I always make my notch filters by direct design in the z domain.
--

Tim Wescott
Wescott Design Services
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

That is similar to a notch filter. I posted a link to a notch filter design not too long ago ftp://ftp.deltacompsys.com/public/PDF/Mathcad%20-%20Notch.pdf Jerry didn't believe it worked so I changed the frequency ftp://ftp.deltacompsys.com/public/PDF/Mathcad%20-%20Notch49Hz.pdf
You should be able to take this example a change it to fit your needs.
Peter Nachtwey
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Looks like your answering your own questions there timmie

Zeta1
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload
Setanta wrote:

-- snip --
Interestingly enough, there are more than one person in the world named 'Tim', and more than one with last names ending in 'W'.
A good statistician could probably even calculate approximately how many.
--

Tim Wescott
Wescott Design Services
  Click to see the full signature.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

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.
Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

So would a shameless self promoter ;-)

Add pictures here
<% if( /^image/.test(type) ){ %>
<% } %>
<%-name%>
Add image file
Upload

Polytechforum.com is a website by engineers for engineers. It is not affiliated with any of manufacturers or vendors discussed here. All logos and trade names are the property of their respective owners.