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

###
Matlab listing: `mps.m`

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.

Note that `mps.m` must be given a *whole spectrum* in
``FFT buffer format''. That is, it must contain dc and positive-frequency
values followed by negative frequency values and be a power of 2 in
length.

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.

###
Matlab listing: `tmps.m`

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))

###
Matlab diary: `tmps.d`

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

**Next Section:**

Signal Plots: swanalplot.m

**Previous Section:**

Matlab listing: clipdb.m