FPGA Central

First of all, thank you a lot for the information about the coefficient
for a A weighting filter.
I tried to implement it with delphi, but the whole filter seems to b
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 weightin
filter.

Tomy

On Tue, 13 Oct 2009 05:09:05 -0500, "TomyN" <[email protected]> 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

Hi Greg,

the samplerate is 48kHz.


Tomy

On Tue, 13 Oct 2009 07:07:55 -0500, "TomyN" <[email protected]> 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

Hi Greg,

thanks for pointing out my mistakes. As you might have guessed, I'm no
familiar with digital filters at all :-(

Am I right if I assume that x(n) are the input samples, where n is th
most current one, and that y(n) is the latched result, where y(n-1) is th
result of the last calculation?

Thanks again

Tomy

On Tue, 13 Oct 2009 12:53:49 -0500, "TomyN" <[email protected]> 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

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 becaus
the software I'm trying to use asumes a -oo dB Value at the Nyquis
frequency, which makes the slope steeper towards half the samplin
frequency.
Is there any 'basic reading' which I should/could do?

Regards

Tomy

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.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Jerry Avins <[email protected]> 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

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 <[email protected] 7>,
[email protected] says...
>
>
>Jerry Avins <[email protected]> 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

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 <[email protected]> wrote in
news:[email protected]:

> 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 <[email protected] 7>,
> [email protected] says...
>>
>>
>>Jerry Avins <[email protected]> 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
>
>

In article <[email protected] 0>,
[email protected] says...
>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?

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.