Minimum Phase Impulse Response

Does anyone know how to find the precise minimum phase discrete-time impulse
response corresponding to a given magnitude spectrum, when the magnitude
spectrum does *not* correspond to rational polynomial transfer function?

Matt:

[snip]
"Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in message
news:edZnb.2748$Nz5.331406@news20.bellglobal.com...
> Does anyone know how to find the precise minimum phase discrete-time
impulse
> response corresponding to a given magnitude spectrum, when the magnitude
> spectrum does *not* correspond to rational polynomial transfer function?
[snip]

The magnitude and phase of a minimum phase function are related by the
Hilbert Transform Integrals
[if they exist].  Depending upon your function you may be able to compute
this using numerical quadrature.

--
Peter
Consultant
Indialantic By-the-Sea, FL
[peter dot brackett at ieee dot org]

"Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in message
news:edZnb.2748$Nz5.331406@news20.bellglobal.com...
> Does anyone know how to find the precise minimum phase discrete-time
impulse
> response corresponding to a given magnitude spectrum, when the magnitude
> spectrum does *not* correspond to rational polynomial transfer function?

Matt,

Hmmmmm.... working backwards and using a counter-example:
Starting with the magnitude spectrum
Calculate the magnitude spectrum squared
Now, assuming that there *was* a rational polynomial transfer function, then
the magnitude spectrum squared will have zeros that are arrayed in conjugate
pairs around the unit circle in the z-plane, and will have poles that are
arrayed in duals inside the unit circle, right?
Then, one of the square roots of this rational function will have
approximately 2x fewer poles and 2x fewer zeros.  The zeros will all lie
inside the unit circle for the minimum phase version.
Well, I'm not so sure about the poles but this is certainly correct for the
zeros.

So far, so good .....

Now, if the magnitude spectrum doesn't have a corresponding family of
rational functions that it would result from, then what the heck kind of
"system" does the impulse response you're looking for represent?
In general, it's not a system that can be defined with a set of ordinary,
linear, difference equations with constant coefficients (OLDEWCC) I do
believe.  Otherwise, the difference equations would result in a rational
transfer function wouldn't it?

So, what part of OLDEWCC isn't met?
Ordinary difference equations -> partial difference equations
Linear difference equations -> nonlinear difference equations
Constant coefficients -> time-varying coefficients

If the coefficients aren't constant, then there is no single impulse
response - there's an infinity of them.

If the equations aren't linear then there is no single impulse response in
the normal sense - because the response doesn't scale with the size of the
"unit sample" amplitude input.  The response varies.

If the difference equations are partial difference equations then the notion
of an impulse response doesn't apply because the equations represent a field
and not a single-input, single-output system.  In that case you might look
at the value of each state variable in response to a particular unit sample
"input" somewhere in the field and you will get "impulse responses" for each
of the variables.

So .... are you *sure* that a rational function doesn't describe your
system?

That much said: perhaps you have a measured magnitude spectrum that doesn't
lend itself to the nice mathematics?  In that case you may want to develop
an approximation that's adequate.  To do that, Google on "system
identification" and "deconvolution" .... that sort of thing.  Once you have
a system defined you will have an impulse response.

Now I'm really curious.  Did I cover all the bases or is there such a thing
as described and I missed how that could be?  Where did the magnitude
spectrum come from??

Fred

Thanks, Fred.

comments interleaved below

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:T_mdnbL_NpwQtjyiRVn-sg@centurytel.net...
> Now, if the magnitude spectrum doesn't have a corresponding family of
> rational functions that it would result from, then what the heck kind of
> "system" does the impulse response you're looking for represent?

Ah well, that's just it -- I don't know.  What I'm trying to do is make a
special type of low-stop filter.  I want to set a constraint that no more
than x% of the signal power can be in the stop band, and then find the
filter that:

a) satisfies this condition; and

b) packs as much energy as possible into the initial samples.

From (b), I can immediately tell that the best filter will be the minimum
phase filter for its magnitude response, but the constraint leaves a *lot*
of freedom in the magnitude response.

There is no reason to assume that the best filter will have a rational
polynomial transfer function.  Perhaps the best frequency response is
derived from the Gaussian, or perhaps the transfer function will have zeros
bigger than points, like a sinc.

