>Plot shows that the digital filter is bent a bit.Plot - https://drive.google.com/file/d/0B2orKJuZJ5YOcE9VZll6X1ZiRzA/view?usp=sharing --------------------------------------- Posted through http://www.DSPRelated.com
>Plot shows that the digital filter is bent a bit.Plot - https://drive.google.com/file/d/0B2orKJuZJ5YOcE9VZll6X1ZiRzA/view?usp=sharing --------------------------------------- Posted through http://www.DSPRelated.com
>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.Sorry bringing this old thread alive. ATM, I'm working with non-RIAA de-emphasis and trying to find solution to get accurate playback curve especially for low sample rates (44.1/48 kHz). I'm using MZT transformation. While tweaking the optimization method I'm using, I did some comparisons against coefficients you calculated using this method. By data given few posts before this, case 2nd order filter: 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 Now, I'm wondering if this information is valid in reality? When compared against analog model, error is more. Here's my magnitude comparison: fs = 44100Hz p1 = 0.9924091 p2 = 0.7083149 z1 = 0.9233820 z2 = -0.2014898 Normalized (0dB @ 1kHz) b = 1.0 -1,700724 0,70293815 a = 0,23806122 -0,17185454 -0,04429178 Magnitude comparison: Fc Analog Digital Error 20 19,274148 19,347428 0,0732799702691 30 18,592575 18,734328 0,1417527973581 40 17,791989 18,002063 0,210074022786 60 16,100613 16,419205 0,3185918991332 120 11,847648 12,28887 0,441222223504 250 6,676967 7,065834 0,3888670946689 500 2,647603 2,812765 0,1651626143418 1000 0 0 0 2000 -2,588541 -2,436851 0,1516895410375 4000 -6,605229 -6,194421 0,410808126892 6000 -9,598882 -9,156793 0,4420892984055 8000 -11,894116 -11,521667 0,3724489269895 10000 -13,734342 -13,477005 0,257337832836 12000 -15,263689 -15,132886 0,1308032495094 15000 -17,156907 -17,153695 0,0032121877763 16000 -17,707729 -17,704712 0,0030176528852 17000 -18,226207 -18,188482 0,0377247338064 18000 -18,715877 -18,599022 0,1168555511977 19000 -19,179733 -18,929519 0,250213973018 20000 -19,620332 -19,173155 0,4471763603379 As this comparison is done only at those 20 points you probably don't see the max error peak position there. IMO, the real maximum error compared to analog is close 0.45dB. Am I wrong? Plot shows that the digital filter is bent a bit. Why I'm asking this is because of I did prepare same filter without scientific optimization methods: fs = 44100Hz p1 = 0,99289463 p2 = 0,72462291 z1 = 0,93117565 z2 = -0,3638376 Normalized (0dB @ 1kHz) b = 1.0 -1,72474828 0,72665356 a = 0,23063639 -0,17280589 -0,03906943 Magnitude comparison: Fc Analog Digital Error 20 19,274148 19,22913 -0,0450181850077 30 18,592575 18,547585 -0,044990382342 40 17,791989 17,747038 -0,044951469567 60 16,100613 16,055773 -0,0448403594181 120 11,847648 11,803406 -0,0442421344389 250 6,676967 6,635353 -0,0416136297316 500 2,647603 2,615705 -0,0318976557532 1000 0 0 0 2000 -2,588541 -2,521626 0,0669152280886 4000 -6,605229 -6,500148 0,1050812334451 6000 -9,598882 -9,540168 0,0587146975408 8000 -11,894116 -11,921155 -0,0270393790294 10000 -13,734342 -13,863248 -0,1289056471034 12000 -15,263689 -15,487502 -0,2238126905102 15000 -17,156907 -17,437916 -0,2810087918641 16000 -17,707729 -17,96218 -0,2544506635893 17000 -18,226207 -18,419286 -0,1930787911619 18000 -18,715877 -18,804655 -0,0887782220308 19000 -19,179733 -19,113065 0,0666683158919 20000 -19,620332 -19,339311 0,281020405759 In my case the error is -0.281...dB to +0.281...dB which (IMO) means that the maximum error against analog mode is 0.281...dB. Am I right. Juha --------------------------------------- Posted through http://www.DSPRelated.com
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. (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 > >I am also playing with digital implementation of RIAA de-emphasis filters, and have stumbled on this thread (a little late, from the look of it). I wonder are Robert Orban's optimizations for fs=176k4Hz and fs=192kHz available anywhere? Thanks
On Apr 4, 9:47 am, Robert Orban <donotre...@spamblock.com> wrote:> In article <1175613530.446199.291...@o5g2000hsb.googlegroups.com>, > jerry_cq...@yahoo.com says... > > > > >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]. > > Does Constantinides' algorithm create minimax magnitude error curves > over a user-specified frequency range?I don't think so. H1(z)=H0(-z), this is simply delay the magnitude response curve 180 degree. Since the given filter is a hpf, the resulted filter is a lpf.
In article <1175613530.446199.291610@o5g2000hsb.googlegroups.com>, jerry_cq_cn@yahoo.com says...> > >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].Does Constantinides' algorithm create minimax magnitude error curves over a user-specified frequency range?
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].
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.
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 functionHmm, 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 <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.
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