First of all, thank you a lot for the information about the coefficients
for a A weighting filter.
I tried to implement it with delphi, but the whole filter seems to be
shifted towards higher frequencies.
My code is:
Function TdbFilter.DoAFilter(Const pIn: Pointer; const pOut: Pointer;
const Points: integer): boolean;
const
B0 = 0.01147155239724;
b1 = -0.0248548824653166;
b2 = 0.0323052801494283;
b3 = -0.0555571372813964;
b4 = 0.233784933167266;
b5 = 0.392882400411297;
b6 = -0.633249911806281;
b7 = -0.479494344932334;
b8 = 0.322061493940389;
b9 = 0.197659563091056;
b10 = 0.00299879461389451;
a0 = 1;
a1 = -0.925699454182466;
a2 = -0.992471193943543;
a3 = 0.837650096562845;
a4 = 0.22307303912603;
a5 = -0.158404327757755;
a6 = 0.0184103295763937;
var i: integer;
pi: pDouble;
po: pDouble;
oben: double;
unten: double;
j: integer;
begin
pI:= pIn;
po:= pOut;
for j:= 1 to points do begin
for i:= 9 downto 1 do nal[i]:= Nal[i-1]; //Schieben
nal[0]:= pI^;
oben:= b0 + b1*nal[0] + b2*nal[1] + b3*nal[2] + b4*nal[3] +
b5*nal[4] + b6*nal[5] + b7*nal[6] + b8*nal[7] +
b9*nal[8] + b10*nal[9];
unten:= a0 + a1*nal[0] + a2*nal[1] + a3*nal[2] + a4*nal[3] +
a5*nal[4] + a6*nal[5];
if unten <> 0 then po^:= oben/unten else po^:= 0;
inc(pi);
inc(po);
end;
end;
Any hints? Furthermore I'm looking for the values for a C type weighting
filter.
Tomy
Problem with A-weighting filter
Started by ●October 13, 2009
Reply by ●October 13, 20092009-10-13
On Tue, 13 Oct 2009 05:09:05 -0500, "TomyN" <Tomy@Tomy-pa.de> wrote:>First of all, thank you a lot for the information about the coefficients >for a A weighting filter. >I tried to implement it with delphi, but the whole filter seems to be >shifted towards higher frequencies.The filter that you used was designed for 48 kHz. What is your sampling rate?>Any hints? Furthermore I'm looking for the values for a C type weighting >filter.The C-weighting filter spec and frequency response are available in the literature. You can use any of a number of design techniques to approximate it, including FDLS. Greg
Reply by ●October 13, 20092009-10-13
Reply by ●October 13, 20092009-10-13
On Tue, 13 Oct 2009 07:07:55 -0500, "TomyN" <Tomy@Tomy-pa.de> wrote:>the samplerate is 48kHz.Okay, then let's look at your implementation. I don't speak delphi, but from context it appears that you are loading the time domain values into the frequency domain transfer function. If your nal[n] equals my input value x(n), then you have implemented: b0 + b1x(n) + b2x(n-1) + ... + b10x(n-9) ---------------------------------------- a0 + a1x(n) + a2x(n-1) + ... + a6x(n-6) That's just plain wrong in so many ways. With the coefficients provided, the time domain difference equation that implements the filter is: y(n) = b0x(n) + b1x(n-1) + ... + b10x(n-10) - a1y(n-1) - ... - a6y(n-6) Or the frequency response can be computed directly from the transfer function by substituting z=exp(sT) and setting s=jw: b0 + b1exp(-jwT) + b2exp(-j2wT) + ... + b10exp(-j10wT) ------------------------------------------------------ a0 + a1exp(-jwT)) + a2exp(-j2wT) + ... + a6exp(-j6wT) Greg
Reply by ●October 13, 20092009-10-13
Hi Greg, thanks for pointing out my mistakes. As you might have guessed, I'm not familiar with digital filters at all :-( Am I right if I assume that x(n) are the input samples, where n is the most current one, and that y(n) is the latched result, where y(n-1) is the result of the last calculation? Thanks again Tomy
Reply by ●October 13, 20092009-10-13
On Tue, 13 Oct 2009 12:53:49 -0500, "TomyN" <Tomy@Tomy-pa.de> wrote:>Am I right if I assume that x(n) are the input samples, where n is the >most current one, and that y(n) is the latched result, where y(n-1) is the >result of the last calculation?Correct. Greg
Reply by ●October 16, 20092009-10-16
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? Regards Tomy
Reply by ●October 16, 20092009-10-16
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 -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●October 16, 20092009-10-16
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? > > JerryBLT 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
Reply by ●December 17, 20092009-12-17
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