> That much said: perhaps you have a measured magnitude spectrum that
doesn't
> lend itself to the nice mathematics?  In that case you may want to develop
> an approximation that's adequate.  To do that, Google on "system
> identification" and "deconvolution" .... that sort of thing.  Once you
have
> a system defined you will have an impulse response.

Unfortunately, I don't know how to determine when an approximation is
adequate.  I can generate any number of minimum phase filters that satisfy
the constraint on the frequency response, but I don't know how to determine
how close they come to the best filter as defined above.

If I knew how to make the minimum phase impulse response for an arbitrary
frequency response, then I could probably find a way to choose the best
filter from the set that satisfies the constraint.  I'm trying to find out
more about the "real" definition of, instead of the one about poles and
zeros.

Perhaps you can help with another question, from which I can derive an
answer to the original one:

Given the entire complex frequency response for a discrete-time signal, how
can you tell whether or not the signal is causal?

Again, I have no rational polynomials.

Matt,

Comments interspersed below:

"Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in message
news:P9kob.4987$Nz5.422060@news20.bellglobal.com...
> Thanks, Fred.
>
> comments interleaved below
>
> "Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
> news:T_mdnbL_NpwQtjyiRVn-sg@centurytel.net...
> > Now, if the magnitude spectrum doesn't have a corresponding family of
> > rational functions that it would result from, then what the heck kind of
> > "system" does the impulse response you're looking for represent?
>
> Ah well, that's just it -- I don't know.  What I'm trying to do is make a
> special type of low-stop filter.  I want to set a constraint that no more
> than x% of the signal power can be in the stop band, and then find the
> filter that:

Please define "low-stop" for me.  Do you mean time domain lowpass lifter.
(lifter means filter but in the time domain)
>
> a) satisfies this condition; and
>
> b) packs as much energy as possible into the initial samples.\\

Well, it appears you want a lowpass filter that is also a lowpass lifter,
right?
I think you will find that an approximation to a Gaussian shape has this
character - being minimum width in both time and frequency at the same time.
Now all you have to do is define a FIR filter that approximates the Gaussian
shape - an the frequency response should do the same.

>
> From (b), I can immediately tell that the best filter will be the minimum
> phase filter for its magnitude response, but the constraint leaves a *lot*
> of freedom in the magnitude response.

I don't know about this one way or another....it seems possible but not a
rigorous statement.

>
> There is no reason to assume that the best filter will have a rational
> polynomial transfer function.  Perhaps the best frequency response is
> derived from the Gaussian, or perhaps the transfer function will have
zeros
> bigger than points, like a sinc.

What do you mean by "zeros bigger than points"?
Yes, there is reason to assume that the best filter will have a rational
transfer function because of all those conditions I enumerated before.  A
filter is a system and is described by ordinary linear difference equations
with constant coefficients.
Now, if you are assuming something else, then perhaps not.  For example, we
can define Fourier transform pairs that don't have anything to do with
"systems".  But, in the end, we are usually interested in systems - and
perhaps systems that approximate ideal Fourier transform pairs.
>
> > That much said: perhaps you have a measured magnitude spectrum that
> doesn't
> > lend itself to the nice mathematics?  In that case you may want to
develop
> > an approximation that's adequate.  To do that, Google on "system
> > identification" and "deconvolution" .... that sort of thing.  Once you
> have
> > a system defined you will have an impulse response.
>
> Unfortunately, I don't know how to determine when an approximation is
> adequate.  I can generate any number of minimum phase filters that satisfy
> the constraint on the frequency response, but I don't know how to
determine
> how close they come to the best filter as defined above.

I don't know why the term "minimum phase" keeps coming up so readily here.
Do you perhaps mean something like "minimum delay" or "maximum integral of
the magnitude squared over a minimum span"?
My recommendation would be to drop the use of that term because it has very
specific meaning or usage - rather than a more general meaning.  Later, when
we agree that we're dealing with system functions and not just abstract
transform pairs, we can talk about minimum phase.

>
> If I knew how to make the minimum phase impulse response for an arbitrary
> frequency response, then I could probably find a way to choose the best
> filter from the set that satisfies the constraint.  I'm trying to find out
> more about the "real" definition of, instead of the one about poles and
> zeros.

Poles and zeros are very "real" in ordinary linear difference equations with
constant coefficients (i.e. "systems")

