Why construct a minimum phase filter from measurements?

Question

Sorry if this is a trivial question, but I am not doing signal processing everyday. I will try to express what I think I have understood as best as possible.

Suppose I apply a (linear) sweep signal to a system (like a loudspeaker in a room) and measure the response. To the input and output of that measurement I apply a finite Fourier transform each. By dividing the FT of the output by the FT of the input, I obtain an estimate of the transfer function $A(f)$ (within the limits of the finite measurement interval). The inverse Fourier transform of this "transfer function" gives me the impulse response of the system, and hence, the coefficients of the FIR filter that represents the system.

If I calculate the complex reciprocal of the transfer function in frequency space (possibly with some regularization to deal with zeroes), i.e. $B(f)=1/A(f)$, I get a transfer function, the inverse Fourier transform of which gives me an FIR filter that reverses the system's action on the signal as much as possible (i.e. it equalizes the room for example). I might then put the coefficients into a DSP board to carry it over into reality, i.e. a real room equalizer.

Now I have seen a YouTube video of someone working with the tool REW for room equalization, who claims (around timestamp 20:10) that one has to turn the reciprocated transfer function $B(f)$ of the system into a "minimum phase" version first, before using it as an "equalizer" of the system, i.e. before transforming it back to an inpulse response/FIR coefficients. Otherwise, so he claims, the filter's response would result in "massive ringing".

What is the theory behind such a claim and how is the construction of minimum phase actually done? Isn't it the case that the system response is trivially causal (because it is a physical system) and stable, and so its inverse must also be causal stable? So, shouldn't it be minimum phase automatically?

And wouldn't it be more natural to use the impulse response ($FT^{-1}(A)$) directly as the coefficients of an IIR filter, which would possibly provide a much better equalization with the same number of coefficients, simply because IIR filters reproduce damped oscillating systems (like sound in air and speaker membranes) better? In order to simulate an exponentially damped sine function, one would need a huge number of FIR coefficients if damping is slow.

Update (2023.02.25): I found that my preferred DSP (Analog Devices ADAU1701) has an "AutoEQ" block (in Sigma Studio), where I can load an impulse response file and let the tool equalize it automatically. It works by fitting a user-defined number of simple filters to the spectrum (or a range thereof). I guess the simple filters are all minimum-phase by design, so the result is also. One can see the optimizer at work in the UI graph window. This is exactly what I meant in my comment to Dan Boschen's answer: a more direct and efficient way to find the equalizer (as opposed to polynomial factorization).

Of course it is already reduced in filter order, and so it is not the same as minimizing phase over the full order of the acquired spectrum. But actually this is not a drawback but a feature, because it already realizes the necessary smoothing mentioned by Hilmar. Plus, the individual filters are IIR and hence, also work well in the low-frequency range. Anyway, it fits my purpose very well. And it is good to have a rough understanding as to how it is working under the hood.

Answer 1

The reason IIR filters are a poor choice is because for all cases except a minimum phase system, the IIR filter will be unstable. Unless we can be assured that the channel is minimum phase, the IIR solution directly as the inverse of the channel’s transfer function will not work. This point is explained in more detail in this related post.

The bottom line reason for this is as follows: A discrete time minimum phase system is characterized as having all poles and zeros inside the unit circle on the z plane. For a causal system to be stable, all poles must be inside the unit circle. (Similarly a continuous time minimum phase system has all poles and zeros in the left half plane on the s plane and for the causal system to be stable all poles need to be in the left half plane). When you invert a system, the poles become zeros and the zeros become poles, so only a causal minimum phase system has a stable causal inverse.

The other choices are linear phase, maximum phase or mixed phase- and for all those cases, which are causal and stable, they will have some zeros outside the unit circle. Thus if they were to be inverted, the zeros become poles and would be therefore unstable.

A very simple example is to compare the following impulse response given as coefficients: [1, 0.5], and [0.5, 1]. The associated z transforms for each of these two cases (from which we can determine the poles and zeros) are $H_1(z)= 1+0.5z^{-1}$ and $H_2(z)=0.5+z^{-1}$. The first case $H_1(z)$ is minimum phase, specifically proven by confirming all of its poles and zeroes are inside the unit circle:

$$H_1(z) = 1+0.5z^{-1} = \frac{z+0.5}{z}$$

Pole at $z=0$ and zero at $z=-.5$ Stable! 😊

Since all poles are inside the unit circle, it is a stable causal system. It's inverse also has all of its poles and zeros inside the unit circle:

$$H_1^{-1}(z) = \frac{z}{z+0.5}$$

Pole at $z=-0.5$, zero at $z=0$. Stable! 😊

The second case as the reverse filter (coefficients are reversed) will have the exact same magnitude response but in this case is a maximum phase filter (all zeros outside the unit circle) so will not have a stable causal inverse.

