InterludeCepstral Analysis
A variety of signal processing applications use the collection of nonlinear techniques known as cepstral analysis. These are based around the core concept of the complex cepstrum. Loosely defined, the complex cepstrum of a signal is defined in terms of its Z transform. , the Z transform of the cepstrum, is defined as the logarithm of , the Z transform of the sequence . Alternately, we can write the definition for the cepstrum directly.
One of the more important properties of the cepstrum is that it is a homomorphic transformation. A homomorphic system is one in which the output is a superposition of the input signals, i.e., the input signals are combined by an operation that has the algebraic characteristics of addition. Under a cepstral transformation, the convolution of two signals becomes equivalent to the sum of the cepstra of the signals .
These properties are described in full detail in Oppenheim and Schafer (1989). The example that follows was drawn from Section 12.8 of that work.
The complex cepstrum can be very difficult to compute analytically. If we put appropriate stability conditions on the signals involved, however, we can define the cepstra in terms of the discrete Fourier transform. Some aliasing compared to the regular cepstrum can be expected. This can be compensated for by wellchosen sampling ranges.
To begin the example, we first need to load the package stubs.
This makes routines from Signals and Systems available for use.
In[1]:=Needs["SignalProcessing`"]
Let us define a signal that will be convolved with a short sequence of impulses. This is equivalent to creating a signal with echoes.
Here is a secondorder transfer function that will define our input signal.
In[2]:=vsystem = (b0 + b1 z^(1))/ ((1  r Exp[I theta] z^(1)) * (1  r Exp[I theta] z^(1)));
Here is the timedomain representation of the input.
In[3]:=vsignal = InverseZTransform[vsystem, z, n]
Out[3]=
Let these be the parameter values we will use in our signal.
In[4]:=params = {b0 > 0.98, b1 > 1, r > 0.9, theta > Pi/6};
This is what the input signal looks like.
In[5]:=DiscreteSignalPlot[vsignal/.params, {n, 0, 80}]
Out[5]=
As an example of a cepstrum that can be computed analytically, we can look at the cepstrum of this signal. We can compute the cepstrum from the transfer function. It is useful to make the transfer function analytic on the unit circle. We can do this by multiplying the transfer function by ; this is equivalent to shifting the sequence by one unit, a feature we will make use of when working with the discrete Fourier transform.
Here is a complex cepstrum computed analytically.
In[6]:=InverseZTransform[ Log[z vsystem], z, n]
Out[6]=
To work with the cepstrum via the discrete Fourier transform, it is useful for us to manipulate the signal as a sequence of data. We can generate the sequence by sampling the function.
We will work with the first samples of the signal.
In[7]:=vdata = ToSampledData[ vsignal/.params, {n, 0, 1023} ];
The pulses with which we will convolve the signal should also be sampled.
This sequence of pulses will be convolved with the input signal. It might represent the impulse response of any system that generates echoes at a spacing of and , where here we let be .
In[8]:=kdata = ToSampledData[ DiscreteDelta[n] + 0.9 DiscreteDelta[n  15] + 0.81 DiscreteDelta[n  30], {n, 0, 30} ];
The signal we will work with is the convolution of the two signals previously defined.
In[9]:=data = DiscreteConvolution[vdata, kdata, n];
Here, we plot the convolution of the signals.
In[10]:=DiscreteSignalPlot[data, {n, 0, 80}]
Out[10]=
The actual computations are fairly straightforward. In a few cases, we will find it useful to have auxiliary numbers such as half the number of data points.
This value will be used in several computations.
In[11]:=hlen = Floor[Length[First[data]]/2];
In the computation of the complex cepstrum, we use the complex logarithm of the sequence, defined as . For uniqueness, it is necessary that the phase be "unwrapped", which eliminates the jumps as the phase passes between and . This careful definition causes the complex cepstrum of a real sequence to also be a real sequence. Following is a utility function that will perform the phase unwrapping for discrete numeric data. This technique may fail for extremely oscillatory functions, and is not appropriate for very noisy data.
This function attempts to unwrap the phase of discrete data, by adding integer increments of 2 Pi to the data whenever it jumps by more than Pi. (These values can also be set by the user.)
In[12]:=unwrapphase[data_?VectorQ, tol_:Pi, inc_:2 Pi] := data + inc FoldList[Plus, 0, Sign[ Chop[ Apply[Subtract, Partition[data, 2, 1], {1} ], tol ] ] ]
We also need to remove a linear component of the phase of the discrete Fourier transform so the phase can be continuous at and at (scaled values). We do this by rotating the data set, making the transform equivalent to that of the sequence instead of . This can be compensated for later if necessary.
For these computations, we will find it convenient to use the vector form of DiscreteFourierTransform. In the case of general sequences, the syntax that involves functions is better.
In[13]:=ftrans = DiscreteFourierTransform[ RotateLeft[First[data]] ];
Given the Fourier transform, we can find the real cepstrum (as opposed to the complex cepstrum) of the data sequence quite easily. While not as useful for deconvolution applications, the real cepstrum (usually simply referred to as the cepstrum) is applied where the energy in various parts of the signal needs to be computed.
This computes the real cepstrum of the data set. Note the application of Chop to remove any small complex components that may remain due to rounding effects.
In[14]:=rcep = InverseDiscreteFourierTransform[ Log[Abs[ftrans]] ]//Chop;
Here is the real cepstrum of the data set. Convention places in the center of the plot, so we rotate the data set and offset it to the left.
In[15]:=DiscreteSignalPlot[ SampledData[RotateLeft[rcep, hlen], hlen, n], {n, 100, 100}, Frame > True, PlotRange > All ]
Out[15]=
Continuing on to the computation of the complex cepstrum, we first examine the phase unwrapping of the utility function defined previously. As it happens, in this particular example, the phase does not need to be unwrapped. The linear offset introduced by rotating the data set did the work for us. In general, however, you will still need to unwrap the phase, so we will quickly examine the wrapped and unwrapped phase of the unshifted data.
First, generate the phases of the discrete Fourier transform of the unshifted data.
In[16]:=unshifted = Arg[DiscreteFourierTransform[First[data]]];
It is only necessary to display the first half of the data set due to the usual symmetry in the transform. We use ListPlot to display the data so that we can join the points into a continuous line.
In[17]:=ListPlot[ unshifted, PlotRange > {{1, hlen}, All}, PlotJoined > True ]
Out[17]=
Here is the unwrapped phase.
In[18]:=ListPlot[ unwrapphase[unshifted], PlotRange > {{1, hlen}, All}, PlotJoined > True ]
Out[18]=
We now have all we need to find the complex cepstrum.
The full complex cepstrum is computed with the complex logarithm.
In[19]:=ccep = InverseDiscreteFourierTransform[ Log[Abs[ftrans]] + I unwrapphase[Arg[ftrans]] ]//Chop;
Here is what the complex cepstrum looks like, using the same plotting convention as used with the cepstrum.
In[20]:=DiscreteSignalPlot[ SampledData[RotateLeft[ccep, hlen], hlen, n], {n, 100, 100}, Frame > True, PlotRange > All ]
Out[20]=
Given the complex cepstrum, we can use techniques similar to frequencydomain filtering methods to deconvolve the signal and its echoes. In this example, the "lowtime" portion of the cepstrum corresponds to the input signal, so by windowing the signal with an appropriate pulse, we can eliminate the "hightime" portion due to the convolved impulses. (See, for example, Section 12.8.4 of Oppenheim and Schafer for more details.) This technique is called homomorphic deconvolution.
Here is a vector of data that may be used to filter the cepstrum. Recall that the original data vector wraps the negative times to vector indices greater than , so it is convenient to use two pulses to filter the times near zero.
In[21]:=tfilt = Table[DiscretePulse[14, t] + DiscretePulse[14, t  Length[ccep] + 14], {t, 0, Length[ccep]  1}];
Computation of the inverse cepstrum is somewhat simpler than the cepstrum, since special care is not required with respect to the phase.
The inverse cepstrum is computed by taking the inverse transform of the exponent of the Fourier transform. We use Chop to eliminate any small complex components due to rounding, and RotateRight compensates for the shifting of the data that made the phase continuous.
In[22]:=filtsig = InverseDiscreteFourierTransform[ Exp[DiscreteFourierTransform[ ccep tfilt ]] ]//Chop//RotateRight;
Here is a plot of the resulting data. Notice how close it is to our original input signal.
In[23]:=DiscreteSignalPlot[ SampledData[filtsig, n], {n, 0, 80} ]
Out[23]=
The complex cepstrum has found broad application in speech processing, seismic analysis, and many other fields. Signals and Systems and Mathematica can be used to investigate this and other transforms with great ease.
