This is very cool Robert, I have used a magnitude squared method which is not nearly as good as your minimax solution. Do you have C or Matlab code that you could share for the minimax calculations? Al Clark www.danvillesignal.com Robert Orban <spambucket2413@earthlink.net> wrote in news:20091217-040135.523.0@Robert-Orban.news.giganews.com:> Here is a solution that uses the technique suggested by Al Clark of > splitting the A-weighting filter into lowpass and highpass sections, > which is a very useful technique for approximating filters like this. (I > got the s-plane pole and zero locations from here: > http://www.cross-spectrum.com/audio/weighting.html) > > The highpass section can be approximated using the matched-z transform > to an accuracy of about +/-0.001 dB, 0-20 kHz. (Ignore the 1E-29; it's a > cheap trick that my software uses to prevent a divide by 0.) > > s-plane zeros: > Zero # Real Imag. > 1 0.1000000E-29 0.000000 > 2 0.1000000E-29 0.000000 > 3 0.1000000E-29 0.000000 > 4 0.1000000E-29 0.000000 > s-plane poles: > Pole # Real Imag. > 1 -20.60000 0.000000 > 2 -20.60000 0.000000 > 3 -107.7000 0.000000 > 4 -737.9000 0.000000 > > Z-PLANE POLES/ZEROS @ 48 kHZ SR USING MATCHED-Z TRANSFORM > > Zero # Real Imag. > 1 1.000000 0.000000 > 2 1.000000 0.000000 > 3 1.000000 0.000000 > 4 1.000000 0.000000 > Pole # Real Imag. > 1 0.9973071 0.000000 > 2 0.9973071 0.000000 > 3 0.9860010 0.000000 > 4 0.9079274 0.000000 > > The lowpass section can be approximated by an almost-minimax > optimization as follows: > > s-plane poles: > Pole # Real Imag. > 1 -12200.00 0.000000 > 2 -12200.00 0.000000 > > SUPPLY DIGITAL SAMPLING FREQUENCY (kHz): 48 > > SUPPLY BEGINNING & END OF APPROXIMATION BAND (Hz): 0 20000 > > Zero # Real Imag. > 1 -0.6064625 0.000000 > 2 -0.2606820 0.000000 > 3 0.6582298E-01 0.000000 > Pole # Real Imag. > 1 -0.5328594 0.000000 > 2 0.2181849 0.4849923E-01 > 3 0.2181849 -0.4849923E-01 > > MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0002765dB > > Note that a true, stable 3rd order minimax approximation is unrealizable > because it requires a pole pair on the unit circle in the z-plane. > However, this approximation is very close to minimax. > > The 4th order lowpass approximation is fully minimax, but it's overkill: > > Zero # Real Imag. > 1 -0.7005693 0.000000 > 2 -0.4281911 0.000000 > 3 -0.1869575 0.000000 > 4 0.1364491E-01 0.000000 > Pole # Real Imag. > 1 -0.6638332 0.000000 > 2 -0.3389376 0.000000 > 3 0.2041116 0.1438447E-01 > 4 0.2041116 -0.1438447E-01 > MAXIMUM ERROR FROM 0.00 Hz TO 20000.00 Hz IS 0.0000056dB > > When the 3rd order lowpass and highpass functions are cascaded, they > approximate the A weighting curve +/- 0.0013 dB or so. > > > > > > > > In article <Xns9CA65A2CCCFA1aclarkdanvillesignal@69.16.185.247>, > aclark@danvillesignal.com says... >> >> >>Jerry Avins <jya@ieee.org> wrote in news:IYXBm.164314$nL7.111940 >>@newsfe18.iad: >> >>> TomyN wrote: >>>> Thanks again, >>>> >>>> I got it working. I tried to build the filter for C-Weighting >>>> >>>> s^2 / ((s + 129.4)^2 * (s + 76655)^2) >>>> >>>> by myself, but failed on both ends of the spectrum. I think thats >>>> because the software I'm trying to use asumes a -oo dB Value at the >>>> Nyquist frequency, which makes the slope steeper towards half the >>>> sampling frequency. >>>> Is there any 'basic reading' which I should/could do? >>> >>> Did you account for frequency warping? >>> >>> Jerry >> >>BLT will not work for the A or C filter. The double poles at 12k are >>much too close to the nyquist frequency which as pointed out will insert >>zeros at Pi. >> >>You can split the function in half and use BLT for the low frequency >>poles and zeros. You do need good bit precision in the implementation. >> >>This leaves you with a fitting problem for the 12k poles which is >>easier. >> >>There are also other approaches like Greg Berchin's FDLS method. This >>method is explained in Streamlining Digital Signal Processing - Lyons >> >>Al Clark >>www.danvillesignal.com > >
Problem with A-weighting filter
Started by ●October 13, 2009
Reply by ●December 17, 20092009-12-17
Reply by ●December 17, 20092009-12-17
In article <Xns9CE46D122E1B2aclarkdanvillesignal@69.16.185.250>, aclark@danvillesignal.com says...>This is very cool Robert, > >I have used a magnitude squared method which is not nearly as goodas your>minimax solution. > >Do you have C or Matlab code that you could share for the minimax >calculations?I wrote the program in Fortran 90. Unfortunately, the rights belong to my employer so I can't share it. Some months back I posted in comp.dsp a brief explanation of how the algorithm works in the context of a thread on accurate magnitude approximations to the RIAA playback curve. The even briefer explanation :-) is that it is a pure magnitude approximation that uses a Remez rational function optimization on a warped frequency axis [which gets rid a lot of the trig terms in the frequency repsonse formulas], where the starting point for the Remez optimization is a least-squares optimization using singular value decomposition. But the actual implementation took quite a bit of time and work; there are several thousand lines of Fortan involved.