$$H_2(z) = 0.5+z^{-1} = \frac{0.5z+1}{z}$$

pole at $z=0$ and zero at $z=2$ Stable! 😊

$$H_2^{-1}(z) = \frac{z}{0.5z+1}$$

pole at $z=2$ and zero at $z=0$ NOT STABLE! 🙁

For this reason, if we wanted to pursue a true inverse filter to equalize a channel, we must remove the maximum possible delay (while still being causal) from the estimated channel response first. If the poles and zeros of the channel response are known, then creating the minimum phase response is done by simply reflecting the zeros outside the unit circle to inside the unit circle ($z = 1/z_p^*$). Doing this on the last example above would mean taking the pole that is at $z=2$ and moving it (reflecting it) to $z=0.5$, which brings us back to the minimum phase example above that. Note that both of these systems $H_1(z)$ and $H_2(z)$ have the exact same magnitude response it is just the phase that is different!

Why this works is made clearer by first understanding how any channel response that isn't minimum phase (meaning mixed phase or maximum phase) can be decomposed into a minimum-phase system cascaded with an all-pass system. The all-pass system has a flat magnitude response over all frequencies, so only modifies the phase, thus with that we pull out the maximum possible delay from the system leaving only the minimum phase system with the given magnitude response (for any given magnitude response, there is only one minimum phase solution). An all-pass system is realized by having all poles inside the unit circle, with a complex conjugate reciprocal zero associated with each pole, such as a demonstrate below showing the poles and zeros on the z-plane for a 2nd order all-pass with a complex pole at $z=z_1$:

The decomposition for an example maximum phase system shows how reflecting any zeros outside the unit circle to inside the circle at their reciprocal locations is the minimum phase system, cascaded with an all-pass. The all pass is done such that the new zeros in the minimum phase system are cancelled by the poles in the all-pass, leaving just the zeros outside the unit circle:

The minimum phase system will have the exact same magnitude response as the maximum phase system, but unlike the maximum phase system it has the desirable property of being invertible (has a stable causal inverse). The above demonstration was done with a maximum phase system, but also works with a mixed phase system. A mixed phase system has some zeros inside and some zeros outside the unit circle, so in this case the minimum phase system extracted is created by reflecting just the zeros that are outside the unit circle to be inside the unit circle at their conjugate reciprocal locations (which means same angle and $1/r$ where $r$ is the radius). Alternatively the Hilbert Transform can be used to determine the minimum phase system from the magnitude response alone, since for a minimum phase system, the phase is the Hilbert transform of the natural log of the magnitude response.

Both of these approaches are detailed further in these links:

Minimum Phase - All Pass Decomposition For Large Linear Phase Filters

Get minimum phase from function

However, an alternate approach that doesn’t require reducing the channel to minimum phase is using a least squared equalizer as detailed in this post: Compensating Loudspeaker frequency response in an audio signal

Creating filters using the "frequency sampling" approach, meaning deriving the impulse response from the inverse FFT of a desired frequency response as suggested by the OP is not recommended over alternate approaches such as least squares outlined above, and firls in MATLAB, Octave and Python, and windowing approaches to filter design. This is because the inverse FFT suffers from time domain aliasing, requiring additional zero-padding which results is a much longer filter to achieve the same performance that the other approaches listed would provide. This is demonstrated in this link. The frequency sampling approach to filter design results in frequency response that has an exact match at the frequency samples chosen, but greater error everywhere else compared to the other approaches given the same complexity (number of coefficients). To mention however, the case where "frequency sampling" design excels is OFDM (orthogonal frequency division multiplexing) since in that case we are only interested in discrete frequencies.

Answer 2

Room Equalization for listening purposes is NOT the same as equalization for, say, acoustic echo cancellation or data transmission purposes.

Inverting the measured impulse response or transfer function and applying it to the playback signal would sound absolutely terrible (for a variety of reasons, ask a different question if interested).

The standard process for listening equalization is therefore quite different

1. Measure transfer function, ideally at multiple representative locations
2. Apply logarithmic spectral smoothing that's mimics roughly how human spectral perception works. A good starting point is 3rd octave smoothing.
3. Energy average over all measurements
4. Invert over the frequency range of interest
5. Apply minimum phase
6. Design FIR or IIR filter through least square or direct inverse Fourier Transform

Energy averaging and smoothing destroys any phase information (on purpose) so the resulting spectral target has no phase information. The best choice of phase here is indeed a minimum phase. A linear phase would result in excessive pre-ringing and increase the number of filter coefficients and wouldn't work for IIR filters.

IIR filters are indeed the most popular choice for a room EQ filter. You can get a decent fit with a very small number of biquads (< 20) while an FIR would require 1000s if not 10s of thousands of coefficients.