DSPRelated.com
Forums

Lost with direct-form filter equations (inverting)

Started by jtp_1960 March 7, 2007
jtp_1960  wrote:
...
> So, as those 88.2/96 filters freezes every software I've tried this far,
You definitely tried the wrong software! Those are probably the ones that use Erik's method (http://groups.google.com/group/comp.dsp/ browse_frm/thread/92f44f954dfc3645/0d3721014782cd54?#0d3721014782cd54) to compute the frequency response :-).
> it must be either those coefficient values or hardware in question that is > giving those issues (since both, 44.1/48kHz filters works well).
Are you talking about the coefficients that you mentioned at the beginning, ie. b = [ 0.04872204977233 -0.09076930609195 0.04202280710877] a = [ 1.00000000000000 -0.85197860443215 -0.10921171201431] ? You cannot invert this filter by simply swapping "b" and "a" because "b" has a zero outside the unit circle. This means that the filter does not have a causal inverse (contrary to common belief it has an acausal and stable inverse that can be approximated using an FIR, for example). Using Greg Berchin's useful FDLS IIR design method [1], I came up with the following (barely) stable inverse filter: B = [ 12.0774464957845 -11.8251514992744 0.339824392108838] A = [ 1.0000000000000 -1.25660022712513 0.257312844988975] which is close but not quite what you want (max error is 15dB at 1 kHz). Tweaking the weights and playing with the orders could get you closer, but I don't have the time.
>If this is the culprit for these issues I'm having here, are there any >'mathematical' methods to get this fixed w/o loosing the accuracy Umminger >got into those original coefficients or, are those 88.2/96kHz coefficients >just unusable when inverted (least the one works well :) )?
Yes, that is exactly the problem. Given a second order section denominator a = [1 a1 a2], then the filter is stable if the coefficients lie inside the "stability triangle". If you think of a coordinate system with a1 as the x-axis and a2 as the y-axis, the stability triangle is bounded by the points (-2,-1), (0,1) and (2,-1). The coefficients may also lie on the triangle boundary except at points (-2,-1) and (2,-1). Regards, Andor [1] Berchin, G: "Precise Filter Design", IEEE Signal Processing Magazine, January 2007.
Oops. I should first learn how to use Greg's method before posting.
After using it correctly, it came up with the exact but unfortunately
unstable inverse ...

Another option would be to use one of the methods for IIR design
posted in another thread ("Generating 1/f^2 and 1/f^3 noise"). Or to
reflect the unstable pole about the unit circle. This should give
exact magnitude but result in a phase error.

Regards,
Andor

Andor wrote:
> Oops. I should first learn how to use Greg's method before posting. > After using it correctly, it came up with the exact but unfortunately > unstable inverse ... > > Another option would be to use one of the methods for IIR design > posted in another thread ("Generating 1/f^2 and 1/f^3 noise"). Or to > reflect the unstable pole about the unit circle. This should give > exact magnitude but result in a phase error.
Does phase error have any significance for noise? Jerry -- Engineering is the art of making what you want from things you can get. ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Andor wrote:

> jtp_1960 wrote: > ... >> So, as those 88.2/96 filters freezes every software I've tried >> this far, > > You definitely tried the wrong software! Those are probably the > ones that use Erik's method (http://groups.google.com/ > group/comp.dsp/browse_frm/thread/92f44f954dfc3645/ > 0d3721014782cd54?#0d3721014782cd54) to compute the frequency > response :-).
Is there more to that statement than the joke?
> b = [ 0.04872204977233 -0.09076930609195 0.04202280710877] > a = [ 1.00000000000000 -0.85197860443215 -0.10921171201431]
> You cannot invert this filter by simply swapping "b" and "a" > because "b" has a zero outside the unit circle.
Oops, I didn't even consider that Frederick's filters might not be minimum-phase. But since the original designs don't attempt to match the phase response of an analog RIAA filter there's no need to reapproximate. Jtp -- this means that you can fix the problem as follows: factor the numerator polynomial b(1)*z^2+b(2)*z+b(3), replace any roots of modulus greater than unity by their reciprocals (make sure that any conjugate pairs stay conjugate), expand the polynomial again and grab the new coefficients. In the present case the result is b2 = [.04872204977233 -.090422059905446592 .041724372373403292] The filter (b2,a) has the same magnitude response as (b,a) but with a stable inverse. Proceed as before to get the inverse. Again, reflecting the large zeros into the unit circle is justified here only because the phase response was unconstrained to begin with. Whether compensating the magnitude alone will give you good sound from an LP is of course another matter that also pertains to Frederick's designs. If you find it doesn't, perhaps try the method Andor cited on the analog inverse RIAA spec complex response directly. Martin -- What is it: is man only a blunder of God, or God only a blunder of man? --Friedrich Nietzsche
On 15 Mrz., 15:51, Jerry Avins <j...@ieee.org> wrote:
> Andor wrote: > > Oops. I should first learn how to use Greg's method before posting. > > After using it correctly, it came up with the exact but unfortunately > > unstable inverse ... > > > Another option would be to use one of the methods for IIR design > > posted in another thread ("Generating 1/f^2 and 1/f^3 noise"). Or to > > reflect the unstable pole about the unit circle. This should give > > exact magnitude but result in a phase error. > > Does phase error have any significance for noise?
Why noise? The filters are supposed to be used for de-emphasis of an analogue emphasis system. In that case, phase matters.
On 15 Mrz., 16:26, Martin Eisenberg <martin.eisenb...@udo.edu> wrote:
> Andor wrote: > > jtp_1960 wrote: > > ... > >> So, as those 88.2/96 filters freezes every software I've tried > >> this far, > > > You definitely tried the wrong software! Those are probably the > > ones that use Erik's method (http://groups.google.com/ > > group/comp.dsp/browse_frm/thread/92f44f954dfc3645/ > > 0d3721014782cd54?#0d3721014782cd54) to compute the frequency > > response :-). > > Is there more to that statement than the joke?
You mean, do I actually know whether the mentioned software(s) actually use the FFT of the impulse response to compute and display the frequency response? No, of course I don't know. I'm guessing.
Andor wrote:
> On 15 Mrz., 15:51, Jerry Avins <j...@ieee.org> wrote: >> Andor wrote: >>> Oops. I should first learn how to use Greg's method before posting. >>> After using it correctly, it came up with the exact but unfortunately >>> unstable inverse ... >>> Another option would be to use one of the methods for IIR design >>> posted in another thread ("Generating 1/f^2 and 1/f^3 noise"). Or to >>> reflect the unstable pole about the unit circle. This should give >>> exact magnitude but result in a phase error. >> Does phase error have any significance for noise? > > Why noise? The filters are supposed to be used for de-emphasis of an > analogue emphasis system. In that case, phase matters.
Deleted within seconds of posting. I had been thinking of another thread. Jerry -- Engineering is the art of making what you want from things you can get. &macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;&macr;
On Mar 15, 12:13 pm, "Andor" <andor.bari...@gmail.com> wrote:
> On 15 Mrz., 15:51, Jerry Avins <j...@ieee.org> wrote: > > > Andor wrote: > > > Oops. I should first learn how to use Greg's method before posting. > > > After using it correctly, it came up with the exact but unfortunately > > > unstable inverse ... > > > > Another option would be to use one of the methods for IIR design > > > posted in another thread ("Generating 1/f^2 and 1/f^3 noise"). Or to > > > reflect the unstable pole about the unit circle. This should give > > > exact magnitude but result in a phase error. > > > Does phase error have any significance for noise? > > Why noise? The filters are supposed to be used for de-emphasis of an > analogue emphasis system. In that case, phase matters.
yeah, but if the pre-emphasis was not minimum-phase (i.e. zeros outside the unit circle which would become nasty-ass poles outside the unit circle in the de-emph filter), this pre-emphasis has no realizable inverse (meaning, in my book, that it was not well designed). now if phase is that important (i've argued both sides of this "is phase audible" issue, BTW), then, if you're willing to put up with some throughput delay, a decent approximation to the inverse filter can be done with an FIR using any of these FIR design methods (and that trick in the dspguru tricks page about using PM to design for magnitude *and* phase). otherwise (if the phase relationship between frequencies is *not* all that critical), just reflect the nasty poles into the interior of the unit circle. dunno if i can add anything else. r b-j
Andor wrote:

> On 15 Mrz., 16:26, Martin Eisenberg <martin.eisenb...@udo.edu> > wrote: >> Andor wrote: >> > jtp_1960 wrote: >> > ... >> >> So, as those 88.2/96 filters freezes every software I've >> >> tried this far, >> >> > You definitely tried the wrong software! Those are probably >> > the ones that use Erik's method (http://groups.google.com/ >> > group/comp.dsp/browse_frm/thread/92f44f954dfc3645/ >> > 0d3721014782cd54?#0d3721014782cd54) to compute the frequency >> > response :-). >> >> Is there more to that statement than the joke? > > You mean, do I actually know whether the mentioned software(s) > actually use the FFT of the impulse response to compute and > display the frequency response? No, of course I don't know. I'm > guessing.
Of course you are -- let me clarify. While I mainly asked for jtp's benefit, I also wondered if you might have commonly had graphing libraries crash on you when fed with humongous numbers, infinities, or NaN's as might result from transforming a divergent response? Martin -- Quidquid latine scriptum est, altum videtur.
Thanks a lot for everyone for your guidance and suggestions.

Yes, those modified coefficients turned the filter stable.
I have the well working 48kHz filter here as a reference (red curve)
against the 'modified' 88.2kHz filter (blue curve). Curve form is not
exact 'enough' yet (though, you can't hear the difference if the volume is
set to equal level and w/ standard EQ plugins I couldn't get even close to
that result).

http://img261.imageshack.us/img261/1102/new88coeffsvsummingers4oi0.png

(NOTE: value 2.929 showing up there is just for compensation to get the
1kHz dropped to ~0dB level ...).


As expected, the difference is not meaningful in practice and this could
easily be compensated w/ channel panning.

http://img254.imageshack.us/img254/6946/filterinvsthostda5.png

(NOTE: curve colouring opposite in SPAN. Playback source ia a wav file in
this case (though ripped earlier from LP)).

Results from 48kHz filter are excellent already but, I would like to get
it as good for all standard samplerates so I could automate the plugin a
bit (filter selected by the samplerate as for an example). (Maybe it's the
method but,) It looks like when lower samplerate settings/data is used w/ a
filter intend for higher samplerate, lets say, 44.1kHz audio through
88.2kHz filter, there are some strange noises added then.

Perhaps I need to look that 'reflecting' method (need to take a deeper
look over it 1st ...) first but, I think, I'll try to get the compensation
curve fit better by preparing a VST plugin w/ sliders for to adjust those
a's manually as close to the 48kHz reference curve as possible and just
then grab those resulting coefficients from there and use them for 88.2
and 96 (maybe 176.4 and 192 too) kHz :).


jtp

Some additional mesurement data for this modified 88.2kHz filter:
Phase : http://img474.imageshack.us/img474/221/phase882wv6.png
Harmonic distortion :
http://img108.imageshack.us/img108/5606/hdistortion882ox9.png