[music-dsp] Re: A Collection of C++ Classes for Digital Signal Processing

[Vinnie]

> From: "Martin Eisenberg"

Let me just open up by saying that I am by no means an expert on signal processing. I just started into the field from scratch a couple of weeks ago. Most of the algorithms in the code I published came from elsewhere. Paul Falstad's Java applet in particular (see the credits).

My "value-add" was in putting it all together in C++, making it easy to integrate (single header file), use of templates for flexibility in the underlying type and number of channels, and some optimizations (including SSE3), and testing in a professional (unreleased) application.

I was frustrated at what I perceived to be a lack of production quality C++ filtering code published under a commercial-use friendly license so I threw this together.

> | template<typename To>
> | ComplexT operator+ ( const ComplexT<To> &c )
> const;
> Having these operators makes the non-template overloads
> superfluous since their To == T instantiations will do
> exactly the same thing.

Yep I see that. I added the other operators as an afterthought and never cleaned out the old ones. However, I did provide a compiler switch for doing away with my Complex class completely. When I put in the trigonometric approximations I will see if I can get rid of the ComplexT class. I haven't figured out how the API for the fast setup of PoleFilter parameters is going to work yet. The Biquads were easy.

> recip(): Both versions miss a minus sign in the result
> imaginary
> part.

Are you saying it should read
	return std::complex<T>( n*c.real(), - n*c.imag() );

Apparently it doesn't matter because I tried it that way and it produces the same results. I tested the Chebyshev 2 filters (those are the only ones that use recip()).


> BrentMinimize():
> | static const CalcT SQRT_DBL_EPSILON =
> |   ::sqrt(DBL_EPSILON);
> 
> Being unfamiliar with Brent's method I'm not sure
> how important it is, but I suppose if you want to allow CalcT == float
> then the above line ought to use
> std::numeric_limits<CalcT>::epsilon().

Brent's method is a function minimization algorithm that is similar to gradient descent, but using parabolas. I use it to find the top of the ripple on chebyshev filters to normalize them for unity gain in the passband. I already do HintPassband() for as many cases as possible but for the bandpass chebyshev type I there is no "safe" place to sample the magnitude response so I must resort to a numerical technique.

> Also, a quick comparison with the treatment in [1] seems to
> indicate that you have switched the roles of that constant
> and parameter epsilon in the line below:
> 
> | tol = SQRT_DBL_EPSILON*::fabs(x) + epsilon;

I have to admit I am fuzzy on the math. I just took Brent's method from the Internet. I have Numerical Recipes but the code is restricted so I am prevented from using it. I don't understand the math behind Brent's method.

> That's not only a needless hack but also makes the
> runtime of Setup() overriders quadratic in the number of stages. Just
> pass an index into SetA/BStage.

Quadratic performance is unacceptible, especially after my efforts at implementing SetupFast()! Major oversight on my part.

I didn't understand exactly what Set#Stage() were doing but now that you've explained it, I understand better. I will fix it!

> I know that interferes with your trick of folding
> successive 1st-order sections (a bug: s->b[0]=x0 should use *= --
> but it's moot, read on)

Its not "my" trick and I didn't know what it was doing until you explained it. So there's a bug?

> which is kinda neat as long as coefficients
> are fixed. But when they vary in time such that section folding
> is incidentally triggered, all downstream history states will
> suddenly get operated upon by different coef sets.

I don't understand what you are saying. If coefficients vary then all stages need to be recalculated.

> The only way to ensure smooth state evolution is to forgo
> the optimization of section folding in the time-varying case.
> But you could have a separate function FoldStages() that also calls
> Clear() and invoking which constitutes a promise that stage
> orders won't change before the next call to Clear() or
> FoldStages().

What is state evolution? What is the time-varying case?

> xxBandPass::PassbandHint(): That should be a geometric, not
> arithmetic, mean; and except at bandwidths up to about pi/6
> you need to form the mean of BLT-unwarped values and warp that
> back. So the body would be:
>   return 2*atan(sqrt(tan(m_wc*0.5) * tan(m_wc2*0.5)));

I'm guessing BLT doesn't mean Bacon Lettuce and Tomato. For PassbandHint() all I wanted to do was choose a value for the angular frequency that was as far away as possible from anything that deviates from unity gain. I figured the middle of the passband was ideal since it is the farthest from the rolloffs at either side.

Are you saying that the formula you provided can be used in place of Brent's method for finding the top of the ripple in the case of Chebyshev Type I band-pass filters?

> xxBandStop::PassbandHint(): Isn't the magnitude the
> same at 0 and pi anyway?

Yes but just in case there are artifacts caused by extreme parameter values and floating point quantization, I choose the side that is farther away from the roll off area.

> And BrentHint(), all versions: I suspect
> it's superfluous to contract the interval by some delta as you
> do.

I did that just in case the evaluation of the magnitude response at 0 and PI is undefined.

Once again thanks for the tips.