Phase and group delay
Contents
10.4. Phase and group delay#
In the previous section, we saw how to investigate the properties of a convolutional filter \(\red{h}\) by examining its frequency response, that is, the magnitude of its DFT \(|\red{H}|\). Examining the magnitude tells us an important part of the story, since it directly encapsulates the gain applied to each frequency when the filter is applied to a signal. However, it’s not the entire story: phase is also important.
Recall that the convolution theorem converts time-domain convolution into frequency-domain multiplication:
Since each \(\red{H[m]}\) is a complex number, we can express it as \(\red{H[m]} = A \cdot e^{\mathrm{j}\phi}\) for magnitude \(A\) and phase \(\phi\). This tells us that the output DFT component \(\magenta{Y[m]}\) will be derived from \(\darkblue{X[m]}\) by scaling (multiply by \(A\)) and rotation by angle \(\phi\):
From the DFT shifting theorem, we know that applying delay to a signal results in a rotation of each DFT component \(\darkblue{X[m]}\), so we can interpret the rotation by \(\phi\) as implementing a delay of the \(m\)th sinusoid. The shifting theorem states precisely how this works: a delay of \(\red{k}\) samples induces a rotation of \(-2\pi\cdot \red{k \cdot m} / \blue{N}\) for frequency index \(\red{m}\) with signal length \(\blue{N}\), so that each frequency is rotated by an angle that depends on both the delay \(\red{k}\) and the frequency index \(\red{m}\).
However, there’s no guarantee that an arbitrary filter \(\red{H}\) will adhere to this structure. Depending on the phase structure of \(\red{H}\), the resulting signal could sound completely distorted. Analyzing the phase structure of \(\red{H}\) can reveal the presence of these distortions.
10.4.1. Examples#
Before going into analyzing filter phase, it will help to have some intuitive grasp of how important phase is to a signal.
The examples below demonstrate filters which do not affect the magnitude spectrum at all: \(|\magenta{Y[m]}| = |\darkblue{X[m]}|\), but the output signal \(\purple{y} = \red{h}*\blue{x}\) is substantially different from the input signal \(\blue{x}\).
Example 10.1 (Discarding phase)
To demonstrate the importance of phase, we can listen to what happens when we discard all phase information in a signal, so that
This is equivalent to applying a filter with unit amplitude and phase exactly opposite of \(\darkblue{X[m]}\):
This can be implemented in Python as follows:
# Take the DFT
X = np.fft.rfft(x)
# Discard phase
Y = np.abs(X)
# Invert the DFT
y = np.fft.irfft(Y)
and the results are demonstrated below.
'Original signal: https://freesound.org/s/50712'
'After removing phase'
Discarding phase information preserves some aspects of the signal: the same frequencies are generally present, but the temporal coherence of the signal has been destroyed.
Example 10.2 (Random phase)
The example above is particularly bad, but not unique. As a second example, we can imagine creating a filter \(\red{H}\) with random phases between 0 and \(2\pi\).
# get the length and DFT of a signal x
N = len(x)
X = np.fft.rfft(x)
# random(N) makes N random numbers between 0 and 1
# multiplying by 2pi gives us random phases
phase = np.random.random(len(X)) * 2 * np.pi
H = np.exp(1j * phase)
# Rotate each X[m] by a random amount
Y = X * H
# recover the time-domain signal
y = np.fft.irfft(Y)
The result of this transformation is demonstrated below. In this case, randomizing the phase has rendered the original signal completely unintelligible.
Phase is important.
'Random phase'
10.4.2. The phase spectrum#
Let \(\red{H}\) denote the DFT spectrum of an impulse response \(\red{h}\). Rather than plot the magnitude spectrum \(|\red{H}|\), the phase spectrum is the sequence of phases \(\phi[m]\), where
In Python, the phase spectrum can be obtained as follows:
# Take the DFT of the filter
H = np.fft.fft(h)
# Extract the angle: between -pi and +pi
phase = np.angle(H)
Plotting \(\phi[m]\) as a sequence can help us understand the delay behavior of the filter for each analysis frequency \(m\). However, correctly interpreting the phase spectrum requires a bit more care than interpreting the magnitude spectrum.
If we interpret delay as a convolutional filter, we can plot the phase spectrum of different delays to see what happens. Some examples are displayed below in Fig. 10.9.
Fig. 10.9 Left: three different delay filters (\(\blue{1}\), \(\red{4}\), and \(\green{7}\) samples) are plotted with length \(N=64\). Right: for each filter, the phase spectrum \(\phi[m]\) is plotted as a function of frequency.#
From Fig. 10.9, we can observe that the different delay values make zig-zag shapes with different slopes. This is exactly what is predicted by the DFT shifting theorem, which implies
for a \(\red{k}\)-sample delay. When plotted as a function of \(\red{m}\), \(\phi[m]\) should produce a line with slope \(-2\pi \cdot \red{k} / N\), but this is not exactly what we see.
The zig-zag effect is due to angles wrapping around at \(\pm\pi\).
The np.angle function extracts the angle from a complex number, but there are ambiguities in how angle is interpreted.
Remember that \(\theta \equiv \theta + 2\pi\) for any angle \(\theta\), so np.angle will always take the equivalent angle that lies within the range \([-\pi, +\pi)\).
This effect can be undone by using the np.unwrap function.
Angle unwrapping takes a sequence of angles \(\phi[m]\) and produces a sequence of equivalent angles \(\phi'[m]\) that avoid large discontinuities (greater than \(\pi\) radians) between \(\phi[m]\) and \(\phi[m+1]\) by adding multiples of \(2\pi\).
The unwrapped angles will land outside the \(\pm\pi\) range, but satisfy \(\phi'[m] = \phi[m] + c\cdot 2\pi\) (for some integer \(c\)).
In code, this is done as follows:
# Extract the angle: between -pi and +pi
phase = np.angle(H)
# Unwrap the phase
phase_unwrap = np.unwrap(phase)
The results of phase unwrapping on the delay filters in Fig. 10.9 are illustrated below.
Fig. 10.10 Top: the phase spectra \(\phi[m]\) of three different delay filters (\(\blue{1}\), \(\red{4}\), and \(\green{7}\) samples) appear as zig-zag lines, bounded between \(-\pi\) and \(+\pi\). Bottom: the unwrapped phase spectra \(\phi'[m]\) appear as straight lines with slope depending on the amount of delay.#
10.4.3. Group delay#
The delay examples above illustrate a general principle. If a filter’s (unwrapped) phase spectrum is linear — i.e., \(\phi'[m] = a\cdot m\) for some constant \(a\) — then the filter implements a delay that preserves the relative positioning of each sinusoidal component. In other words, filters with linear phase preserve the phase coherence of any input signal.
Conversely, a filter with non-linear phase, like the examples above, will delay different sinusoidal components by different amounts, and can significantly distort the signal. For this reason, linear phase is generally a desirable property.
For delay filters, the delay parameter can be recovered from the phase spectrum by rearranging (10.2) to solve for \(k\):
and this value would be the same for all frequency indices \(m>0\). The sample delay \(\red{k}\) can be converted into delay time \(t\) (measured in seconds) by dividing by the sampling rate:
Note that the delay \(t\) is proportional to the negative slope of the phase spectrum.
More generally, if someone hands you an arbitrary filter \(\red{h}\), you wouldn’t know in advance whether it has linear phase or not, and it would be incorrect to apply (10.3) because it is built on the assumption of linear phase. However, we can generalize the idea to relax this assumption.
Instead of assuming linear phase — that \(\phi[m] \propto \red{m}\), or equivalently, that the slope of the phase response is constant — we can estimate the negative slope separately for each frequency \(\red{m}>0\). If the filter truly has linear phase, then all estimates should agree, and give us the desired delay parameter. If the filter does not have linear phase, we can detect this by seeing disagreements in the slope estimates.
The simplest way to estimate slope is the “rise-over-run” method, that is, dividing the change in neighboring (unwrapped) phase measurements (the rise)
by the change in the corresponding frequencies (the run, in Hz):
The result is a sequence of negative slope measurements:
Applying dimensional analysis to (10.5), the numerator \(\Delta_\phi[m]\) is a difference of angles measured in [radians], while the denominator \(\Delta_f[m]\) is a difference of frequencies measured in [cycles / second].
The ratio therefore has units [radian-second / cycles], which is not entirely easy to understand.
However, if we divide by \(2\pi\) [radians / cycle], the result will have units of [seconds], which nicely corresponds to our intuitive notion of delay.
Putting this all together gives us the formal definition of group delay.
Definition 10.2 (Group delay)
Let \(\red{h}\) be an impulse response of length \(\blue{N}\) samples, and let \(\phi'[m]\) denote its unwrapped phase spectrum.
The group delay of \(\red{h}\) is the given by the sequence:
and is measured in [seconds].
For FIR filters, group delay can be calculated in Python as follows:
# Take the DFT
H = np.fft.fft(h)
# Extract angles
phase = np.angle(H)
# Unwrap the angles
phase_unwrap = np.unwrap(phase)
# Calculate group delay from the differences in
# successive phase measurements
# We tack the last phase onto to the beginning so that the m=0 case
# is handled correctly, and compares to m=-1.
G = - len(h) / (2 * np.pi * fs) * np.diff(phase_unwrap, prepend=phase_unwrap[-1])
Alternatively, the scipy module provides a function for computing group delay for both FIR and infinite impulse response (IIR) filters, which we will see in the next chapter:
# First get our frequencies
frequencies = np.fft.fftfreq(len(h), 1/fs)
# Compute group delay (in samples) using scipy
# - (h, 1) is how we indicate that this is an FIR filter
# - w=frequencies is how we specify which frequencies to measure (DFT analysis frequencies)
frequencies, G_samples = scipy.signal.group_delay((h, 1), w=frequencies, fs=fs)
# Convert delay measurements from samples to seconds
G = G_samples / fs
10.4.4. Linear and non-linear phase#
As mentioned above, not all filters will have completely linear phase. To demonstrate this, let’s create two filters, one low-pass and one high-pass, using the window method:
fs = 22050
f_cutoff = 3000
# 4 cycles of our cutoff frequency
order = 4 * fs // f_cutoff
h_low = scipy.signal.firwin(order, f_cutoff, window='hann', pass_zero='lowpass', fs=fs)
h_high = scipy.signal.firwin(order, f_cutoff, window='hann', pass_zero='highpass', fs=fs)
Fig. 10.11 illustrates the impulse response, frequency response, and unwrapped phase spectrum of these two filters.
Fig. 10.11 Two FIR filters are shown (left column), along with their frequency response \(|H[m]|\) (center column) and unwrapped phase spectrum \(\phi'[m]\) (right column). Shaded regions correspond to the stop-bands of the filters in the frequency domain. Top row: a low-pass filter with \(f_c=3000\) Hz. Bottom row: a high-pass filter with \(f_c=3000\) Hz. Neither filter has completely linear phase response, but both are linear when restricted to the pass band, as indicated by the solid line.#
Unlike the delay filters depicted in Fig. 10.10, these filters do not have a strictly linear response. In both cases, the unwrapped phase \(\phi'[m]\) has a sharp bend near the cutoff frequency \(f_c\), where a strictly linear phase response would continue in the same direction indicated in the pass-band.
However, if we restrict attention to the pass band of the filters, they both have linear phase response, as indicated by the solid lines in Fig. 10.11. The non-linear part of the phase response is limited to the stop band, and this is fine. Frequencies in the stop band should be attenuated (and ideally, inaudible), and even if phase distortion happens there, it won’t matter because we won’t hear it.
For this reason, it is typically sufficient to have linear phase only in the pass band of a filter, and it is common to refer to such filters as having linear phase, even if their entire phase response is technically non-linear.
10.4.5. Summary#
In this section, we’ve seen how to analyze the phase spectrum of FIR filters. This lead to the definition of group delay, which allows us to determine A) if a filter is non-linear, and B) how much delay a filter will apply to a signal.
Practically speaking, the filters we’ve covered in this chapter (window method and Parks-McClellan) are designed to have linear phase in their pass-bands, so subsequent analysis is not strictly necessary to determine linearity. That said, the method demonstrated above can still be useful if you need to infer the delay of a given filter.
In the next chapter, we will introduce infinite impulse response filters, which generally will not have linear phase. Understanding group delay in that setting will be much more important, so it’s a good idea to get some practice with the method in the FIR setting first.