Does anyone have experience with an alternative for the bilinear transform?

[FLWrd]

Hi,

The other day, I stumbled across this thread (http://dsp.stackexchange.com/questions/ ... -transform, bottom reply), where the poster remarks that the traditional s <- 2/T (z-1) / (z+1) is actually trapezoidal integration, and that s <- 3/T (z^2 - 1) / (z^2+4z+1) gives Simpom's rule. The derivation seems valid. However, when I try it, I get unstable filter coefficients.

E.g., an RC filter has a transfer function H(s) = 1/(1+RCs), right? The traditional substitution (R=100 Ohm, C=1uF) leads to a0 = a1 = T / (0.00032+T) }
b1 = (1-3125*T) / (-1-3125*T).

Using in the alternative substitition yields a1 = T / (0.00048+T), a2 = T / (0.00012+0.25*T) a3 = T / (0.00048+T), b1 = T / (0.00012+0.25*T), b2 = (1-2083.3333333333333*T) / (-1-2083.3333333333333*T), or {0.04511,0.18044,0.04511} and {1,0.18044,-0.90978} for T = 1/44100, which is unstable.

Does anyone know what is wrong?

[Music Engineer]

i didn't check the math, but chances are, that there's nothing wrong at all. simpson's rule and (AFAIK) also all the other common higher order numerical integration schemes do not in general preserve the stability of an analog prototype filter. the trapezoidal rule is kinda special in this regard.

[FLWrd]

Thanks. I thought I intuitively understood the idea of some equivalence between integration and transformation, but clearly not.

[DaveGamble]

Thanks. I thought I intuitively understood the idea of some equivalence between integration and transformation, but clearly not.

Oooh.. Well... there is... sorta...

I mean... Are you familiar with the phrase "conformal mapping"?

There IS a connection between integration and "transformation" (which is a term I'm afraid I can't permit you to use, sorry), better described as mapping, but to describe it as "intuitive" would indicate that you've not got it right.

The mathematical objects at play here are:
- The space of functions from N->C (and, of course, the space from C_[-pi,p]->C, which is the space of Fourier transforms of such functions)
- The space of functions from R->C (and hence C->C as Fourier transforms)
(you're welcome to argue that the target space is R, but I counterargue that that's an irrelevant simplification and wasting time.)

I'm going to be a leeeeeetle bit messy, because in reality we have much more than just these spaces; we have a C* algebra on the spaces (of which the set of characters form the basis for the "Fourier" domains, but really this is by-the-by), but because the nomenclature is comfortable, I'm going to claim Z and S for the N->C algebra and R->C algebra respectively.

Now, differentiation is an operator in S. It's usually denoted "s", and it implies a differentiation in the most comfortable abstract sense of an operator on a Hilbert space. Your typical calculus derivative will turn out to work just fine here almost all the time (and I say that in the technical sense, as you'd expect).

Now, your term "transformation" which I'm mechanically rewriting as "conformal mapping", is a specific example of the notion of a mapping from S to Z. You mean it to imply some kind of mapping from the line i*f in the Fourier domain of S to end up having something to do with the mapping on the unit circle e^i*f in the Fourier domain of Z.

So, perhaps what you're imagining is the simulation of the linear operator s in S by some function in Z.
s IS differentiation in a comfortable sense, and conformal mapping CAN often begin by attempting to build some approximation to s within Z's algebra.

