laplace tranform convert to code

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

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.

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

Tim Wescott wrote:

   ...


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

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

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

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

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

Looks like your answering your own questions there timmie

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

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.

So would a shameless self promoter ;-)