1

I learned in school that any real-valued, even-symmetric sequence can be expressed as the convolution of a minimum-phase sequence and its time-reverse, which is maximum-phase. In my case, the real-valued, even-symmetric sequence is $sinc(x)$. What causal, real-valued, minimum-phase sequence convolved with its time-reverse produces $sinc(x)$?

I am familiar with the traditional approach of computing the frequency-domain zeros and constructing the minimum-phase sequence from those inside the unit circle and half of the pairs on the unit circle. However, in this case, the sequence has infinite length. I would prefer a closed-form expression for the minimum-phase sequence function.

In practice, I will be using a windowed version of the sequence -- an FIR filter. I would also like to learn how to deconvolve any real-valued, even-symmetric, finite sequence (FIR filter) into its minimum-phase component, whose frequency response magnitude is the square-root of the symmetric FIR. I would like a time-domain method for doing this because there are precision problems with computing frequency domain zeros in very long sequences, such as a million taps.

lennon310
  • 3,590
  • 19
  • 24
  • 27
hiccup
  • 13
  • 2

2 Answers2

0

For approximating a minimum phase sequence from a finite length linear phase sequence (such as a windowed Sinc), here are some existing stackoverflow and dsp.stackexchange answers regarding a matlab minimum phase approximation script which uses cepstrum deconvolution:

https://stackoverflow.com/questions/16753508/cepstrum-deconvolution-matlab-code

Derive minimum phase from magnitude

And here's an earlier post on the same matlab script, with more comments, from the music-dsp mailing list:

http://www.music.columbia.edu/pipermail/music-dsp/2004-february/059372.html

And the code:

x = [x; zeros(size(x))];
m = size(x,2);
n = size(x,1);
wn = [ones(1,m); 2*ones(n/2-1,m); ones(1,m);  zeros(n/2-1,m)];
y = real(ifft(exp(fft(wn.*real(ifft(log(abs(fft(x)))))))));
y = y(1:n/2;);

I've found that a greater amount of zero=padding of the original linear phase sequence may help improve the approximation.

Community
  • 1
hotpaw2
  • 35,346
  • 9
  • 47
  • 90
  • good to show cepstrum method for others who reference this topic. I, however, need good precision when deconvolving large sequences where FFTs are impossible or add too much noise. That is why I requested a time-domain approach. – hiccup Jan 30 '15 at 20:43
  • First, you might have to ask if there can even exist a (finite) time-domain approach which can produce the appropriate Hilbert transform pair. – hotpaw2 Jan 30 '15 at 21:17
  • In cases where I can use shorter impulse responses, the code and references you provided are working quite well. Still hoping for a closed-form solution to help with longer sequences. – hiccup Feb 04 '15 at 21:41
0

It is not true that "any real-valued, even-symmetric sequence can be expressed as the convolution of a minimum-phase sequence and its time-reverse, which is maximum-phase". You can see this as follows: if the symmetric sequence is $r[n]$, and the minimum-phase and maximum-phase components are $s[n]$ and $s[-n]$, respectively, then you have:

$$r[n]=s[n]*s[-n]$$

where $*$ denotes convolution. Assuming that $s[n]$ is real-valued, the spectrum of $r[n]$ is given by

$$R(\omega)=S(\omega)S^*(\omega)=|S(\omega)|^2\ge 0$$

i.e. $R(\omega)$ must be non-negative, but this is generally not the case for all even symmetric sequences. Their spectrum is just real-valued and symmetric, but not necessarily non-negative. So you have to require the additional condition of $r[n]$ having a non-negative spectrum, otherwise it can't be factored.

Now for the algorithms: apart from methods based on the Fourier transform and the Cepstrum, there are iterative methods such as the Wilson-Burg algorithm, and related methods, which are all based on Newton's method for computing a square root. This paper gives a nice overview and has further references.

Matt L.
  • 89,963
  • 9
  • 79
  • 179
  • Your point is well taken about the requirement for a non-negative spectrum. I couldn't figure out how to apply an iterative method, although the topic is quite interesting. Since I am trying to design filters that are hundreds of thousands of taps, an iterative method sounds really slow. Thanks for the help! – hiccup Feb 04 '15 at 21:49