[futurebasic] Re: [FB] off-topic: algorithmic help

Message: < previous - next > : Reply : Subscribe : Cleanse
Home   : September 2010 : Group Archive : Group : All Groups

From: Robert Purves <listrp@...>
Date: Wed, 22 Sep 2010 09:47:16 +1200
Waverly wrote:

>> The implementation of IIRFilter() below seems entirely satisfactory, at least for a sampling rate of 44.1 kHz.
>> From the numerical standpoint it is vital to use double precision for the calculation and the recursively-used elements of y().
>> I seem to remember from a previous thread that Waverly's data type is not double, but short integer. Some care would be needed in adapting IIRFilter().

> Excellent demo!  I'm hoping to have time to adapt it to my existing code in the next day or two.

Please use this in your trials:

/*
7th-order IIR A-weighting filter.
http://dev.aubio.org/browser/src/temporal/a_weighting.c
http://dev.aubio.org/browser/src/temporal/a_weighting.h
In this simple implementation, inData() and outData() must not be the same array.
*/
local mode
local fn A_WeightingFilter( n as long, inData(_maxLong) as double, outData(_maxLong) as double )
'~'1
dim as long    j, k
_filterOrder = 7
dim as double  a(_filterOrder - 1), b(_filterOrder - 1) 
// coefficients for 44.1 kHz sample rate
b(0) =  2.557411252042575133813784304948057979345321655273437500e-01
b(1) = -5.114822504085150267627568609896115958690643310546875000e-01
b(2) = -2.557411252042575133813784304948057979345321655273437500e-01
b(3) =  1.022964500817030053525513721979223191738128662109375000e+00
b(4) = -2.557411252042575133813784304948057979345321655273437500e-01
b(5) = -5.114822504085150267627568609896115958690643310546875000e-01
b(6) =  2.557411252042575133813784304948057979345321655273437500e-01
a(0) =  1.000000000000000000000000000000000000000000000000000000e+00
a(1) = -4.019576181115832369528106937650591135025024414062500000e+00
a(2) =  6.189406442920693862674852425698190927505493164062500000e+00
a(3) = -4.453198903544116404873420833609998226165771484375000000e+00
a(4) =  1.420842949621876627475103305187076330184936523437500000e+00
a(5) = -1.418254738303044160119270600262098014354705810546875000e-01
a(6) =  4.351177233495117681327801761881346465088427066802978516e-03

for j = 0 to n - 1
outData(j) = b(0)*inData(j)
for k = 0 to _filterOrder - 1
long if ( j >= k )
outData(j) += b(k)*inData(j - k) - a(k)*outData(j - k)
end if
next
next 
end fn



> I was able to locate a functioning filter and started putting into a test rig yesterday.

>   http://users.skynet.be/solaris/linuxaudio/downloads/aweight.tar.bz2

There doesn't seem to be any documentation of the filter design. From inspection of the code, it consists of multiple 2nd-order sections. That is a commonly-recommended way to design IIR filters, but completely different from A_WeightingFilter() above.

Robert P.