Fit biquad to 5 frequency response points?

On Friday, December 2, 2016 at 9:01:32 AM UTC-5, Greg Berchin wrote:
> On Friday, December 2, 2016 at 7:31:18 AM UTC-6, thejoh...@gmail.com wrote:
> > Is it possible to fit a biquad filter to a curve given only 5 (complex) frequency response values?
> > 
> > I'm a relative newcomer. Any help would be greatly appreciated.
> 
> In general, at comp.dsp we don't solve homework problems for you, but since you asked so nicely I'll give you some hints.
> 
> 1. You have 5 unknowns; a1, a2, b0, b1, b2. You have 5 measurements. You have enough information to create 5 equations in 5 unknowns.
> 2. Remember that z = exp(jwT) when dealing with frequency response on the unit circle.
> 3. Use Euler's Identity.
> 4. You have an expression for the transfer function H(z) in terms of z. You know the values of H(z) at several frequencies. You can compute the values of z are at those same frequencies.
> 
> Turn the crank.

after Euler's Identity and step 4, you have to compute the magnitude |H(e^(jwT)| to fit 5 *magnitude* measurements.  BTW, when you do that for a biquad you get the equation at the bottom of this: http://dsp.stackexchange.com/questions/16885/how-do-i-manually-plot-the-frequency-response-of-a-bandpass-butterworth-filter-i/16911#16911 

this is a *real*-valued equation with 5 unknowns.  fit that to 5 constraints (but 5 arbitrary points sound to me to be 10 constraints) and perhaps you can solve for the unknown coefficients, but they might not be real and, even if they are real, you might not have stability (but you should be able to reflect the unstable poles to inside the unit circle and get the same magnitude).

there are 5 constraints in this EQ design problem: http://dsp.stackexchange.com/questions/19225/audio-eq-cookbook-without-frequency-warping/19253#19253 it is not 5 arbitrary points.  there is one arbitary point (the peak frequency and gain, that's two constraints), the gain at DC and the gain at Nyquist.  that's 4 constraints.  when you add a Q constraint to that (and a definition of Q that everyone accepts), you have a full specification for the 5 biquad coefficients.

r b-j

Thanks guys.

