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?

# Minimum Phase Impulse Response

Started by ●October 29, 2003

Reply by ●October 30, 20032003-10-30

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-timeimpulse> 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]

Reply by ●October 30, 20032003-10-30

"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-timeimpulse> 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

Reply by ●October 30, 20032003-10-30

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 thatdoesn'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 youhave> 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.

Reply by ●October 31, 20032003-10-31

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 havezeros> 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 todevelop> > 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 todetermine> 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

Reply by ●October 31, 20032003-10-31

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

Reply by ●October 31, 20032003-10-31

"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

Reply by ●October 31, 20032003-10-31

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

Reply by ●November 2, 20032003-11-02

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

Reply by ●November 2, 20032003-11-02

"Rune Allnor" <allnor@tele.ntnu.no> wrote in message news:f56893ae.0311020443.7b0bf86b@posting.google.com...> allnor@tele.ntnu.no (Rune Allnor) wrote in messagenews:<f56893ae.0310310733.7b32304f@posting.google.com>...> > "Matt Timmermans" <mt0000@sympatico.nospam-remove.ca> wrote in messagenews:<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