f



Fit biquad to 5 frequency response points?

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.

I have a gentle curve which I am looking to find a closed form solution for the coefficients in terms of the 5 known frequency response values.

More detailed problem description is here: https://dl.dropboxusercontent.com/u/9600561/question.pdf

The best attempt I've had is to use Mathematica's Solve which has given me something like...
https://dl.dropboxusercontent.com/u/9600561/biquadfit.m
....however it gives a division by zero and NaN's.

Many thanks in advance for your time,
John.
0
thejohnflynn
12/2/2016 1:31:15 PM
comp.dsp 20333 articles. 1 followers. allnor (8509) is leader. Post Follow

24 Replies
248 Views

Similar Articles

[PageSpeed] 57

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.
0
Greg
12/2/2016 2:01:24 PM
Thanks for the swift reply Greg!

Actually it's not homework, I typed up that pdf to make the question as clear as possible. I guess this is what I do for fun.

I am on the right track, I think...

I have these five equations:
https://dl.dropboxusercontent.com/u/9600561/q2.pdf

Using Mathematica I can rearrange the equation for appropriate back-substitution which has given me:
https://dl.dropboxusercontent.com/u/9600561/biquadfit.m

Which doesn't work :(

When you say 'turn the crank', is there something easier that I'm missing??

Are there pen and paper simplifications I should be making first?

Many thanks again, John.




On Friday, December 2, 2016 at 2:01:32 PM UTC, 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.
0
thejohnflynn
12/2/2016 2:23:43 PM
On Friday, December 2, 2016 at 8:23:45 AM UTC-6, thejoh...@gmail.com wrote:

> When you say 'turn the crank', is there something easier that I'm missing??
> 
> Are there pen and paper simplifications I should be making first?

Use Euler's Identity to resolve the "exp(3pi/4)" terms. After that, it's algebra.
0
Greg
12/2/2016 3:09:47 PM
On Friday, December 2, 2016 at 3:09:52 PM UTC, Greg Berchin wrote:
> On Friday, December 2, 2016 at 8:23:45 AM UTC-6, thejoh...@gmail.com wrote:
> 
> > When you say 'turn the crank', is there something easier that I'm missing??
> > 
> > Are there pen and paper simplifications I should be making first?
> 
> Use Euler's Identity to resolve the "exp(3pi/4)" terms. After that, it's algebra.

Thanks Greg,

Yes for these I have the complex value -1/sqrt(2) + j1/sqrt(2) or -0.70711 + 0.70711i. So yes I know these e^ terms are 'knowns'.

My question then, is more about the algebra...

Do I just start grinding with the DC equation to get

b0 = -b2 - b1 + H1 * a2 + H1 * a1 + H1

Then substitute "-b2 - b1 + H1 * a2 + H1 * a1 + H1" for every b0 in the next equation.
Rearrange that new equation so b1 = ...some complicated expression...
Rearrange that new equation so b2 = ...some even more complicated expression...
Rearrange...
....
And on till the end?

Is there an easier way?!

Sorry if this is really obvious, I've been away from maths for 17 years so rusty isn't the word.

Thanks!
0
thejohnflynn
12/2/2016 3:51:50 PM
On Friday, December 2, 2016 at 9:51:53 AM UTC-6, thejoh...@gmail.com wrote:

> Is there an easier way?!

Try putting the equations into matrix form.
0
Greg
12/2/2016 4:32:40 PM
Embarrassing, my code actually did work :/

Typo in the call...
biquadfit(hDc, hPi2, h3Pi4, h3Pi4Neg, hPi2)
....should have been...
biquadfit(hDc, hPi2, h3Pi4, h3Pi4Neg, hPi2Neg) <--

Unfortunately I end up with a pole outside the unit circle and only the magnitude but not the phase is matched.

Is the solution unique?

Since the unknowns are not squared/cubed/etc., I would have presumed that there is only one solution.

I'm wondering why the phase is not matched...

Thanks for all your input, it's been very helpful,
John.


0
thejohnflynn
12/2/2016 6:06:04 PM
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.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

I'm looking for work -- see my website!
0
Tim
12/2/2016 6:10:52 PM
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.

    Best wishes,
    --Phil

