[music-dsp] Re: Revisiting C++ Filtering Classes
thevinn at yahoo.com
Tue Jun 2 01:22:55 EDT 2009
> Date: Sat, 23 May 2009 19:57:26 +0200
> From: "Martin Eisenberg" <martin.eisenberg at udo.edu>
> So it seems your best options are:
> a) leave the tests in Prepare() alone and make it clear to
> subclass authors that the index passed to GetPole/Zero()
> over both members of any complex root pairs.
> b) have Prepare() count only up to "stages" and let
> specify each stage as a unit (two Complex values which
> combines with their conjugates if nonreal), rather than
> poles and zeros separately. Subclasses get to decide
> whether real
> roots are combined into 2nd-order sections or are turned
> stages of their own. I think this raises no technical
> once the cleverness in Set*Stage is dropped, and it saves
> currently wasted half of Prepare()'s computation.
I've been studying the Orfanidis book, also Digital Filters... by Andreas Antoniou, a couple of research papers by Orfanidis, and the DSP Filters cookbook series by John Lane et. al.
I have a much better grasp of the transfer function, the s/z plane, and what is going on inside these filters. I have identified the formula in the Orfanidis book that corresponds to the Butterworth poles in the code, as well as the highpass/bandpass/bandstop transformations.
I finally invented my own filter by taking H(s)=A(s)/B(s) bandpass and transforming it into a peaking filter by setting H'(s)=1+A(s)/B(s) and solving the resulting equation, getting it back into a rational polynomial.
At first, I thought that the computation of the conjugate poles was a waste. It was embarassing to compute them, only to have SetAStage() and SetBStage() throw them away. The current interface also suffers from the problem that constant intermediate calculations are re-computed for every pole and/or zero when they don't have to be. A better model would be a functional programming style, implemented by removing the state variables and switching over to having the derived class drive the setup of stages rather than going through the base.
It turns out that the calculation of the conjugate poles are in fact necessary. Witness:
for( int i=0;i<poles;i++ )
c=GetPole( i );
if( ::std::abs(c.imag())<1e-6 )
c=Complex( c.real(), 0 );
if( c.imag()==0 )
SetAStage( c.real(), 0 );
else if( c.imag()>0 )
SetAStage( 2*c.real(), -std::norm(c) );
In the case where the imaginary part of the pole or zero is sufficiently close to the i axis, the coordinate and its conjugate are merged together and get factored into the stages twice. SetAStage() is called twice with the same value c.real(), with the imaginary part set to 0, for the case where conjugate poles or zeroes are sufficiently close to the i axis (in this case within 1e-6 units).
Is my analysis correct? Is there an easy way to ignore the conjugates?
More information about the music-dsp