FYI closed form because my application will be a realtime piece of audio software, so the target curve will be changing. (However I'm looking at the Maxima expressions and they're verrry long: https://dl.dropboxusercontent.com/u/9600561/biquadfit.m Some repetitions could be removed there.)

Not to mention unstable results for my inputs.

Is this why people tend to match by magnitude only, then choose the minimum phase solution?

Cheers,
John.

On Friday, December 2, 2016 at 7:13:57 PM UTC, Tim Wescott wrote:
> On Fri, 02 Dec 2016 13:39:39 -0500, Phil Martel wrote:
> 
> > On 12/2/2016 13:10, Tim Wescott wrote:
> >> On Fri, 02 Dec 2016 06:01:24 -0800, Greg Berchin wrote:
> >>
> >>> On Friday, December 2, 2016 at 7:31:18 AM UTC-6, thejoh...@gmail.com
> >>> wrote:
> >>>> Is it possible to fit a biquad filter to a curve given only 5
> >>>> (complex)
> >>>> frequency response values?
> >>>>
> >>>> I'm a relative newcomer. Any help would be greatly appreciated.
> >>>
> >>> In general, at comp.dsp we don't solve homework problems for you, but
> >>> since you asked so nicely I'll give you some hints.
> >>>
> >>> 1. You have 5 unknowns; a1, a2, b0, b1, b2. You have 5 measurements.
> >>> You have enough information to create 5 equations in 5 unknowns.
> >>> 2. Remember that z = exp(jwT) when dealing with frequency response on
> >>> the unit circle.
> >>> 3. Use Euler's Identity.
> >>> 4. You have an expression for the transfer function H(z) in terms of
> >>> z. You know the values of H(z) at several frequencies. You can compute
> >>> the values of z are at those same frequencies.
> >>>
> >>> Turn the crank.
> >>
> >> Given five complex frequency-response values and five real coefficient
> >> values, you have _ten_ equations in five unknowns.  That may be a
> >> tougher nut to crack.
> >>
> >> I don't know where the OP is trying to go with this, but sometimes even
> >> when there is a closed-form solution, it's better to just numerically
> >> optimize.  As a for-instance: when was the last time that you factored
> >> a 5-th order polynomial using the formula?  Do you even know where to
> >> FIND the formula?  I think the CRC handbook only goes up to 3rd order.
> >>
> > In general, you can't factor a 5th order polynomial.  I'm not sure about
> > CRC (mine is from '67), but Maxima will crank out the solutions for
> > general 3rd and 4th order polynomials.  They're huge, messy expressions.
> 
> I thought 5th was the highest you could crack, not the lowest that you 
> couldn't.  My bad.
> 
> -- 
> 
> Tim Wescott
> Wescott Design Services
> http://www.wescottdesign.com
> 
> I'm looking for work -- see my website!

On 2.12.16 21:13, Tim Wescott wrote:
> On Fri, 02 Dec 2016 13:39:39 -0500, Phil Martel wrote:
>
>> On 12/2/2016 13:10, Tim Wescott wrote:
>>> On Fri, 02 Dec 2016 06:01:24 -0800, Greg Berchin wrote:
>>>
>>>> On Friday, December 2, 2016 at 7:31:18 AM UTC-6, thejoh...@gmail.com
>>>> wrote:
>>>>> Is it possible to fit a biquad filter to a curve given only 5
>>>>> (complex)
>>>>> frequency response values?
>>>>>
>>>>> I'm a relative newcomer. Any help would be greatly appreciated.
>>>>
>>>> In general, at comp.dsp we don't solve homework problems for you, but
>>>> since you asked so nicely I'll give you some hints.
>>>>
>>>> 1. You have 5 unknowns; a1, a2, b0, b1, b2. You have 5 measurements.
>>>> You have enough information to create 5 equations in 5 unknowns.
>>>> 2. Remember that z = exp(jwT) when dealing with frequency response on
>>>> the unit circle.
>>>> 3. Use Euler's Identity.
>>>> 4. You have an expression for the transfer function H(z) in terms of
>>>> z. You know the values of H(z) at several frequencies. You can compute
>>>> the values of z are at those same frequencies.
>>>>
>>>> Turn the crank.
>>>
>>> Given five complex frequency-response values and five real coefficient
>>> values, you have _ten_ equations in five unknowns.  That may be a
>>> tougher nut to crack.
>>>
>>> I don't know where the OP is trying to go with this, but sometimes even
>>> when there is a closed-form solution, it's better to just numerically
>>> optimize.  As a for-instance: when was the last time that you factored
>>> a 5-th order polynomial using the formula?  Do you even know where to
>>> FIND the formula?  I think the CRC handbook only goes up to 3rd order.
>>>
>> In general, you can't factor a 5th order polynomial.  I'm not sure about
>> CRC (mine is from '67), but Maxima will crank out the solutions for
>> general 3rd and 4th order polynomials.  They're huge, messy expressions.
>
> I thought 5th was the highest you could crack, not the lowest that you
> couldn't.  My bad.


My grand-uncle, professor Kalle Väisälä catalogued in his PhD 
disertation all algebraically solvable 5th order equations nearly 100 
years ago (1917). He also proved that they are all, so that there is a 
great number of 5th order equations that are not solvable in closed form.

-- 

-TV

On Friday, December 2, 2016 at 12:11:00 PM UTC-6, Tim Wescott wrote:

> Given five complex frequency-response values and five real coefficient=20
> values, you have _ten_ equations in five unknowns.  That may be a tougher=
=20
> nut to crack.

Hmmm ... Tim, you are correct. I didn't think about it long enough. So the =
algebraic problem is overspecified.

BTW; FDLS solves this problem in a least-squares sense, but uses a differen=
t method. FDLS assumes only causality. (The most common formulation of FDLS=
 also assumes real coefficients, but it's not absolutely necessary.) FDLS t=
hen finds the difference equation corresponding to the transfer function, r=
econstructs the input and output waveforms, and uses them to fill the matri=
ces.

The roots of a fourth-order polynomial can be found via closed form.
This polynomial has five coefficients, and so (usually) corresponds
to a problem where one is fitting to five datapoints.

e.g.

One datapoint -- polynomial is a constant (0th order)
Two datapoints -- polynomial is a line (1st order)
etc.

Tangentially -- over a Galois field, there is also a closed form
for fifth-order polynomials.  I might have this information in
my notes from way back.


Steve

On Friday, December 2, 2016 at 2:28:20 PM UTC-6, Greg Berchin wrote:

> Hmmm ... Tim, you are correct. I didn't think about it long enough. So the algebraic problem is overspecified.

I thought about this some more, and now I'm not certain whether the problem is overspecified or redundantly-specified.

If the coefficients are constrained to be real, then the constraint that the poles and zeroes have conjugate symmetry is redundant; there are 5 real unknowns and 5 measurements.

If the coefficients are allowed to be complex, then the constraint that the poles and zeroes have conjugate symmetry causes the imaginary parts of the coefficients to be zero; there are 10 unknowns (5 coefficients, each with a real part and an imaginary part) and 10 measurements.

Hi all, yes, you clearly know a lot more about the theory than me. This is helping me learn!

Through a trial-and-error approach I chose conjugate symmetric match points as I wanted to end up with real coefficients (I thought this was the way it needed to be done). In my implementation this seems to work, all imaginary parts are zero.

How else does one constrain the coefficients to be real?

Thanks, John.



> I thought about this some more, and now I'm not certain whether the problem is overspecified or redundantly-specified.
> 
> If the coefficients are constrained to be real, then the constraint that the poles and zeroes have conjugate symmetry is redundant; there are 5 real unknowns and 5 measurements.
> 
> If the coefficients are allowed to be complex, then the constraint that the poles and zeroes have conjugate symmetry causes the imaginary parts of the coefficients to be zero; there are 10 unknowns (5 coefficients, each with a real part and an imaginary part) and 10 measurements.

On Saturday, December 3, 2016 at 9:10:02 AM UTC-6, thejoh...@gmail.com wrote:

> How else does one constrain the coefficients to be real?

There are a couple of interrelated concepts at work here.

Probably back in algebra class, we learned that polynomials with real coefficients have either real roots or conjugate-symmetric roots. And there's an old method by which one can determine constraints upon the number of positive and negative real roots by counting the number of transitions of the algebraic signs of the coefficients of the polynomial.

In signal processing we learn that a real-valued impulse response has a frequency response that is conjugate-symmetric about DC. If it is a discrete-time impulse response, then its frequency response is also symmetric about Fs/2 because it is periodic. 

Can we say that a real-valued discrete-time impulse response can only be generated by a z-domain difference equation that has only real coefficients? I can't say so for certain off the top of my head, but it seems intuitively satisfying. And a z-domain difference equation can be mapped directly to a z-domain transfer function by applying a causality constraint, so if the difference equation has real coefficients, then so does the transfer function.

On Sat, 03 Dec 2016 05:42:01 -0800, Greg Berchin wrote:

> On Friday, December 2, 2016 at 2:28:20 PM UTC-6, Greg Berchin wrote:
> 
>> Hmmm ... Tim, you are correct. I didn't think about it long enough. So
>> the algebraic problem is overspecified.
> 
> I thought about this some more, and now I'm not certain whether the
> problem is overspecified or redundantly-specified.
> 
> If the coefficients are constrained to be real, then the constraint that
> the poles and zeroes have conjugate symmetry is redundant; there are 5
> real unknowns and 5 measurements.
> 
> If the coefficients are allowed to be complex, then the constraint that
> the poles and zeroes have conjugate symmetry causes the imaginary parts
> of the coefficients to be zero; there are 10 unknowns (5 coefficients,
> each with a real part and an imaginary part) and 10 measurements.

Hmm.

Going back to the original problem, the OP wants to specify a biquad 
(with, presumably, five available coefficients) at five complex frequency 
points.

You've got five degrees of freedom in specifying a single-stage biquad 
with real coefficients (one master gain, two poles and two zeros).  But 
five complex frequency points have ten degrees of freedom.

You could expand this to two biquads and a single-order stage, and then 
you'd have 11 degrees of freedom (one master gain, five poles, and five 
zeros).  That might be better.

However, I'm not sure even this is a good approach, because the fit to 
the curve in between those specified points would be totally 
unconstrained, and there's no guarantee that Really Bad Things wouldn't 
happen (at worst, the problem as stated does nothing to guarantee a 
stable filter, and even if that doesn't happen the frequency response 
could be wildly different from what's desired).

I think the OP needs to step back a bit and ask what he _really_ wants to 
accomplish (equalization, effects, what?) and then go forward along some 
other path.

-- 
Tim Wescott
Control systems, embedded software and circuit design
I'm looking for work!  See my website if you're interested
http://www.wescottdesign.com

> there's no guarantee that Really Bad Things wouldn't 
> happen

Indeed this is what I've discovered!

> I think the OP needs to step back a bit and ask what he _really_ wants to 
> accomplish

I have an audio EQ filter that deviates from some ideal.

This deviation is gentle but arbitrary. I have an mag/phase expression for the deviation (so I can take 'match point' measurements anywhere I want).

I'm looking to correct this deviation, yielding a filter that matches the ideal.

If I go down the FIR route, what options do I have for matching arbitrary responses? Both magnitude and phase?

Ultimately this will be realtime parameter controlled C++ software.

Thanks all, for the enlightening discussion!
John.