>
> Perhaps you can help with another question, from which I can derive an
> answer to the original one:
>
> Given the entire complex frequency response for a discrete-time signal,
how
> can you tell whether or not the signal is causal?
>
> Again, I have no rational polynomials.


Well, maybe you do and maybe you don't.  It depends on where you are in the
process of finding a solution to the problem you've posed.  I can understand
that there are optimum functions that aren't system functions but there are
no systems that aren't described by system functions.

First, strictly speaking, there is no such thing as a causal "signal" - only
a causal system response.
A system response is causal if the poles of the frequency response lie
inside the unit circle.
A system is described by an ordinary linear difference equation with
constant coefficients - so it does have a rational system function.
A minimum phase system has all the poles inside the unit circle and all the
zeros inside or on the unit circle.

Take a look at page 199 of Temes, Barcilon and Marshall "The Optimization of
Bandlimited Sytems" Proc IEEE Vol 61 No. 2 Feb 1973 pp 196-234.  Here,
Victor Barcilon discusses time spread and frequency spread, deltaT and
deltaOmega respectively, and shows that the Gaussian is optimum, that
deltaT*deltaOmega>=1/2, etc......  These are optimum functions and not
system functions.
Later in the paper we discuss the realization of system functions - well
things like FIR filters and window functions which are discrete and
continuous versions of the same thing.

Fred

Hi Fred,

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:Fr2dnUZI8ZJShD-iRVn-iA@centurytel.net...
> Please define "low-stop" for me.  Do you mean time domain lowpass lifter.
> (lifter means filter but in the time domain)

It means I just want to attenuate low frequencies, and don't care much about
the high frequencies -- I don't require any limits on passband ripple or
other passband errors.

> > From (b), I can immediately tell that the best filter will be the
minimum
> > phase filter for its magnitude response, but the constraint leaves a
*lot*
> > of freedom in the magnitude response.
>
> I don't know about this one way or another....it seems possible but not a
> rigorous statement.
> [...]
> I don't know why the term "minimum phase" keeps coming up so readily here.
> Do you perhaps mean something like "minimum delay" or "maximum integral of
> the magnitude squared over a minimum span"?

See the final paragraphs of this link.  The characteristic "maximum decay"
result for minimum phase filters is exactly what I'm looking for, and
defines the term without reference to poles or zeros.

http://ccrma-www.stanford.edu/~jos/filters/Definition_Minimum_Phase.html

> Now, if you are assuming something else, then perhaps not.  For example,
we
> can define Fourier transform pairs that don't have anything to do with
> "systems".  But, in the end, we are usually interested in systems - and
> perhaps systems that approximate ideal Fourier transform pairs.

Yes.  What I exepect to get at the end is a long FIR filter that truncates
or approximates an infinite response that wasn't generated by differential
equations.  I want to find the ideal so I can tell when the approximation is
good.

> First, strictly speaking, there is no such thing as a causal "signal" -
only
> a causal system response.

One-sided, if you prefer.

> A system response is causal if the poles of the frequency response lie
> inside the unit circle.

An signal is one-sided if it is zero at t<0

By causal freqency response, I mean exactly the Fourier transform of such a
signal.

Hmmm... Now that I put it that way, surely it can't be too hard to state the
necessary condition.

"Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in message news:<P9kob.4987$Nz5.422060@news20.bellglobal.com>...

> Again, I have no rational polynomials.

Actually, there may be a chance that you do. The Z transform of a FIR filter 
*is* a rational polynomial, though with only a numerator and no denominator.

Rune

Fred Marshall wrote:
   ...



> I think you will find that an approximation to a Gaussian shape has this
> character - being minimum width in both time and frequency at the same time.
> Now all you have to do is define a FIR filter that approximates the Gaussian
> shape - an the frequency response should do the same.

Going off half cocked -- that's "detente" in French -- I think that that
shape is binomial. I mean by that the binomial coefficients on x and y
in (x + y)^n and scaled by 1/2^n, For n = 5, that gives
                    1/32[1, 5, 10, 10, 5, 1]
and so on. I've used these filters before I knew that DSP was a branch
of learning. They can all be derived by repeating the average of two
adjacent values, which is the filter 1/2[1, 1]. Would anyone care to say
that a Pascal filter is the discrete form of an Gaussian? (A la Poisson
and Rayleigh.)

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