-- 
Best wishes,
--Phil
pomartel At Comcast(ignore_this) dot net
0
Phil
12/2/2016 6:39:39 PM
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!
0
Tim
12/2/2016 7:13:48 PM
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 wrot=
e:
> > Is it possible to fit a biquad filter to a curve given only 5 (complex)=
 frequency response values?
> >=20
> > I'm a relative newcomer. Any help would be greatly appreciated.
>=20
> In general, at comp.dsp we don't solve homework problems for you, but sin=
ce you asked so nicely I'll give you some hints.
>=20
> 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 =3D exp(jwT) when dealing with frequency response on t=
he unit circle.
> 3. Use Euler's Identity.
> 4. You have an expression for the transfer function H(z) in terms of z. Y=
ou know the values of H(z) at several frequencies. You can compute the valu=
es of z are at those same frequencies.
>=20
> 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 biqua=
d 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=20

this is a *real*-valued equation with 5 unknowns.  fit that to 5 constraint=
s (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, ev=
en if they are real, you might not have stability (but you should be able t=
o reflect the unstable poles to inside the unit circle and get the same mag=
nitude).

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

r b-j
0
robert
12/2/2016 7:50:55 PM
Thanks guys.

FYI closed form because my application will be a realtime piece of audio so=
ftware, so the target curve will be changing. (However I'm looking at the M=
axima expressions and they're verrry long: https://dl.dropboxusercontent.co=
m/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:
>=20
> > 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 =3D 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 comput=
e
> >>> 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 eve=
n
> >> 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 abou=
t
> > CRC (mine is from '67), but Maxima will crank out the solutions for
> > general 3rd and 4th order polynomials.  They're huge, messy expressions=
..
>=20
> I thought 5th was the highest you could crack, not the lowest that you=20
> couldn't.  My bad.
>=20
> --=20
>=20
> Tim Wescott
> Wescott Design Services
> http://www.wescottdesign.com
>=20
> I'm looking for work -- see my website!
0
thejohnflynn
12/2/2016 7:53:46 PM
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

0
Tauno
12/2/2016 7:54:47 PM
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.
0
Greg
12/2/2016 8:28:17 PM
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
0
spope33
12/2/2016 9:29:59 PM
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 th=
e 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 th=
e poles and zeroes have conjugate symmetry is redundant; there are 5 real u=
nknowns 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.
0
Greg
12/3/2016 1:42:01 PM
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 i=
t needed to be done). In my implementation this seems to work, all imaginar=
y 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 probl=
em is overspecified or redundantly-specified.
>=20
> 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.
>=20
> If the coefficients are allowed to be complex, then the constraint that t=
he poles and zeroes have conjugate symmetry causes the imaginary parts of t=
he coefficients to be zero; there are 10 unknowns (5 coefficients, each wit=
h a real part and an imaginary part) and 10 measurements.

0
thejohnflynn
12/3/2016 3:09:59 PM
On Saturday, December 3, 2016 at 9:10:02 AM UTC-6, thejoh...@gmail.com wrot=
e:

> 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 coeff=
icients have either real roots or conjugate-symmetric roots. And there's an=
 old method by which one can determine constraints upon the number of posit=
ive and negative real roots by counting the number of transitions of the al=
gebraic signs of the coefficients of the polynomial.

In signal processing we learn that a real-valued impulse response has a fre=
quency response that is conjugate-symmetric about DC. If it is a discrete-t=
ime impulse response, then its frequency response is also symmetric about F=
s/2 because it is periodic.=20

Can we say that a real-valued discrete-time impulse response can only be ge=
nerated 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 d=
ifference equation has real coefficients, then so does the transfer functio=
n.
0
Greg
12/3/2016 3:55:21 PM
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
0
Tim
12/3/2016 10:14:45 PM
> 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.

0
thejohnflynn
12/4/2016 11:18:01 AM
On Sun, 4 Dec 2016 03:18:01 -0800 (PST), thejohnflynn@gmail.com wrote:

>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?

If the deviation is minimum phase, then you can compute an inverse
transfer function that is stable. Of course, you have to have the
transfer function (poles and zeroes) before you can invert it. Note that
in general this will not be a FIR filter.

If the deviation is nonminimum phase, then you can use a FIR filter to
approximate the inverse, but you will introduce delay in the process.
That delay may or may not be a problem for you.

Another possibility is the aforementioned FDLS (Frequency Domain Least
Squares) approximation technique. There are no FDLS constraints upon
stability, but if the "deviation inverse" frequency response turns out
to be unstable, then the approximation will be unstable, too. Oh, for
FDLS the "deviation inverse" must also be causal.

Send me a PM if you want FDLS info and sample Matlab or C program.
Change "chatter" to "charter" and remove "invalid". Others are welcome
to do the same.
0
Greg
12/4/2016 3:49:32 PM
Greg, it's very reassuring to leave this forum, study the literature and experiment, come back to a post like this and see that I've been working along the right lines. Thank you.

> Of course, you have to have the
> transfer function you can invert it.

Apologies, I don't have the transfer function of the error curve. The ideal filter is in the S domain.

> Another possibility is the aforementioned FDLS (Frequency Domain Least
> Squares) approximation technique.

Now this is really helping.

Using FDLS with 8 measurements to give a 7 zero correction filter I get _extremely_ good results. (-60dB magnitude error. And very useable results with lower settings).

Are there any other FIR methods I should be looking into?

For low-zero correction FDLS seems to have it covered. But also, for example if I wanted to do an 'total overkill' match. 64, 128, 256 zeros and higher.

Is there a method where I could 'pin down' the response exactly at 64, evenly spaced points?

(That also behaves pretty well between the points).

Thanks,
John.
0
thejohnflynn
12/4/2016 6:50:23 PM
On Sunday, December 4, 2016 at 12:50:28 PM UTC-6, thejoh...@gmail.com wrote=
:

> Using FDLS with 8 measurements to give a 7 zero correction filter I get _=
extremely_ good results. (-60dB magnitude error. And very useable results w=
ith lower settings).

Understand that FDLS is more typically, and more preferably, used with hund=
reds or thousands of frequency response measurements. In addition to MA (FI=
R) models, it can create AR (all pole) models or ARMA (IIR) models. In fact=
, ARMA is the most commonly used configuration.

> Are there any other FIR methods I should be looking into?
>=20
> For low-zero correction FDLS seems to have it covered. But also, for exam=
ple if I wanted to do an 'total overkill' match. 64, 128, 256 zeros and hig=
her.
>=20
> Is there a method where I could 'pin down' the response exactly at 64, ev=
enly spaced points?

Yeah; the inverse DFT.
0
Greg
12/4/2016 11:37:54 PM
On Sunday, December 4, 2016 at 6:18:09 AM UTC-5, thejoh...@gmail.com wrote:
> 
> I have an audio EQ filter that deviates from some ideal.

what is the ideal?  do you have a prototype EQ or target response defined?

....

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

if causality isn't needed, you can do whatever magnitude and phase you want.  Eric Jacobsen suggested a quick and dirty idea in the previous millennium. http://dspguru.com/dsp/tricks/using-parks-mcclellan-to-design-non-linear-phase-fir-filter 

to make it causal, add delay.

r b-j
0
robert
12/5/2016 6:13:42 AM
On Sun, 04 Dec 2016 03:18:01 -0800, thejohnflynn wrote:

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

Assuming that you can add arbitrary amounts of delay to your signal, the 
brute force way to do this is to do an inverse DFT (which ends up being 
an inverse FFT in most reasonable systems -- an FFT is just a fast DFT).  
It may not be the shortest-time way to do the job, but it'll always take 
the same amount of time, which is desirable in a real-time system.

The algorithm becomes: (1) determine your desired frequency response, (2) 
window it, symmetrically around f = 0, to limit the extent of the time 
response, (3) do an IFFT, (4) shift the IFFT so that your filter is 
causal, (5) use the result as the set of taps in an FIR filter.

The hardest part will be automagically figuring out the proper amount of 
shift, and making sure you have enough taps so that the impulse response 
has mostly died off before the end of the FIR (if you don't do this, then 
the tail end of the filter response will bleed into the beginning, 
because of the nature of the DFT).

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

I'm looking for work -- see my website!
0
Tim
12/5/2016 6:24:48 PM
Reply: