Return to Table of Contents

Appendix A. Inverting a stable filter that is not minimum phase

Statement of problem

Oftentimes a signal processor is presented with a situation that is modeled (appropriately or not) by an unstable digital filter. The challenge is then to take that unstable filter and modify it to get a stable filter without disrupting either the magnitude or phase response or both. This can be difficult if not impossible to do and retain causality, so many times the resulting filter is noncausal but a stable and excellent approximation to the unstable filter. Noncausality is a small price to pay if magnitude and phase information are critical to performance.

Example of a phase-critical situation

In my case it was necessary to build a model of the GEMS's filter characteristics, as the lower 3 dB frequency was low enough (about 70 Hz) to result in noticeable distortion of signal components below about 300 Hz. As a significant amount of energy is contained in the GEMS signal below this frequency (especially for males), it is important to know what the "raw" (undistorted by filters) GEMS signal was so that we could correlate it to the physical position of reflecting interfaces. I therefore used a sweeping frequency generator to inject artificial signals into the GEMS’s receiver and from that generated a magnitude and phase versus frequency plot. From that, I built a digital model of the GEMS’s filters with the System Identification Toolbox in Matlab using a method described in section 2.1.2. Figure 2.8 shows the resulting plots. Notice that after about 300 Hz, there is not much further change in phase and the magnitude is pretty flat after about 300 Hz. What we want to do, then, is build a filter that is the perfect inverse of the filters in the GEMS so we may remove their effects from the signal.

Unfortunately, the best fit using the output error model was perfectly stable but not minimum phase, which meant that upon inversion there was a single pole outside the unit circle, denoting instability. Since the filter response has to be inverted so that the effects of the GEMS filters can be removed, we must stabilize this inverse filter. Correcting the phase is important since (for once) we do care what the signal looks like and are not just concerned with its power spectrum.

A brief primer on the z plane

An introduction to the concept of magnitude and phase response in the z plane is warranted to make the following discussion more easily understood.

The poles and zeros of a filter are by definition the roots of the denominator and numerator of the transfer function, respectively. Their values determine everything about the characteristics of the filter – its magnitude, phase response, and stability. The values of the poles and zeros may be plotted on either an x-y (real-imaginary) plot (common in control engineering) or a z-plane plot (more common in signal processing, where z is , r always positive). Both display the same information in a slightly different format. As I am primarily concerned with the z plane, I will explain the procedure using that as my basis.

As an example I will use a filter with a pole at z = 3/5 and a zero at z = 5/4. This is a rather simple case, but is representative of the type of problem and illustrates the points sufficiently. The case of multiple poles and zeros is easily handled by superposition. The frequency response H(z) is represented by:

A.1

This is the form of a highpass filter and is stable, as the only pole is inside the unit circle. The normalized frequency response is shown in figure A.1. Although it is stable, it is not minimum phase, which requires all poles and zeros be located inside the unit circle. When it is inverted, the pole exchanges positions with the zero and is now on the outside, causing instability. The inverse of the filter H(z) is

A.2

A representative z-plane plot of the inverted, unstable filter H^-1(z) is shown in figure A.2.

Figure A.1. Normalized frequency response for the example highpass filter.

In the z plane, the unit circle is special in that it represents the points where r = 1, and z reduces to the complex exponential . On these points, and the z transform is equal to the Fourier transform. It is from points in the z plane on the unit circle that we determine the frequency and phase response of the system.

At , we are on the unit circle and at the point corresponding to DC. As we make our way around the circle, we go through all possible frequencies until we make it to , which corresponds to the highest possible frequency attainable in a digital system, the Nyquist frequency. This is equal to ½ the sampling frequency, and the analog signal sampled by the digital system must have no components equal to or greater than this frequency or the resulting digital signal will be distorted. The lower half of the unit circle adds no new information; it is the mirror of the upper and contains the response for "negative" frequencies, a mathematical artifact that has little bearing on the real world. Somewhat like an old professor of my acquaintance.

Anyway, determining the magnitude response of the filter at a certain frequency is simply a matter of determining the distance from its position on the unit circle to all the poles and zeros of the system. The magnitude response at is given by

A.3

with n the number of zeros and m the number of poles.

Similarly, the phase response of the system is given by

A.4

Figure A.2 illustrates how these lengths and angles are calculated. From this example it is possible to get an idea how an allpass filter has a flat magnitude response. For an allpass, the zero is located at z = B (which may be complex), and the pole is located at z = 1/B. Thus the distance from any point on the unit circle to the pole/zero pair varies, but the ratio of the distances remains constant.

Now that we have a handle on the z-plane, we can begin the explanation of how we can repair unstable filters using one or more anticausal allpass filters.

Stabilizing unstable filters using allpass noncausal filters

The first step is to build an allpass (AP) filter to cancel out the offending exterior pole while preserving the magnitude response. An allpass filter has by definition unity magnitude response, but the phase is not zero, resulting in a 2-stage filter that now has a different phase response than desired. In Figure A.3, the inverted filter is shown with the pole outside the unit circle at z = 5/4, and the zero inside at z = 3/5.

This is an unstable filter, and so we build an AP filter with a zero at z = 5/4 and a pole at z = 4/5. The offending pole at 5/4 is cancelled out and a new, minimum phase filter results. As seen in figure A.4, the magnitude response of the filter is inverted, but the phase response is not and is quite distorted. The frequency response of this inverted, stable filter H_s(z) is

A.5

How badly is the phase distorted? Well, it is quite distorted near w = 0, as the pole is now contributing ~ 0 degrees to the phase where it was contributing ~ p when it was outside the unit circle. It soon becomes clear that in order to cancel out the phase of the original filter we need a filter with negative phase at w = 0. As you can see in Figure A.2, a negative phase at w = 0 can only exist if the pole lies outside the unit circle, rendering the filter unstable. Therefore, in order to get negative phase at w = 0, we must utilize a noncausal filter.

So now that we have decided on a noncausal filter, how do we go about designing it? It must be an allpass filter so as not to disturb the magnitude response, and the phase should be chosen to be the negative of what we need to fix the distortion. Thus, we design for a certain phase response which, when subtracted from the present phase, will give us something that is closer to the phase response we are trying to get. It may not be perfect, especially for complex filters, but arbitrarily many of these allpass filters may be used so that we can get as close as necessary to the phase we desire.

To begin, we look for a frequency where the phase is important and the distortion is large. In Figure A.4, we see that the maximum error occurs at 0 Hz, but as the lowest voiced frequency in males is normally around 100 Hz, we will try there first. At 100 Hz we determine the correction needed, in this case about 75 degrees or 1.31 radians. As our sampling frequency in this example is 1 kHz, p is equivalent to 500 Hz so that our correction will take place at . We need to determine the location of the pole and zero of the allpass filter so that the phase of the allpass is 75 degrees at 100 Hz. In Figure A.5 the new allpass (to be anticausal) filter is plotted along with the inverted, stabilized filter with poor phase response. To the right are three triangles we will use in our calculations. They are magnified and not to scale for clarity.

We have determined that we need 1.31 radians to correct the filter phase at f = 100 Hz (w = 0.2p). From the top triangle in figure A.5, we get the values for x and y:

A.5

Figure A.4. The inverted, stable filter H_s(z). The magnitude is perfectly inverted but the phase is not.

A.6

which we can then use to find f_p and f_z:

A.7

A.8

Now, since the phase contributed by the new zero and pole is measured from a horizontal line (parallel to the w = 0 line, see figure A.2), the total phase of the new AP filter will be

. A.9

Since we want this to be equal to 1.31, we have the following equation:

A.10

or

A.11

and since

A.12

we get

. A.13

Using the values for x and y in equations A.5 and A.6, we solve for B and find 1.664 to be the only positive answer. This leads us to an AP filter with the transfer function of

where the normalization factor (1/B) in the numerator is necessary to make the magnitude equal to one. It comes from the general AP filter form:

A.14

where in this case a is real. The phase response for the allpass filter (the magnitude is not shown as it is unity for all frequencies) is shown in Figure A.6. The phase response for 100 Hz is very nearly 75 degrees. It is not exactly 75 degrees due to rounding errors.

The new combined filter H_c(z) is thus

A.15

But we must remember that the addition of the causal allpass filter will make the phase of the combination inverse filter/allpass equal to the original filter at 100 Hz. We want the phase to be the inverse of the original, not equal to it. Therefore we must use the filter backwards in time (anticausally) in order to retain the magnitude response and invert the phase response.

Mathematically this is accomplished by replacing z^-1 (a delay) with z (an advance), rendering the filter anticausal.

Figure A.6. Phase response for causal allpass filter designed to have a phase shift of 75 degrees at 100 Hz.

To accomplish this in Matlab we just take the coefficients of H_c(z) form a filter using the filter algorithm. To make it anticausal, we simply take the signal we are going to filter and reverse it in time:

A.15

where N is the length of the signal. Now we filter the reversed signal with the causal filter to get the phase corrected output:

A.16

where h_c is the impulse response of H_c(z). The frequency response of the original filter H(z), the stable, inverted filter H_s(z), and the final noncausal filter H_c(z^-1) are all plotted in Figure A.7. In this simple example, an almost perfect inversion of the phase was possible.

Sadly, for more complex problems this will not be the case. As we are only compensating the phase at a single frequency, at others the phase may still not be a good enough cancellation. Therefore, it is sometimes necessary to try to fit the phase at several different frequencies near the one you chose to begin with to see if you can get a better fit.

If you have very exact specifications, you may try and add additional causal or anticausal AP to further correct the phase in a particular frequency range, but be warned that the fit near w = 0 will suffer. This is because placing poles and zeros on the real axis at w = 0 causes another p to be added to the phase at w = 0. Thus near w = 0 the phase correction gets steadily worse as we add AP filters. If, however, your frequencies of interest are not near w = 0, this technique can work quite well.

Conclusions

After you are finished with the design of your stable, noncausal filter, it is wise to test the result if the inversion is not perfect. This will give you a way to check the distortion still present at different frequencies. A quick way to do this is to build an artificial signal, such as a square wave, filter it using your original model, and then with your inverted and stabilized model. You should get back out the original signal if you have sufficiently compensated for the phase distortion. This is assuming, of course, that your instability resulted from the inversion of a stable filter. If you are trying to compensate for an unstable model, there is no way to compare waveforms directly – you just have to be satisfied that your phase and magnitude plots are accurate enough. We will not perform such a test, as our inversion was perfect and it would not be illustrative.

Figure A.7. The frequency response for the original filter H(z), the stable inverted filter H_s(z), and the noncausal filter H_c(z^-1).

Return to Table of Contents