Hello,

there exist quite a few methods to design an IIR filter for an arbitrary
frequency response with magnitude and phase. 
For example, there is invfreqz in Matlab (based on the so-called 'Levi's
method'), Prony's method, Steiglitz-McBride, Yule-Walker, to name a few.


Now what do I do if I only know the magnitude response, but not the phase?


Common wisdom says: make it minimum phase, for example in [1], slide 53.
But, does this still hold when numerator and denominator order are
different? For example, consider an IIR filter with a 100 tap numerator and
denominator order 1. I'd expect that it behaves essentially like an FIR
filter, where linear phase would be the best choice.
Also, I don't see how to convert a filter specification to minimum phase,
where my frequency response is weighted and has don't-care areas.

[1] also proposes to iterate the phase (slide 54), use the phase of the
previous design round as target. But I believe I have at least one example
where this approach gets stuck in a local minimum. Some kind of grid search
would be manageable, but it adds one iteration loop, and maybe there is a
better way?

There is no shortage of papers on IIR design methods for magnitude
response. For example, [2] looks promising. Reference 33 in [3] (Steiglitz)
reportedly addresses the problem. The final Matlab program in [4] should do
what I need but didn't converge for high numerator orders, and the math
isn't exactly straighforward.

My question is, are there any known "shortcuts" to handle the unknown phase
response, or maybe any generic method that's known to be robust and
reliable, and ideally can be bent to work on non-standard design problems?

Cheers

Markus

[1] http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/00900_OptimalIIR.pdf
[2]
http://www.mendeley.com/research/frequencydomain-steiglitz-mcbride-method-leastsquares-iir-filte
r-design-arma-modeling-periodogram-smoothing/#page-1
[3] http://cnx.org/content/m41681/1.2/
[4] http://cnx.org/content/m41690/latest/?collection=col11383/latest

On Sun, 24 Jun 2012 10:56:07 -0500, mnentwig wrote:

> Hello,
> 
> there exist quite a few methods to design an IIR filter for an arbitrary
> frequency response with magnitude and phase.
> For example, there is invfreqz in Matlab (based on the so-called 'Levi's
> method'), Prony's method, Steiglitz-McBride, Yule-Walker, to name a few.
> 
> 
> Now what do I do if I only know the magnitude response, but not the
> phase?
> 
> 
> Common wisdom says: make it minimum phase, for example in [1], slide 53.
> But, does this still hold when numerator and denominator order are
> different? For example, consider an IIR filter with a 100 tap numerator
> and denominator order 1. I'd expect that it behaves essentially like an
> FIR filter, where linear phase would be the best choice.
> Also, I don't see how to convert a filter specification to minimum
> phase,
> where my frequency response is weighted and has don't-care areas.

The magnitude vs. frequency response is the same whether the filter is 
minimum or non-minimum phase, so keeping the thing minimum phase should 
simplify computation.

And keeping an IIR filter guarantees that it is stable (or altogether 
unstable), because a minimum phase filter has all its poles and zeros on 
the same side of the stability boundary.

> [1] also proposes to iterate the phase (slide 54), use the phase of the
> previous design round as target. But I believe I have at least one
> example where this approach gets stuck in a local minimum. Some kind of
> grid search would be manageable, but it adds one iteration loop, and
> maybe there is a better way?

If I were going to use an algorithm that depended on phase, and I had a 
complete magnitude response, I would use the Hilbert transform.

> There is no shortage of papers on IIR design methods for magnitude
> response. For example, [2] looks promising. Reference 33 in [3]
> (Steiglitz)
> reportedly addresses the problem. The final Matlab program in [4] should
> do what I need but didn't converge for high numerator orders, and the
> math isn't exactly straighforward.

Goodness, that is complex.  Is it trying to find the whole filter as a 
ratio of polynomials?  You don't say what "high" numerator orders are, 
but if you mean more than ten or so then it may simply be running afoul 
of polynomial coefficient sensitivity issues.

> My question is, are there any known "shortcuts" to handle the unknown
> phase response, or maybe any generic method that's known to be robust
> and reliable, and ideally can be bent to work on non-standard design
> problems?

Here's what I'd try.  It's rough and ready, which means that it'll never 
find it's way into a hoity-toity journal:

Use some dumb-thumb general-purpose optimization program.  Scilab comes 
with gradient and genetic optimization built in, and I think they'll both 
do constrained optimisation with little effort.  I assume that Matlab has 
something, too.

Express your filter as a cascade of 2nd-order IIR filters which are 
described in terms of pole and zero angle and magnitude (or natural 
frequency and damping, or some other 'natural' means, but only two for 
each).  Use _one_ overall gain for the whole filter.

If your overall filter is supposed to be lowpass, constrain each filter 
section to have a DC gain of one.  If your overall filter is supposed to 
be bandpass, constrain each filter section to be bandpass, with a peak 
gain of 1 (some cleverness may be required here, but do this as a first 
try).

Choose your optimization critera to match your desires (i.e., if you need 
closest approach or "do not exceed" or whatever).

Etc.

Then just plug the whole thing into your blind optimization routine, go 
get a good night's sleep, and in the morning see what it's coughed up.

> Cheers

Cheers to you, too -- I hope this helps.

> Markus
> 
> [1]
> http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/00900_OptimalIIR.pdf
> [2]
> http://www.mendeley.com/research/frequencydomain-steiglitz-mcbride-
method-leastsquares-iir-filter-design-arma-modeling-periodogram-smoothing/
#page-1
> [3] http://cnx.org/content/m41681/1.2/
> [4] http://cnx.org/content/m41690/latest/?collection=col11383/latest





-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com

Tim, 

thanks for the comments. 
I should mention that this is not related to any "real-life" problem - at
least, not yet. 

>> and I had a complete magnitude response
I think that's the catch. The interesting part - transition from
pass-to-stopband is unknown.

I've tried minimum phase using both Hilbert transform and reflection of the
cepstrum for linear phase, and alternatively linear phase. The
minimum-phase starting point doesn't get even close to the filter resulting
from a linear phase starting point.
It could be a problem in my implementation, but I suspect it's conceptually
a different problem from (for example) system identification, where the
full frequency response is known, and the numerator order is usually lower
or equal to the denominator.

>> Is it trying to find the whole filter as a ratio of polynomials?

well... I had only a quick look what it does exactly. I think it linearizes
the equation-error problem into a solution error problem. Details can be
found here, page 82:
http://scholarship.rice.edu/handle/1911/16483

One numerical issue I've spotted with the codes on the page is the
ramping-up of the order to get a minimax solution. I understand "desired P
= 400" means it ultimately tries to calculate something like |H-D|^400, not
good. 

And yes, a general-purpose numerical solver like "fminsearch" (Matlab) or
"fmincg" should probably do the job. Those could indeed run over night, and
if needed, I could even get partial derivatives with Maxima CAS without too
much hassle. 

-markus

PS: If anybody is interested, here is a link to the Octave code I've got so
far (probably works on Matlab, too). 
This is by no means production code, just something I've hacked up over the
weekend. Expect bugs.
The problem is arbitrary (6th order recursive and 21st order non-recursive
filter. Read it a 3x2nd order IIR cascaded with a 15th order FIR, and
re-factor the resulting numerator with roots(...); poly(...) ).
http://www.dsprelated.com/blogimages/MarkusNentwig/comp.dsp/IIRtest2.m

On 6/24/12 11:56 AM, mnentwig wrote:
>
> there exist quite a few methods to design an IIR filter for an arbitrary
> frequency response with magnitude and phase.
> For example, there is invfreqz in Matlab (based on the so-called 'Levi's
> method'), Prony's method, Steiglitz-McBride, Yule-Walker, to name a few.

don't forget Greg Bechin's FDLS.


-- 

r b-j                  r...@audioimagination.com

"Imagination is more important than knowledge."

>> don't forget Greg Bechin's FDLS.
Thanks, I'll need to look that up.

One solution that doesn't look too complicated is here:
http://www.nt.tuwien.ac.at/fileadmin/users/gerhard/diss_Lang.pdf
section 2.1.4, least-squares with independent weighting for magnitude and
phase error. 
It's written for FIR, and I wonder whether the same would still work with
IIR.

Another practical way could be to design an initial filter for linear
phase, then look at the numerator. If non-minimum zeros are found, reflect
them, calculate the phase, set as target phase and redesign. 

On Jun 24, 2:37 pm, robert bristow-johnson <r...@audioimagination.com>
wrote:

> don't forget Greg Bechin's FDLS.

FDLS works really well for non-arbitrary frequency responses, that is,
frequency responses that are known to have originated from either
Laplace-domain transfer functions or z-domain transfer functions. But
Vladimir experimented with arbitrary frequency responses and found
that even constraining them to minimum-phase was not sufficient to
guarantee good convergence. So as a frequency response modeler, FDLS
is good. But as a frequency response synthesizer, it is not so good.

Greg

>> But Vladimir experimented with arbitrary frequency responses and found
>> that even constraining them to minimum-phase was not sufficient to
>> guarantee good convergence.

Hello,

maybe I'm just running into the same problems. But, I suspect the problem
is not FDLS as such but misuse of the "minimum phase" concept when I've got
only a vague idea on how my frequency response is going to look.

I know at least of one other filter design program, "lslevin" in
http://www.nt.tuwien.ac.at/fileadmin/users/gerhard/diss_Lang.pdf page 20
that also uses cosine waves (via exp() and real part) on an arbitrary
frequency grid. It works great for FIR filter _synthesis_, if this means
anything. 
I reverse-engineered the program here, if anybody is interested: Follow the
download link at the bottom of
http://www.dsprelated.com/showarticle/165.php , the function run_LS() is
the relevant part). 

> Hello,
>
> there exist quite a few methods to design an IIR filter for an arbitrary
> frequency response with magnitude and phase.
> For example, there is invfreqz in Matlab (based on the so-called 'Levi's
> method'), Prony's method, Steiglitz-McBride, Yule-Walker, to name a few.
>
> Now what do I do if I only know the magnitude response, but not the phase?

1. Find frequency area with maximum deviation of desired response from 
response that you currently have.
2. Add biquad section designed to rectify that area.
3. Make an optimization pass.
4. Repeat 1-2-3 until all sections are used.

If LS approximation is suitable, then generate noise with desired frequency 
response by FFT, then design ARMA filter to that noise.

 Vladimir Vassilevsky
 DSP and Mixed Signal Consultant
 www.abvolt.com

>   Vladimir Vassilevsky
>   DSP and Mixed Signal Consultant
>   www.abvolt.com


nice to see you around Vlad.


-- 

r b-j                  r...@audioimagination.com

"Imagination is more important than knowledge."