Matlab listing: mps.m and test program
This section lists and describes the Matlab/Octave function mps which estimates a minimum-phase spectrum given only the spectral magnitude. A test program for mps together with a listing of its output are given in §J.11 below. See §11.7 for related discussion.
function [sm] = mps(s) % [sm] = mps(s) % create minimum-phase spectrum sm from complex spectrum s sm = exp( fft( fold( ifft( log( clipdb(s,-100) )))));The clipdb and fold utilities are listed and described in §J.10 and §J.9, respectively.
The mps function works well as long as the desired frequency response is smooth. If there are any zeros on the frequency axis (``notches''), the corresponding minimum-phase impulse response will be time aliased because the corresponding exponentials in the cepstrum never decay. To suppress time-aliasing to some extent, the desired frequency response magnitude is clipped to 100 dB below its maximum. Time aliasing can be reduced by interpolating the desired frequency response s to a higher sampling density (thereby increasing the time available for exponential decay in the cepstral domain). However, for pure notches (zeros right on the unit circle), no amount of oversampling will eliminate the time aliasing completely. To avoid time aliasing in the cepstrum, such a desired spectrum must be smoothed before taking the log and inverse FFT. Zero-phase smoothing of the spectral magnitude is a typical choice for this purpose. When greater accuracy is required, all notch frequencies can be estimated so that terms of the form can be effectively ``divided out'' of the desired spectrum and carried along as separate factors.
Below is the test script (for ease of copy/paste extraction for online viewers). Following is the test script along with its output.
Note that this is a toy example intended only for checking the code on a simple example, and to illustrate how a spectrum gets passed in. In particular, the would-be minimum-phase impulse response computed by this example is clearly not even causal (thus indicating an insufficient spectral sampling density [FFT length]). Moreover, the starting spectrum (a sampled ideal lowpass-filter frequency response) has zeros on the unit circle, so there is no sampling density that is truly sufficient in the frequency domain (no sufficiently long FFT in the time domain).
spec = [1 1 1 0 0 0 1 1]'; % Lowpass cutting off at fs*3/8 format short; mps(spec) abs(mps(spec)) ifft(spec) ifft(mps(spec))
Below is the output of the test script with echo on in Matlab.
spec = [1 1 1 0 0 0 1 1]'; % Lowpass cutting off at fs*3/8 format short; mps(spec) ans = 1.0000 - 0.0000i 0.3696 - 0.9292i -0.2830 - 0.9591i 0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -0.2830 + 0.9591i 0.3696 + 0.9292i abs(mps(spec)) ans = 1.0000 1.0000 1.0000 0.0000 0.0000 0.0000 1.0000 1.0000 ifft(spec) ans = 0.6250 0.3018 - 0.0000i -0.1250 -0.0518 - 0.0000i 0.1250 -0.0518 + 0.0000i -0.1250 0.3018 + 0.0000i ifft(mps(spec)) ans = 0.1467 - 0.0000i 0.5944 - 0.0000i 0.4280 - 0.0000i -0.0159 - 0.0000i -0.0381 - 0.0000i 0.1352 - 0.0000i -0.0366 - 0.0000i -0.2137 - 0.0000i(The last four terms above would be zero if the result were causal.)
Signal Plots: swanalplot.m
Matlab listing: clipdb.m