DSPRelated.com
Forums

Lost with direct-form filter equations (inverting)

Started by jtp_1960 March 7, 2007
jtp_1960 wrote:

> 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
> http://img261.imageshack.us/img261/1102/ > new88coeffsvsummingers4oi0.png
Are you expecting the curves to coincide? They're not supposed to; rather, if you normalize the blue curve for unit gain at what this plot calls 500 Hz then it should appear shifted down one octave from the normalized 44.1 kHz filter. And the red curve should hit 0 dB in this plot not at 1 kHz but at 44.1/48 = 919 Hz.
> As expected, the difference is not meaningful in practice and > this could easily be compensated w/ channel panning.
I don't understand. Do you mean stereo panning?
> 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.
What noises? Martin -- Research is what I'm doing when I don't know what I'm doing. --Wernher von Braun
In article <1173875134.989477@localhost>, martin.eisenberg@udo.edu 
says...
> > >jtp_1960 wrote: > >> hmm.. yes, the curve 'form' looks good in TobyBears >> FilterExplorer when use those coefficients: >> >> http://img76.imageshack.us/img76/2259/riaa88errornt7.png > >I don't think that's what you meant to show us, so I can't say much >there. But note that the Filter Explorer uses the reverse naming >convention for the coefficients from Frederick Umminger's submission. >So you need to put your b vector in its a fields and vice versa. > >>>I notice that your plot of the original curve only goes up to >>>about -25 dB while the true value is almost -20 dB. > >> Yes, I noticed that too ... maybe this is related to the >> C.W.Buddes "VST Plugin Analyzer" which one I'm using (it's >> working @ 44.1kHz all the time and there are no samplerate >> setting option available in software or in soundcard settings >> panel... and there are some other 'strange' issues w/ it as >> well)? > >It may be that the program's spectrum analyzer has some top-end >droop, I don't know. But it's most probably nothing to do with the >sampling rate because the same coefficients will give the same curve >at any rate, just with different numbers along the frequency axis.
I can't see the first part of this thread, but if you are trying to do an IIR simulation of the RIAA phono de-emphasis curve (assuming s-plane poles at 50.5 and 2122 Hz and an s-plane zero at 5005. Hz), here are some good minimum-phase magnitude approximations. (The RIAA de-emphasis is minimum-phase in the analog domain.) 44.1 kHz: SUPPLY # POLES IN Z-PLANE (<=10):2 Zero # Real Imag. 1 -0.2014898 0.000000 2 0.9233820 0.000000 Pole # Real Imag. 1 0.7083149 0.000000 2 0.9924091 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.2239207dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 30 degrees where the "phase error" is computed after a constant delay is added or subtracted to make the phase error equiripple SUPPLY # POLES IN Z-PLANE (<=10):3 Zero # Real Imag. 1 -0.5374877 0.000000 2 -0.9884768E-01 0.000000 3 0.9307718 0.000000 Pole # Real Imag. 1 -0.4600918 0.000000 2 0.7371564 0.000000 3 0.9928704 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0113530dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 23 degrees SUPPLY # POLES IN Z-PLANE (<=10):4 Zero # Real Imag. 1 -0.6929308 0.000000 2 -0.3386390 0.000000 3 -0.5961402E-01 0.000000 4 0.9311520 0.000000 Pole # Real Imag. 1 -0.6658348 0.000000 2 -0.2488928 0.000000 3 0.7389606 0.000000 4 0.9928295 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0005780dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 21 degrees 48 kHz: SUPPLY # POLES IN Z-PLANE (<=10):2 Zero # Real Imag. 1 -0.1766069 0.000000 2 0.9321590 0.000000 Pole # Real Imag. 1 0.7396325 0.000000 2 0.9931330 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.1395898dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 24 degrees SUPPLY # POLES IN Z-PLANE (<=10):3 Zero # Real Imag. 1 -0.4646165 0.000000 2 -0.8130194E-01 0.000000 3 0.9364602 0.000000 Pole # Real Imag. 1 -0.3741387 0.000000 2 0.7568389 0.000000 3 0.9934040 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0037544dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 16 degrees SUPPLY # POLES IN Z-PLANE (<=10):4 Zero # Real Imag. 1 -0.6147195 0.000000 2 -0.2762502 0.000000 3 -0.4728733E-01 0.000000 4 0.9365811 0.000000 Pole # Real Imag. 1 -0.5753610 0.000000 2 -0.1881046 0.000000 3 0.7574483 0.000000 4 0.9934112 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0000998dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 15 degrees 88.2 kHz SUPPLY # POLES IN Z-PLANE (<=10):2 Zero # Real Imag. 1 -0.1168735 0.000000 2 0.9648312 0.000000 Pole # Real Imag. 1 0.8590646 0.000000 2 0.9964002 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0081862dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 3 degrees SUPPLY # POLES IN Z-PLANE (<=10):3 Zero # Real Imag. 1 -0.3159579 0.000000 2 -0.4655857E-01 0.000000 3 0.9649734 0.000000 Pole # Real Imag. 1 -0.2060289 0.000000 2 0.8597030 0.000000 3 0.9964089 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0000096dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 2 degrees 96 kHz: SUPPLY # POLES IN Z-PLANE (<=10):2 Zero # Real Imag. 1 -0.1141486 0.000000 2 0.9676817 0.000000 Pole # Real Imag. 1 0.8699137 0.000000 2 0.9966946 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0057028dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 2.4 degrees SUPPLY # POLES IN Z-PLANE (<=10):3 Zero # Real Imag. 1 -0.3096394 0.000000 2 -0.4513594E-01 0.000000 3 0.9677730 0.000000 Pole # Real Imag. 1 -0.1992839 0.000000 2 0.8703280 0.000000 3 0.9967002 0.000000 MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0000046dB MAXIMUM PHASE ERROR FROM 0.00 Hz TO 20000.00 Hz IS ~+/- 1.6 degrees
>Are you expecting the curves to coincide? They're not supposed to; >rather, if you normalize the blue curve for unit gain at what this >plot calls 500 Hz then it should appear shifted down one octave from >the normalized 44.1 kHz filter. And the red curve should hit 0 dB in >this plot not at 1 kHz but at 44.1/48 = 919 Hz.
No. OK
> >> As expected, the difference is not meaningful in practice and >> this could easily be compensated w/ channel panning. > >I don't understand. Do you mean stereo panning?
Neither do I :) ... maybe I meaned overal gain (but was thinking through the test too much ... there it would be panning L/R.
>> 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. > >What noises?
Maybe it's just a phase effect ... least not distortion. Something like the dynamics becomes flatten and sound muffled (I can host link for an audio example later if needed). Those F. Ummingers coefficicents were equal for both, 44.1 and 48 kHz ... now when R. Orban kindly published his data (unigue values for 44.1kHz too), this issue is there not anymore in my implementation since I can now use unique cofficents for 44.1kHz data. jtp
> >I can't see the first part of this thread, but if you are trying to do >an IIR simulation of the RIAA phono de-emphasis curve (assuming s-plane >poles at 50.5 and 2122 Hz and an s-plane zero at 5005. Hz), here are >some good minimum-phase magnitude approximations. (The RIAA de-emphasis >is minimum-phase in the analog domain.) >
WOW!! Thank you very much. I got 'transferred' those 2P2Z Re/Im values to a/b coefficients and already tried those versions ... results looks excellent (48kHz is equal w/ the one I already have, 44.1 gives ~bout the same frequency response as 48kHz but then, 88.2 and 96 gives a bit different form FR curve in PluginAnalyzer ... I think it was explained already why in this way) even I thought the accuracy suffers because of the accuracy of your data (7 decimals). The 44.1kHz data playback issue is also OK now as I pointed earlier. Need to try those other filter forms as well when I just have some spare time needed for calculations/coding. Time to start learning/preparing the Subsonic filter part, I quess. Thanks to all of you for your kind helping/guidancing. jtp
Robert Orban wrote:

> I can't see the first part of this thread, but if you are trying > to do an IIR simulation of the RIAA phono de-emphasis curve > (assuming s-plane poles at 50.5 and 2122 Hz and an s-plane zero > at 5005. Hz), here are some good minimum-phase magnitude > approximations.
Neat, thanks! How did you make them? Martin -- Teach a man to make fire, and he will be warm for a day. Set a man on fire, and he will be warm for the rest of his life. --John A. Hrastar
jtp_1960 wrote:

> now when R. Orban kindly published his data (unigue values for > 44.1kHz too), this issue is there not anymore in my > implementation since I can now use unique cofficents for 44.1kHz > data.
So you're all set now, congrats! Martin -- Whatever you do will be insignificant, but it is very important that you do it. --Mahatma Gandhi
In article <1174071885.92833@localhost>, martin.eisenberg@udo.edu 
says...
> > >Robert Orban wrote: > >> I can't see the first part of this thread, but if you are trying >> to do an IIR simulation of the RIAA phono de-emphasis curve >> (assuming s-plane poles at 50.5 and 2122 Hz and an s-plane zero >> at 5005. Hz), here are some good minimum-phase magnitude >> approximations. > >Neat, thanks! How did you make them?
I used a program I wrote (in ye olde Fortran :-). The outline goes as follows: Given a desired magnitude response in the z-plane, there exists a response in a frequency-warped u-plane that, when bilinear-transformed to the z-plane, creates the desired z-plane magnitude response. -Compute the [magnitude response]^2 of the s-plane prototype on a grid. This is the square of the desired z-plane response. -Warp the frequency axis by using the bilinear transform, recognizing that we are approximating using omega^2 as our frequency variable. The warp maps Nyquist to infinity. -Make a least-squares rational approximation (i.e., ratio of polynomials) to the values on the frequency grid. (I used the Numerical Recipes routine RATLSQ, which uses Chebychev polynomials.) -Refine the approximation to make the fractional error minimax by using Remez's Second Algorithm [which applies to rational functions; it's not the same as the Remez algorithm used in the classical MPR FIR design program; see Forman S. Acton, Numerical Methods That Work (Revised Edition), Washington D.C., American Mathematical Society, 1990, pp 310- 314] -Transform the result into the z-plane in two steps. The first recognizes that we have been approximating using the magnitude square function, so we must take the square roots of the poles and zeros of the approximated rational function, taking the negative real parts to guarantee a stable, minimum-phase function. The second step is to apply the bilinear transform to the result of the first step. This yields the final z-plane poles and zeros. There are some "interesting" numerical issues in making this procedure work, mainly because the Remez update formulas require solving a system of mildly nonlinear equations that tend be ill-conditioned. The nice thing about the algorithm is that the frequency-warping moves Nyquist to infinity and thus increases the resolution of the approximation close to Nyquist, which is where difficulties often occur.
Robert Orban wrote:

> -Transform the result into the z-plane in two steps. The first > recognizes that we have been approximating using the magnitude > square function, so we must take the square roots of the poles > and zeros of the approximated rational function
Hmm, I don't get this. But I can't really think now either since I'm moving. Thanks for the explanation anyway. Martin -- Quidquid latine scriptum est, altum videtur.
In article <1175161226.968714@localhost>, martin.eisenberg@udo.edu 
says...
> > >Robert Orban wrote: > >> -Transform the result into the z-plane in two steps. The first >> recognizes that we have been approximating using the magnitude >> square function, so we must take the square roots of the poles >> and zeros of the approximated rational function > >Hmm, I don't get this. But I can't really think now either since I'm >moving. Thanks for the explanation anyway.
It's one of those things that's much clearer in a formal exposition with the appropriate equations. But USENET, being text based, is not very good for that sort of thing.
What you need is Spectral Transformation of filters[A. C.
Constantinides, Spectral transformations for digital filters, Proc.
IEE 117, 1585-1590 (Aug. 1970)]. So, The transformation is: z^^-1=-
z'^-1, 'cuz |Hlpf(z)|+|Hhpf(z)|=1. The solution is not so complicated
as yours:
b'=[b0 -b1 b2],
a'=[1 -a1 a2].