allnor@tele.ntnu.no (Rune Allnor) wrote in message news:<f56893ae.0310310733.7b32304f@posting.google.com>...
> "Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in message news:<P9kob.4987$Nz5.422060@news20.bellglobal.com>...
> 
> > Again, I have no rational polynomials.
> 
> Actually, there may be a chance that you do. The Z transform of a FIR filter 
> *is* a rational polynomial, though with only a numerator and no denominator.

I just couldn't let that go, so here's the results of my "sunday morning
mentals"...

We are concerned with a FIR filter with filter coefficients 

   h(n), n=-N/2,...,N/2                                             [1]

and z transform 

            k=N/2
   H(z) =   sum    h(k)z^(-k).                                      [2]
           k=-N/2

We know that a filters on the form 

   H(z) = +/-H(z^(-1))                                              [3]

have linear phase responses, and we know that filters that have 
all their roots inside the unit circle, 

            k=N/2
   H(z)  =  prod    (z-z_k)      |z_k| < 1 for all k                [4]
            k=-N/2

are "minimum phase". 

So the problem is stated as "given the magnitude spectrum of a FIR filter,
how do we find the filter coefficients that are minimum phase."

Since we start out with a magnitude spectrum, I prefer to express it in 
terms of a non-causal linear-phase FIR filter with zero phase lag at all 
frequencies. This is, by [3] a filter with symmetrical impulse response.

Next, I recall the well-known fact that symmetrical FIR filters have 
reciprocal roots, i.e.

  H(z)=G(z)G(z^-1)                                                  [5]

subject to 

  |H(z)| = |G(z)|^2                                               [6.a]
         = |G(z^-1)|^2                                            [6.b]
         = |G(z)G(z^-1)|                                          [6.c]

where G(z) represents the causal part of H(z) and G(z^-1) represents 
the anticausal part of H(z).

By construction, a naive IDFT of the specified magnitude spectrum would 
yield a filter impulse respone that zero-phase, though circular convolved,
representing |H(z)| by [6.c] above.  

So my very naive idea (concieved before 9 AM on a sunday morning, so no 
guarantees aboute any trace of coherence, or even sanity, apply!) is
to try to extract the causal part of the impulse response, and design 
a causal linear phase FIR filter according to [6.a] above. 

One way of doing this would be to formulate the magnitude spectrum H(z), 
inverse transform the spectrum to find the non-causal imulse tranform 
of the filter and then use the coefficients of the spectrum as 
coefficients of the filter difference equation and solve for the roots 
of that equation. Accepting only the roots z_m such that |z_m| < 1, 
we now should be able to find the coefficients of the "half filter"
G(z), which in turn could be cascaded by itself to find H(z). 

This could work if the number of coefficients in the FIR filter is 
"sufficiently small", whatever that may mean, as I think most of the 
work and most of the inaccuracy would be concerned with rooting the 
polynomial.

My 2 c.

Rune

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message
news:f56893ae.0311020443.7b0bf86b@posting.google.com...
> allnor@tele.ntnu.no (Rune Allnor) wrote in message
news:<f56893ae.0310310733.7b32304f@posting.google.com>...
> > "Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in message
news:<P9kob.4987$Nz5.422060@news20.bellglobal.com>...
> >
,,,,,,,,,,,,,,,,,,,,,,,,,,,,

> One way of doing this would be to formulate the magnitude spectrum H(z),
> inverse transform the spectrum to find the non-causal imulse tranform
> of the filter and then use the coefficients of the spectrum as
> coefficients of the filter difference equation and solve for the roots
> of that equation. Accepting only the roots z_m such that |z_m| < 1,
> we now should be able to find the coefficients of the "half filter"
> G(z), which in turn could be cascaded by itself to find H(z).
>
> This could work if the number of coefficients in the FIR filter is
> "sufficiently small", whatever that may mean, as I think most of the
> work and most of the inaccuracy would be concerned with rooting the
> polynomial.


Rune,

There are a multiplicity of filters with the same magnitude response -
depending on which zeros are inside and which zeros are outside the unit
circle.  Any pair of zeros with the same real part and imaginary parts of
opposite sign can be replaced by a pair of zeros (switching to polar
coordinates) with the same value for rho and with distance from the orgin of
1/L where L is the original distance.

To get minimum phase from the same magnitude response, all you have to do is
replace the zeros that are outside the unit circle with zeros (as above)
that are inside the unit circle ... that are conjugate to the ones being
replaced.

Fred