But until you evaluate the strategy (assuming that it is, in fact, what you're using), you don't actually KNOW anything about your conformal map. One cannot simply "assume" when working in spaces of complex linear functionals. A little working is a very good thing.

Dave.

[xoxos]

i totally have no truck with the transfer function math. if it's worth anything at all to say, from my uninformed, utterly naive standpoint, the bilinear transform isn't too far away from euler's equation (at least the one i've met, i expect there's more..) which is experientially and purportedly unstable unless oversampled a few times.

below that point, it works like a charm. of course, the fact that the two equations only remind me of each other may not be much help.

[AXP]

There IS an equivalence between Trapezoidal integration and Biliner "transform" (mapping). You can see that just by fiddling a bit with corresponding expressions, or by realizing that Trapezoidal rule approximates a differential with finite difference at the midpont, while Bilinear transforms maps the "s" operator (differentiation in the S-domain) into the Z-domain ("z" = discrete delay) using the same mid-point technique. Sorry for the layman-ish terms.

You can of course approximate "s" by higher order methods, but just as with numeric integration they may have better accuracy, but might be unstable. Bilinear transform does seem like a sweet spot in most cases. If you are unhappy with accuracy, oversample!
I've spent a little time experimenting in the same direction as you, but got results not unlike yours.

[DaveGamble]

Bilinear transform does seem like a sweet spot in most cases. If you are unhappy with accuracy, oversample!

^ This is what's wrong with everything in the world. This is why we can't have nice things.

Edit: After posting this, it dawned on me that after reading my reply, you still wouldn't know any better. So if I'm to prove of any use, I should explain the problem. Bilinear transform is a poor, poor strategy, and just googling Orfanidis will get you a better result. And if you care about such things (which actually OP didn't, best I can tell), it falls to you to develop good techniques. You can't just oversample, because if you do, you're just adding more problems. To start with, you had a poor response. Now you have a poor response and the need to downsample. So you're paying for bad design with latency, cpu cost, and a poor response from your downsampling filter. Oversampling isn't a silver bullet. Quite the opposite. It can easily leave you with the same, or worse problems than you started with. It's like rope to hang oneself with. With a BLT, you just have a shit EQ/filter. With an oversampled BLT, you now have a shit EQ/filter, latency, extra cpu cost, and the need to implement a respectable resampling algorithm. And I'll wager that if BLT appeals to you, your depth of insight into resampling could do with some work too.

[AXP]

Bilinear transform does seem like a sweet spot in most cases. If you are unhappy with accuracy, oversample!

^ This is what's wrong with everything in the world. This is why we can't have nice things.

Edit: After posting this, it dawned on me that after reading my reply, you still wouldn't know any better. So if I'm to prove of any use, I should explain the problem. Bilinear transform is a poor, poor strategy, and just googling Orfanidis will get you a better result. And if you care about such things (which actually OP didn't, best I can tell), it falls to you to develop good techniques. You can't just oversample, because if you do, you're just adding more problems. To start with, you had a poor response. Now you have a poor response and the need to downsample. So you're paying for bad design with latency, cpu cost, and a poor response from your downsampling filter. Oversampling isn't a silver bullet. Quite the opposite. It can easily leave you with the same, or worse problems than you started with. It's like rope to hang oneself with. With a BLT, you just have a shit EQ/filter. With an oversampled BLT, you now have a shit EQ/filter, latency, extra cpu cost, and the need to implement a respectable resampling algorithm.

Oh thanks for putting effort in writing that, but I'm aware of all the implications you mention and was just answering OP's initial question, who as you've said "didn't care about such things".

And I'll wager that if BLT appeals to you, your depth of insight into resampling could do with some work too.

I don't quite understand the offensive tone of your post, but well, there's never too much depth in any insight, can't argue that.

[DaveGamble]

I don't quite understand the offensive tone of your post, but well, there's never too much depth in any insight, can't argue that.

I am sorry if I've offended, but I felt quite strongly that it was bad advice given perhaps too carelessly.

Good to know you're on my side of the equation! I find I have to be very careful giving advice with this stuff!

Dave.

[FLWrd]

I'll sketch the background of my question. I'm interested in the effect a relatively small electronic circuit in a digital synth has on its sound. I've found the schematics, and have got a transfer function (which sort of seems to fit with the device's output). The circuit consists of a 3rd and 2nd order Sallen-Key filter in series and the second step has a peak followed by a steep roll-off at 15kHz, which is pretty high. And indeed, applying BLT results in a frequency response that differs quite a lot from the original. So I started to look for alternatives, and they all have their own problems. Then I found this, but just couldn't understand why it doesn't work.

I still don't really understand why it doesn't work, but I'm thankful that for the explanations. (I did hear about conformal mapping way back at university, but it was just one of those features of complex functions that didn't get any further treatment; that was for the physics/electronics guys, while I continued with language algebras and complexity; consequently, I couldn't tell if a mapping is conformal if it strips, paints itself pink and starts singing "conformal mappings are here again" while dancing on a harpsichord).

In the mean time, I've found more alternatives (MTZi), and I will definitely read up on Orfanidis. Thanks for that too.

which actually OP didn't, best I can tell

Well, let's hope you're better with maths than with words... (smiley of choice here)

[DaveGamble]

which actually OP didn't, best I can tell

Well, let's hope you're better with maths than with words... (smiley of choice here)

Your original post asked why a particular conformal mapping did not guarantee stability, and later about how this jars with your intuition.
You actually didn't ask about accuracy of match of magnitude response.

I wasn't being snarky or anything, you literally just weren't asking about that. I thought your question was very specific.

Dave.

[FLWrd]

I wasn't being snarky or anything, you literally just weren't asking about that. I thought your question was very specific.

No problem. I didn't assume it either. It's just that phrases not only say things literally, but also implictly, and the implication here *could* have been "he didn't ask for it, but he should have". I guess that's why it stood out. But no offense taken whatsoever. I had my answer plus extra information.

[DaveGamble]

Cool. Let's fix you up with some more info.
You have Orfanidis and MZTi, which are great starts, and neither of them are conformal maps, and they will both get you pretty acceptable results - in some cases, and with tweaking, even excellent results.

So, the root of your first question, if I've understood it right is: What property of a conformal map guarantees preservation of stability?

The shortest answer is this:
If, for all s=a+ib for a<=0, the mapping generates z=c+id where c*c+d*d<=1, then stability is guaranteed.

This formalism means, explicitly, that any point in the left-hand side of the s plane will map to a point inside the unit circle in the z-plane. This is because a stable system in s has poles only on the lefthand side (to make exponentials decay), and a stable system in z has all poles inside the unit circle.

So, BLT (I'm ignoring the 2/T, because it's boring): s=(z-1)/(z+1)
Let's wiggle that:
(z+1)s = (z-1)
zs + s -z +1 =0
z(s-1) + (s+1)=0
z=-(s+1)/(s-1)
z=(1+s)/(1-s)
We want to know about magnitude of Z, so, let's call s=a+ib:
|z|^2=((1+a)^2+b^2) / ((1-a)^2+b^2)
|z|^2=(1+2a+a^2 +b^2) / (1-2a+a^2+b^2)
|z|^2=(1-2a+a^2+b^2+4a) / (1-2a+a^2+b^2)
|z|^2=1 + 4a/(1-2a+a^2+b^2)
|z|^2=1 + 4a/((1-a)^2+b^2)

Now, we're interested in ensuring that |z|^2 <=1.
This means we need the second part (4a/((1-a)^2+b^2)) to be -ve.
By inspection, the denominator (1-a)^2+b^2 is strictly positive.
Hence, we need 4a to be -ve.

This is what we hoped to prove.

So that's a sketch of what such a proof entails. A little algebra, and the constraints we're looking for.
I suspect that Simpson's rule doesn't offer the same guarantee. But you could rearrange in the same way (reverse to get z as a function of s, look at the magnitude of z, see what must be true about s), and see whether you can get a similar constraint (actually, you already empirically know that you can't get as simple a constraint as for trapezoidal, or your filters would be stable!) and see what it looks like.

Hope this is helpful

Dave.

[mystran]

In the mean time, I've found more alternatives (MTZi), and I will definitely read up on Orfanidis. Thanks for that too.

Don't bother with Orfanidis. The Orfanidis approach leads to a crazy amount of computationally expensive and often ill-conditioned math (not to mention you need to handle every different type of filter separately), with results that are usually somewhat better than BLT... which in this case is really just a technical way of saying "tends to suck a bit less."

In contrast MZTi starts from an exact (but aliased) transform of the poles, with the expectation that the spectrum of the error (aliasing) is relatively smooth (which most of the time is entirely reasonable assumption) and so can be corrected with a low-order filter. As a result the whole thing reduces to a simple FIR design problem (and in fact you can use whatever FIR design procedure you like, though the one in the paper is kinda nice if you want to preserve order). Best part is that you can trivially improve it by simply using a better correction filter (or design one for a different error measure, or correct several filters together, or whatever).

I don't really see much value in using anything based on BLT for any LTI problem. It's more useful in non-LTI (especially non-linear) cases where dealing with transfer functions directly is "inconvenient" (so you just want a method that works "as-is" whatever the linearized transfer function happens to be) but if you're not constrained by such consideration then there's not much point in trying to deal with the warping. [Observe that in non-linear cases aliasing considerations usually force some oversampling, which mitigates the warping issue to some extent]