Chapter 3: Time-discretization

Now that we have introduced the basic ideas of analog filter analysis, we will develop an approach to convert analog filter models to the discrete time.

3.1 Discrete-time signals

The discussion of the basic concepts of discrete-time signal representation and processing is outside the scope of this book. We are assuming that the reader is familiar with the basic concepts of discrete-time signal processing, such as sampling, sampling rate, sampling period, Nyquist frequency, analog-to-digital and digital-to-analog signal conversion. However we are going to make some remarks in this respect.

As many other texts do, we will use the square bracket notation to denote discrete-time signals and round parentheses notation to denote continuous-time signals: e.g. $x[n]$ and $x(t)$.

We will often assume a unit sampling rate $f_s=1$ (and, respectively, a unit sampling period $T=1$), which puts the Nyquist frequency at $1/2$, or, in the circular frequency terms, at $\pi$. Apparently, this can be achieved simply by a corresponding choice of time units.

Theoretical DSP texts typically state that discrete-time signals have periodic frequency spectra. This might be convenient for certain aspects of theoretical analysis such as analog-to-digital and digital-to-analog signal conversion, but it’s highly unintuitive otherwise. It would be more intuitive, whenever talking of a discrete-time signal, to imagine an ideal DAC connected to this signal, and think that the discrete-time signal represents the respective continuous-time signal produced by such DAC. Especially, since by sampling this continuous-time signal we obtain the original discrete-time signal again. So the DAC and ADC conversions are exact inverses of each other (in this case). Now, the continuous-time signal produced by such DAC doesn’t contain any partials above the Nyquist frequency. Thus, its Fourier integral representation (assuming $T=1$) is

$$x[n]={\int }_{-\pi }^{\pi }X(\omega )e^{j\omega n}\frac{d\omega }{2\pi }$$

and its Laplace integral representation is

$$x[n]={\int }_{\sigma -j\pi }^{\sigma +j\pi }X(s)e^{sn}\frac{ds}{2\pi j}$$

Introducing notation $z=e^s$ and noticing that

$$ds=d(\log z)=\frac{dz}{z}$$

we can rewrite the Laplace integral as

$$x[n]=\oint X(z)z^n\frac{dz}{2\pi jz}$$

(where $X(z)$ is apparently a different function than $X(s)$) where the integration is done counterclockwise along a circle of radius $e^{\sigma }$ centered at the complex plane’s origin:¹

$$z=e^s=e^{\sigma +j\omega }=e^{\sigma }\cdot e^{j\omega }(-\pi \le \omega \le \pi )\text{\hspace{1em}(3.1)}$$

We will refer the representation (3.1) as the z-integral.² The function $X(z)$ is referred to as the z-transform of $x[n]$.

In case of non-unit sampling period $T\ne 1$ the formulas are the same, except that the frequency-related parameters get multiplied by $T$ (or divided by $f_s$), or equivalently, the $n$ index gets multiplied by $T$ in continuous-time expressions:³

$$x[n]={\int }_{-\pi f_s}^{\pi f_s}X(\omega )e^{j\omega Tn}\frac{d\omega }{2\pi }$$

$$x[n]={\int }_{\sigma -j\pi f_s}^{\sigma +j\pi f_s}X(s)e^{sTn}\frac{ds}{2\pi j}$$

$$z=e^{sT}$$

$$x[n]=\oint X(z)z^n\frac{dz}{2\pi jz}(z=e^{\sigma +j\omega T},-\pi f_s\le \omega \le \pi f_s)$$

The notation $z^n$ is commonly used for discrete-time complex exponential signals. A continuous-time signal $x(t)=e^{st}$ is written as $x[n]=z^n$ in discrete- time, where $z=e^{sT}$. The Laplace-integral amplitude coefficient $X(s)$ in $X(s)e^{st}$ then may be replaced by a z-integral amplitude coefficient $X(z)$ such as in $X(z)z^n$.

3.2 Naive integration

The most “interesting” element of analog filter block diagrams is obviously the integrator. The time-discretization for other elements is trivial, so we should concentrate on building the discrete-time models of the analog integrator.

The continuous-time integrator equation is

$$y(t)=y(t_0)+{\int }_{t_0}^{t}x(\tau )d\tau$$

In discrete time we could approximate the integration by a summation of the input samples. Assuming for simplicity $T=1$, we could have implemented a discrete-time integrator as

$$y[n]=y[n_0-1]+\sum _{\nu =n_0}^{n}x[\nu ]$$

We will refer to the above as the naive digital integrator.

A pseudocode routine for this integrator could simply consist of an accumulating assignment:

// perform one sample tick of the integrator
integrator_output := integrator_output + integrator_input;

It takes the current state of the integrator stored in the integrator_output variable and adds the current sample’s value of the integrator_input on top of that.

In case of a non-unit sampling period $T\ne 1$ we have to multiply the accumulated input values by $T$:⁴

// perform one sample tick of the integrator
integrator_output := integrator_output + integrator_input*T;

3.3 Naive lowpass filter

We could further apply this “naive” approach to construct a discrete-time model of the lowpass filter in Fig. 2.2. We will use the naive integrator as a basis for this model.⁵

Let the $x$ variable contain the current input sample of the filter. Considering that the output of the filter in Fig. 2.2 coincides with the output of the integrator, let the $y$ variable contain the integrator state and simultaneously serve as the output sample. As we begin to process the next input sample, the $y$ variable will contain the previous output value. At the end of the processing of the sample (by the filter model) the $y$ variable will contain the new output sample. In this setup, the input value for the integrator is apparently $(x-y){\omega }_c$, thus we simply have

// perform one sample tick of the lowpass filter
y := y + (x-y)*omega_c;

(mind that ${\omega }_c$ must have been scaled to the time units corresponding to the unit sample period!)

A naive discrete-time model of the multimode filter in Fig. 2.13 could have been implemented as:

// perform one sample tick of the multimode filter
hp := x-lp;
lp := lp + hp*omega_c;

where the integrator state is stored in the lp variable.

The above naive implementations (and any other similar naive implementations, for that matter) work reasonably well as long as ${\omega }_c\ll 1$, that is the cutoff must be much lower than the sampling rate. At larger ${\omega }_c$ the behavior of the filter becomes rather strange, ultimately the filter gets unstable. We will now develop some theoretical means to analyse the behavior of the discrete-time filter models, figure out what are the problems with the naive implementations, and then introduce another discretization approach.

3.4 Block diagrams

Let’s express the naive discrete-time integrator in the form of a discrete-time block diagram. The discrete-time block diagrams are constructed from the same elements as continuous-time block diagrams, except that instead of integrators they have unit delays. A unit delay simply delays the signal by one sample. That is the output of a unit delay comes “one sample late” compared to the input. Apparently, the implementation of a unit delay requires a variable, which will be used to store the new incoming value and keep it there until the next sample. Thus, a unit delay element has a state, while the other block diagram elements are obviously stateless. This makes the unit delays in a way similar to the integrators in the analog block diagrams, where the integrators are the only elements with a state.

A unit delay element in a block diagram is denoted as:

$$⟶\boxed{z^{-1}}⟶$$

The reason for the notation $z^{-1}$ will be explained a little bit later. Using a unit delay, we can create a block diagram for our naive integrator (Fig. 3.1). For an arbitrary sampling period we obtain the structure in Fig. 3.2. For an integrator with embedded cutoff gain we can combine the ${\omega }_c$ gain element with the $T$ gain element (Fig. 3.3). Notice that the integrator thereby becomes invariant to the choice of the time units, since ${\omega }_{c}T$ is invariant to this choice.

Figure 3.1: Naive integrator for T = 1.

Figure 3.2: Naive integrator for arbitrary T.

Figure 3.3: Naive integrator with embedded cutoff.

Now let’s construct the block diagram of the naive 1-pole lowpass filter. Recalling the implementation routine:

// perform one sample tick of the lowpass filter
y := y + (x-y)*omega_c;

we obtain the diagram in Fig. 3.4. The $z^{-1}$ element in the feedback from the filter’s output to the leftmost summator is occurring due to the fact that we are picking up the previous value of $y$ in the routine when computing the difference $x-y$.

Figure 3.4: Naive 1-pole lowpass filter (the dashed line denotes the integrator).

This unit delay occurring in the discrete-time feedback is a common problem in discrete-time implementations. This problem is solvable, however it doesn’t make too much sense to solve it for the naive integrator-based models, as the increased complexity doesn’t justify the improvement in sound. We will address the problem of the zero-delay discrete-time feedback later, for now we’ll concentrate on the naive model in Fig. 3.4. This model can be simplified a bit, by combining the two $z^{-1}$ elements into one (Fig. 3.5), so that the block diagram explicitly contains a single state variable (as does its pseudocode counterpart).

Figure 3.5: Naive 1-pole lowpass filter with just one $z^{-1}$ element (the dashed line denotes the integrator).

3.5 Transfer function

Let $x[n]$ and $y[n]$ be respectively the input and the output signals of a unit delay:

$$x[n]⟶\boxed{z^{-1}}⟶y[n]$$

For a complex exponential input $x[n]=e^{sn}=z^n$ we obtain

$$y[n]=e^{s(n-1)}=e^{sn}{e}^{-s}=z^n{z}^{-1}=z^{-1}x[n]$$

That is

$$y[n]=z^{-1}x[n]$$

That is, $z^{-1}$ is the transfer function of the unit delay! It is common to express discrete-time transfer functions as functions of $z$ rather than functions of $s$. The reason is that in this case the transfer functions are nonstrictly proper⁶ rational functions, similarly to the continuous-time case, which is pretty convenient. So, for a unit delay we could write $H(z)=z^{-1}$.

Now we can obtain the transfer function of the naive integrator in Fig. 3.1. Suppose⁷ $x[n]=X(z)z^n$ and $y[n]=Y(z)z^n$, or shortly, $x=X(z)z^n$ and $y=Y(z)z^n$. Then the output of the $z^{-1}$ element is $y{z}^{-1}$. The output of the summator is then $x+y{z}^{-1}$, thus

$$y=x+y{z}^{-1}$$

from where

$$y(1-z^{-1})=x$$

and

$$H(z)=\frac{y}{x}=\frac{1}{1-z^{-1}}$$

This is the transfer function of the naive integrator (for $T=1$).

It is relatively common to express discrete-time transfer functions as rational functions of $z^{-1}$ (like the one above) rather than rational functions of $z$. However, for the purposes of the analysis it is also often convenient to have them expressed as rational functions of $z$ (particularly, for finding their poles and zeros). We can therefore multiply the numerator and the denominator of the above $H(z)$ by $z$, obtaining:

$$H(z)=\frac{z}{z-1}$$

Since $z=e^s$, the frequency response is obtained as $H(e^{j\omega })$. The amplitude and phase responses are $|H(e^{j\omega })|$ and $\arg H(e^{j\omega })$ respectively.⁸

For $T\ne 1$ we obtain

$$H(z)=T\frac{z}{z-1}$$

and, since $z=e^{sT}$, the frequency response is $H(e^{j\omega T})$.

Now let’s obtain the transfer function of the naive 1-pole lowpass filter in Fig. 3.5, where, for the simplicity of notation, we assume $T=1$. Assuming complex exponentials $x=X(z)z^n$ and $y=Y(z)z^n$ we have $x$ and $y{z}^{-1}$ as the inputs of the first summator. Respectively the integrator’s input is ${\omega }_c(x-y{z}^{-1})$. And the integrator output is the sum of $y{z}^{-1}$ and the integrator’s input. Therefore

$$y=y{z}^{-1}+{\omega }_c(x-y{z}^{-1})$$

From where

$$(1-(1-{\omega }_c)z^{-1})y={\omega }_{c}x$$

and

$$H(z)=\frac{y}{x}=\frac{{\omega }_c}{1-(1-{\omega }_c)z^{-1}}=\frac{{\omega }_{c}z}{z-(1-{\omega }_c)}$$

The transfer function for $T\ne 1$ can be obtained by simply replacing ${\omega }_c$ by ${\omega }_{c}T$.

The respective amplitude response is plotted in Fig. 3.6. Comparing it to the amplitude response of the analog prototype we can observe serious deviation closer to the Nyquist frequency. The phase response (Fig. 3.7) has similar deviation problems.

In principle, the amplitude response deviation can be drastically reduced by correcting the filter’s cutoff setting. E.g. one could notice that the second of the amplitude responses in Fig. 3.6 is occurring a bit too far to the right, compared to the analog response (which is what we’re aiming at). Therefore we could achieve a better matching between the two responses by reducing the cutoff setting of the digital filter by a small amount. Depending on the formal definition of the response matching, one could derive an analytical expression for such cutoff correction. There are two main problems with that, though.

Figure 3.6: Amplitude response of a naive 1-pole lowpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs.

Figure 3.7: Phase response of a naive 1-pole lowpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs.

One problem is that many other filters, e.g. a 2-pole resonating lowpass, have more parameters, e.g. not only cutoff but also the resonance, and we potentially may need to correct all of them, which results in much more involved math. This problem is not as critical though, and there are some methods utilizing this approach.

The other problem, though, is the phase response. Looking at Fig. 3.7 it seems that no matter how we try to correct the filter’s cutoff, the phase response will be always zero at Nyquist, whereas we would desire something close to $-\pi /2$. The effects of the deviation of the filter’s phase response are mostly quite subtle. Therefore it’s somewhat difficult to judge how critical the phase devations might be.⁹ However there’s one absolutely objective and major issue associated with the phase deviations. Attempting to mix outputs of two filters with some deviations in either or both of the amplitude and phase responses may easily lead to unexpected and undesired results. For that reason in this book we will concentrate on a different method which is much more robust in this respect.

Poles and zeros

Discrete-time block diagrams are differing from continuous-time block diagrams only by having $z^{-1}$ elements instead of integrators. Recalling that the transfer function of an integrator is $s^{-1}$, we conclude that from the formal point of view the difference is purely notational.

Now, the transfer functions of continuous-time block diagrams are non- strictly proper rational functions of $s$. Respectively, the transfer functions of discrete-time block diagrams are nonstrictly proper rational functions of $z$.

Thus, discrete-time transfer functions will have poles and zeros in a way similar to continuous-time transfer functions. Similarly to continuous-time transfer functions, the poles will define the stability of a linear time-invariant filter. Consider that $z=e^{sT}$ and recall the stability criterion $\mathrm{Re}s<0$ (where $s=p_n$, where $p_n$ are the poles). Apparently, $\mathrm{Re}s<0⟺|z|<1$. We might there- fore intuitively expect the discrete-time stability criterion to be $|p_n|<1$ where $p_n$ are the discrete-time poles. This is indeed the case, a linear time-invariant difference system¹⁰ is stable if and only if all its poles are located inside the unit circle. We will give more detail about the mechanisms behind this in the discussion of the discrete-time transient response in Sections 3.12 and 7.13.

3.6 Trapezoidal integration

Instead of naive integration, we could attempt using the trapezoidal integration method ($T=1$):

// perform one sample tick of the integrator
integrator_output := integrator_output +
    (integrator_input + previous_integrator_input)/2;
previous_integrator_input := integrator_input;

Notice that now we need two state variables per integrator: integrator_output and previous_integrator_input. The block diagram of a trapezoidal integrator is shown in Fig. 3.8. We’ll refer to this integrator as a direct form I trapezoidal integrator. The reason for this term will be explained later.

Figure 3.8: Direct form I trapezoidal integrator (T = 1).

We could also construct a trapezoidal integrator implementation with only a single state variable. Consider the expression for the trapezoidal integrator’s output:

$$y[n]=y[n_0-1]+\sum _{\nu =n_0}^n\frac{x[\nu -1]+x[\nu ]}{2}\text{\hspace{1em}(3.2)}$$

Suppose $y[n_0-1]=0$ and $x[n_0-1]=0$, corresponding to a zero initial state (recall that both $y[n_0-1]$ and $x[n_0-1]$ are technically stored in the $z^{-1}$ elements). Then

$$\begin{array}{rl}y[n]& =\sum _{\nu =n_0}^n\frac{x[\nu -1]+x[\nu ]}{2}=\frac{1}{2}\left(\sum _{\nu =n_0}^{n}x[\nu -1]+\sum _{\nu =n_0}^{n}x[\nu ]\right)=\\ & =\frac{1}{2}\left(\sum _{\nu =n_0+1}^{n}x[\nu -1]+\sum _{\nu =n_0}^{n}x[\nu ]\right)=\frac{1}{2}\left(\sum _{\nu =n_0}^{n-1}x[\nu ]+\sum _{\nu =n_0}^{n}x[\nu ]\right)=\\ & =\frac{u[n-1]+u[n]}{2}\end{array}$$

where

$$u[n]=\sum _{\nu =n_0}^{n}x[\nu ]$$

Now notice that $u[n]$ is the output of a naive integrator, whose input signal is $x[n]$. At the same time $y[n]$ is the average of the previous and the current output values of the naive integrator. This can be implemented by the structure in Fig. 3.9. Similar considerations apply for nonzero initial state. We’ll refer to the integrator in Fig. 3.9 as a direct form II or canonical trapezoidal integrator. The reason for this term will be explained later.

We can develop yet another form of the bilinear integrator with a single state variable. Let’s rewrite (3.2) as

$$y[n]=y[n_0-1]+\frac{x[n_0-1]}{2}+\sum _{\nu =n_0}^{n-1}x[\nu ]+\frac{x[n]}{2}$$

and let

$$u[n-1]=y[n]-\frac{x[n]}{2}=y[n_0-1]+\frac{x[n_0-1]}{2}+\sum _{\nu =n_0}^{n-1}x[\nu ]$$

Figure 3.9: Direct form II (canonical) trapezoidal integrator (T = 1).

Notice that

$$y[n]=u[n-1]+\frac{x[n]}{2}\text{\hspace{1em}(3.3a)}$$

and

$$u[n]=u[n-1]+x[n]=y[n]+\frac{x[n]}{2}\text{\hspace{1em}(3.3b)}$$

Expressing (3.3a) and (3.3b) in a graphical form, we obtain the structure in Fig. 3.10. We’ll refer to the integrator in Fig. 3.10 as a transposed direct form II or transposed canonical trapezoidal integrator. The reason for this term will be explained later.

Figure 3.10: Transposed direct form II (transposed canonical) trapezoidal integrator (T = 1).

The positioning of the 1/2 gain prior to the integrator in Fig. 3.10 is quite convenient, because we can combine the 1/2 gain with the cutoff gain into a single gain element. In case of an arbitrary sampling period we could also include the $T$ factor into the same gain element, thus obtaining the structure in Fig. 3.11. A similar trick can be performed for the other two integrators, if we move the 1/2 gain element to the input of the respective integrator. Since the integrator is a linear time-invariant system, this doesn’t affect the integrator’s behavior in a slightest way.

Typically one would prefer the direct form II integrators to the direct form I integrator, because the former have only one state variable. In this book we will mostly use the transposed direct form II integrator, because this is resulting in slightly simpler zero-delay feedback equations and also offers a nice possibility for the internal saturation in the integrator.

Figure 3.11: Transposed direct form II (transposed canonical) trapezoidal integrator with “embedded” cutoff gain.

The transfer functions of all three integrators are identical. Let’s obtain e.g. the transfer function of the transposed canonical integrator (in Fig. 3.10). Assuming signals of the exponential form $z^n$, we can drop the index $[n]$, under- standing it implicity, while the index $[n-1]$ will be replaced by the multiplication by $z^{-1}$. Then (3.3) turn into

$$\begin{array}{rl}y& =u{z}^{-1}+\frac{x}{2}\\ u& =y+\frac{x}{2}\end{array}$$

Substituting the second equation into the first one we have

$$\begin{array}{rl}y& =\left(y+\frac{x}{2}\right)z^{-1}+\frac{x}{2}\\ yz& =y+\frac{x}{2}+\frac{x}{2}z\\ y(z-1)& =\frac{x}{2}(z+1)\end{array}$$

and the transfer function of the trapezoidal integrator is thus

$$H(z)=\frac{y}{x}=\frac{1}{2}\cdot \frac{z+1}{z-1}$$

For an arbitrary $T$ one has to multiply the result by $T$, to take the respective gain element into account:

$$H(z)=\frac{T}{2}\cdot \frac{z+1}{z-1}$$

If also the cutoff gain is included, we obtain

$$H(z)=\frac{{\omega }_{c}T}{2}\cdot \frac{z+1}{z-1}$$

One can obtain the same results for the other two integrators.

What is so special about this transfer function, that makes the trapezoidal integrator so superior to the naive one, is to be discussed next.

3.7 Bilinear transform

Suppose we take an arbitrary continuous-time block diagram, like the familiar lowpass filter in Fig. 2.2 and replace all continuous-time integrators by discrete- time trapezoidal integrators. On the transfer function level, this will correspond to replacing all $s^{-1}$ with $\frac{T}{2}\cdot \frac{z+1}{z-1}$. That is, technically we perform a substitution

$$s^{-1}=\frac{T}{2}\cdot \frac{z+1}{z-1}$$

in the transfer function expression.

It would be more convenient to write this substitution explicitly as

$$s=\frac{2}{T}\cdot \frac{z-1}{z+1}\text{\hspace{1em}(3.4)}$$

The substitution (3.4) is referred to as the bilinear transform, or shortly BLT. For that reason we can also refer to trapezoidal integrators as BLT integrators. Let’s figure out, how does the bilinear transform affect the frequency response of the filter, that is, what is the relationship between the original continuous- time frequency response prior to the substitution and the resulting discrete-time frequency response after the substitution.

Let $H_a(s)$ be the original continuous-time transfer function. Then the re- spective discrete-time transfer function is

$$H_d(z)=H_a\left(\frac{2}{T}\cdot \frac{z-1}{z+1}\right)\text{\hspace{1em}(3.5)}$$

Respectively, the discrete-time frequency response is

$$\begin{array}{rl}{H}_d(e^{j\omega T})& =H_a\left(\frac{2}{T}\cdot \frac{e^{j\omega T}-1}{e^{j\omega T}+1}\right)=H_a\left(\frac{2}{T}\cdot \frac{e^{j\omega T/2}-e^{-j\omega T/2}}{e^{j\omega T/2}+e^{-j\omega T/2}}\right)=\\ & =H_a\left(\frac{2}{T}j\tan \frac{\omega T}{2}\right)\end{array}$$

Notice that $H_a(s)$ in the last expression is evaluated on the imaginary axis!!! That is, the bilinear transform maps the imaginary axis in the $s$-plane to the unit circle in the $z$-plane! Now, $H_a\left(\frac{2}{T}j\tan \frac{\omega T}{2}\right)$ is the analog frequency response evaluated at $\frac{2}{T}\tan \frac{\omega T}{2}$. That is, the digital frequency response at $\omega$ is equal to the analog frequency response at $\frac{2}{T}\tan \frac{\omega T}{2}$. This means that the analog fre- quency response in the range $0\le \omega <+\infty$ is mapped into the digital frequency range $0\le \omega T<\pi$ ($0\le \omega <\pi f_s$), that is from zero to Nyquist!¹¹ Denoting the analog frequency as ${\omega }_a$ and the digital frequency as ${\omega }_d$ we can express the argument mapping of the frequency response function as

$${\omega }_a=\frac{2}{T}\tan \frac{{\omega }_{d}T}{2}\text{\hspace{1em}(3.6)}$$

or, in a more symmetrical way

$$\frac{{\omega }_{a}T}{2}=\tan \frac{{\omega }_{d}T}{2}\text{\hspace{1em}(3.7)}$$

Notice that for frequencies much smaller that Nyquist frequency we have $\omega T\ll 1$ and respectively ${\omega }_a\approx {\omega }_d$.

This is what is so unique about the bilinear transform. It simply warps the frequency range $[0,+\infty )$ into the zero-to-Nyquist range, but otherwise doesn’t change the frequency response at all! Considering in comparison a naive inte- grator, we would have obtained:

$$\begin{array}{rl}{s}^{-1}& =\frac{z}{z-1}\\ s& =\frac{z-1}{z}\end{array}\text{\hspace{1em}(3.8)}$$

$$H_d(z)=H_a\left(\frac{z-1}{z}\right)$$

$$H_d(e^{j\omega })=H_a\left(\frac{e^{j\omega }-1}{e^{j\omega }}\right)=H_a\left(1-e^{-j\omega }\right)$$

which means that the digital frequency response is equal to the analog transfer function evaluated on a circle of radius 1 centered at $s=1$. This hardly defines a clear relationship between the two frequency responses.

So, by simply replacing the analog integrators with digital trapezoidal in- tegrators, we obtain a digital filter whose frequency response is essentially the same as the one of the analog prototype, except for the frequency warping. Particularly, the relationship between the amplitude and phase responses of the filter is fully preserved, which is particularly highly important if the filter is to be used as a building block in a larger filter. Very close to perfect!

Furthermore, the bilinear transform maps the left complex semiplane in the $s$-domain into the inner region of the unit circle in the $z$-domain. Indeed, let’s obtain the inverse bilinear transform formula. From (3.4) we have

$$(z+1)\frac{sT}{2}=z-1$$

from where

$$1+\frac{sT}{2}=z\left(1-\frac{sT}{2}\right)$$

and

$$z=\frac{1+\frac{sT}{2}}{1-\frac{sT}{2}}\text{\hspace{1em}(3.9)}$$

The equation (3.9) defines the inverse bilinear transform. Now, if $\mathrm{Re}s<0$, then, obviously

$$\left|1+\frac{sT}{2}\right|<\left|1-\frac{sT}{2}\right|$$

and $|z|<1$. Thus, the left complex semiplane in the $s$-plane is mapped to the inner region of the unit circle in the $z$-plane. In the same way one can show that the right complex semiplane is mapped to the outer region of the unit circle. And the imaginary axis is mapped to the unit circle itself. Comparing the stability criterion of analog filters (the poles must be in the left complex semiplane) to the one of digital filters (the poles must be inside the unit circle), we conclude that the bilinear transform exactly preserves the stability of the filters!

In comparison, for a naive integrator replacement we would have the follow- ing. Inverting the (3.8) substitution we obtain

$$\begin{array}{rl}sz& =z-1\\ z(1-s)& =1\end{array}$$

and

$$z=\frac{1}{1-s}$$

Assuming $\mathrm{Re}s<0$ and considering that in this case

$$\left|z-\frac{1}{2}\right|=\left|\frac{1}{1-s}-\frac{1}{2}\right|=\left|\frac{1-\frac{1}{2}+\frac{s}{2}}{1-s}\right|=\left|\frac{1}{2}\cdot \frac{1+s}{1-s}\right|<\frac{1}{2}$$

we conclude that the left semiplane is mapped into a circle of radius 0.5 cen- tered at $z=0.5$. So the naive integrator overpreserves the stability, which is not nice, since we would rather have digital filters behaving as closely to their analog prototypes as possible. Considering that this comes in a package with a poor frequency response transformation, we should rather stick with trapezoidal integrators.

So, let’s replace e.g. the integrator in the familiar lowpass filter structure in Fig. 2.2 with a trapezoidal integrator. Performing the integrator replacement, we obtain the structure in Fig. 3.12.¹² We will refer to the trapezoidal integrator replacement method as the topology-preserving transform (TPT) method. This term will be explained and properly introduced later. For now, before we simply attempt to implement the structure in Fig. 3.12 in code, we should become aware of a few further issues.

3.8 Cutoff prewarping

Suppose we are using the lowpass filter structure in Fig. 3.12 and we wish to have its cutoff at ${\omega }_c$. If we however simply put this ${\omega }_c$ parameter into the respective integrator gain element ${\omega }_{c}T/2$, the frequency response itself and, specifically, its value at the cutoff will be different from the expected one. Fig. 3.13 illustrates. The $-3$dB level is specifically highlighted in Fig. 3.13, since this is the amplitude response value of the 1-pole lowpass filter at the cutoff, thereby aiding the visual identification of the cutoff point on the response curves.¹³

Apparently, the difference between analog and digital response is occuring due to the warping of the frequency axis (3.6). We would like to estimate the frequency error introduced by the warping. To simplify the further discussion let’s rewrite (3.6) as a mapping function $\mu (\omega )$:

Figure 3.12: 1-pole TPT lowpass filter (the dashed line denotes the trapezoidal integrator).

Figure 3.13: Amplitude response of an unprewarped bilinear-transformed 1-pole lowpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

$${\omega }_a=\mu ({\omega }_d)=\frac{2}{T}\tan \frac{{\omega }_{d}T}{2}\text{\hspace{1em}(3.10)}$$

Now, given some desired analog response, we could take some point ${\omega }_a$ on this response and ask ourselves, where is the same point located on the digital response. According to (3.10), it is located at ${\mu }^{-1}({\omega }_a)$ (where ${\mu }^{-1}$ is the function inverse of $\mu$). Thus the ratio of the actual and desired frequencies is ${\mu }^{-1}({\omega }_a)/{\omega }_a$, or, in the octave scale:

$$\Delta P={\log }_2\frac{{\mu }^{-1}({\omega }_a)}{{\omega }_a}$$

The solid curve in Fig. 3.14 illustrates (note that Fig. 3.14 labels the $\Delta P$ axis in semitones).

Figure 3.14: Bilinear transform’s detuning of analog frequencies plotted against analog frequency (solid curve) or digital frequency (dashed curve). Sampling rate 44.1kHz.

We could also express the detuning of analog frequencies in terms of the corresponding digital frequency. Given the digital frequency response, we take some point ${\omega }_d$ and ask ourselves, where is the same point located on the analog response. According to (3.10), it is located at $\mu ({\omega }_d)$ and thus the frequency ratio is ${\omega }_d/\mu ({\omega }_d)$, respectively

$$\Delta P={\log }_2\frac{{\omega }_d}{\mu ({\omega }_d)}$$

The dashed curve in Fig. 3.14 illustrates. Note that we are not talking about how much the specified digital frequency will be detuned (because digital frequencies are not getting detuned, they are already where they are), it’s still about how much the corresponding analog frequency will be detuned.

Thus, given an analog filter with a frequency response $H_a(j\omega )$, its digital counterpart will have its frequencies detuned as shown in Fig. 3.14. Particularly, the cutoff point, instead of being at the specified frequency $\omega ={\omega }_c$, will be at

$${\omega }_d={\mu }^{-1}({\omega }_c)\text{\hspace{1em}(3.11)}$$

In principle, one could argue that the frequency response change in Fig. 3.13 is not that drastic and could be tolerated, especially since the deviation occurs mostly in the high frequency range, which is not the most audible part of the frequency spectrum. This might have been the case with the 1-pole lowpass filter, however for other filters with more complicated amplitude responses it won’t be as acceptable. Fig. 3.15 illustrates the frequency error for a 2-pole resonating lowpass filter. The resonance peaks (which occur close to the filter’s cutoff) are very audible and so would be their detuning, which according to Fig. 3.14 is in the range of semitones. Particularly, at 16kHz the dashed curve in Fig. 3.14 shows a detuning of ca. 1 octave, meaning that we would have a resonance at this point when it should have been occurring at ca. 32kHz.

Figure 3.15: Amplitude response of an unprewarped bilinear-transformed resonating 2-pole lowpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

Prewarping at cutoff

As a general rule (to which there are exceptions), we would like the cutoff point of the filter to be positioned exactly at the specified cutoff frequency ${\omega }_c$. In this regard we could notice that if we used a different cutoff value

$${\tilde{\omega }}_c=\mu ({\omega }_c)\text{\hspace{1em}(3.12)}$$

then (3.11) would give

$${\omega }_d={\mu }^{-1}({\tilde{\omega }}_c)={\mu }^{-1}(\mu ({\omega }_c))={\omega }_c$$

and the cutoff point would be exactly where we wanted it to be. Fig. 3.16 illus- trates. The cutoff correction (3.12) is a standard technique used in combination with the bilinear transform. It is referred to as cutoff prewarping.

Technically, cutoff prewarping means that we use ${\tilde{\omega }}_c$ instead of ${\omega }_c$ in the gains of the filter’s integrators. However, the integrator gains are not exactly ${\omega }_c$ but rather ${\omega }_{c}T/2$. From (3.12) and (3.10) we have

$$\frac{{\tilde{\omega }}_{c}T}{2}=\frac{2}{T}\tan \frac{{\omega }_{c}T}{2}\cdot \frac{T}{2}=\tan \frac{{\omega }_{c}T}{2}\text{\hspace{1em}(3.13)}$$

Thus, we can directly apply (3.13) to compute the prewarped gains ${\tilde{\omega }}_{c}T/2$. Note that (3.13) is essentially identical to (3.7).

The cutoff prewarping redistributes the frequency error shown in Fig. 3.14. In the absence of prewarping the error was zero at $\omega =0$ and monotonically growing as $\omega$ increases. With the cutoff prewarping the error is zero at $\omega ={\omega }_c$ instead and grows further away from this point.

Indeed, let $H(j\omega )$ be unit-cutoff analog response of the filter in question. And let’s pick up an analog frequency ${\omega }_a$ and find the respective detuning. The correct frequency response at ${\omega }_a$ is

$$H_a({\omega }_a)=H(j{\omega }_a/{\omega }_c)\text{\hspace{1em}(3.14a)}$$

Figure 3.16: Amplitude response of a prewarped bilinear-transformed 1-pole lowpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

On the other hand, given the prewarped cutoff ${\tilde{\omega }}_c=\mu ({\omega }_c)$, the digital frequency response at some frequency ${\omega }_d$ is

$$H_d(e^{j{\omega }_d})=H(j\mu ({\omega }_d)/{\tilde{\omega }}_c)=H(j\mu ({\omega }_d)/\mu ({\omega }_c))\text{\hspace{1em}(3.14b)}$$

We want to find such ${\omega }_d$ that the arguments of $H(j\omega )$ in (3.14a) and (3.14b) are identical:

$$\frac{\mu ({\omega }_d)}{\mu ({\omega }_c)}=\frac{{\omega }_a}{{\omega }_c}\text{\hspace{1em}(3.15)}$$

From where

$${\omega }_d={\mu }^{-1}\left({\omega }_a\frac{\mu ({\omega }_c)}{{\omega }_c}\right)$$

The solid curves in Fig. 3.17 illustrate the respective detuning ${\omega }_a/{\omega }_d$ (in semi- tones) of analog frequencies.

Alternatively from (3.15) we could express ${\omega }_a$ as a function of ${\omega }_d$:

$${\omega }_a=\frac{{\omega }_c}{\mu ({\omega }_c)}\mu ({\omega }_d)$$

thus expressing the analog frequency detuning ${\omega }_a/{\omega }_d$ as a function of ${\omega }_d$. The dashed curves in Fig. 3.17 illustrate.

Apparently, the maximum error to the left of $\omega ={\omega }_c$ is attained at $\omega =0$. Letting ${\omega }_a\to 0$ and, equivalently, ${\omega }_d\to 0$ we have $\mu ({\omega }_d)\sim {\omega }_d$ and (3.15) turns into

$$\frac{{\omega }_d}{\mu ({\omega }_c)}\sim \frac{{\omega }_a}{{\omega }_c}$$

or

$$\frac{{\omega }_d}{{\omega }_a}\sim \frac{\mu ({\omega }_c)}{{\omega }_c}$$

and thus the detuning at $\omega =0$ is

$$\Delta P{|}_{\omega =0}={\log }_2\frac{\mu ({\omega }_c)}{{\omega }_c}\text{\hspace{1em}(3.16)}$$

Figure 3.17: Prewarped bilinear transform’s detuning of analog frequencies plotted against analog frequency (solid curves) or digital frequency (dashed curves). Different curves correspond to prewarping at different cutoff frequencies. Sampling rate 44.1kHz.

Other prewarping points

We have just developed the prewarping technique from the condition that the cutoff point must be preserved by the mapping (3.6). However instead we could have required any other point ${\omega }_p$ to be preserved.

Given the cutoff ${\omega }_c$, the analog response at ${\omega }_p$ is

$$H_a({\omega }_p)=H(j{\omega }_p/{\omega }_c)\text{\hspace{1em}(3.17a)}$$

On the other hand, given the prewarped cutoff ${\tilde{\omega }}_c$ (we want to prewarp at a different point now, therefore we don’t know yet, what is the relationship between ${\omega }_c$ and ${\tilde{\omega }}_c$) the digital frequency response at ${\omega }_p$ is

$$H_d(e^{j{\omega }_p})=H(j\mu ({\omega }_p)/{\tilde{\omega }}_c)\text{\hspace{1em}(3.17b)}$$

We want to find such ${\tilde{\omega }}_c$ that the arguments of $H(j\omega )$ in (3.17a) and (3.17b) are identical:

$$\frac{{\omega }_p}{{\omega }_c}=\frac{\mu ({\omega }_p)}{{\tilde{\omega }}_c}$$

and

$${\tilde{\omega }}_c=\frac{\mu ({\omega }_p)}{{\omega }_p}{\omega }_c\text{\hspace{1em}(3.18)}$$

Equation (3.18) is the generalized prewarping formula, where ${\omega }_p$ is the frequency response point of zero detuning. We will refer to ${\omega }_p$ as the prewarping point.

According to (3.18) prewarping at ${\omega }_p$ simply means that the cutoff should be multiplied by $\mu ({\omega }_p)/{\omega }_p$. At ${\omega }_p={\omega }_c$ this multiplication reduces to (3.12).

In order to find the detuning at other frequencies, notice that equations (3.14) turn into

$$H_a({\omega }_a)=H(j{\omega }_a/{\omega }_c)$$

$$H_d(e^{j{\omega }_d})=H(j\mu ({\omega }_d)/{\tilde{\omega }}_c)=H\left(j\frac{{\omega }_p\mu ({\omega }_d)}{{\omega }_c\mu ({\omega }_p)}\right)$$

from where, equating the arguments of $H(j\omega )$:

$$\frac{{\omega }_p\mu ({\omega }_d)}{{\omega }_c\mu ({\omega }_p)}=\frac{{\omega }_a}{{\omega }_c}$$

we have

$$\frac{\mu ({\omega }_d)}{\mu ({\omega }_p)}=\frac{{\omega }_a}{{\omega }_p}\text{\hspace{1em}(3.19)}$$

Equation (3.19) is identical to (3.15) except that it has ${\omega }_p$ in place of ${\omega }_c$. Thus we could reuse the results of the previous analysis of the analog frequency detuning. In particular Fig. 3.17 fully applies, different curves corresponding to different prewarping points. At the same time, (3.16) simply turns to

$$\Delta P{|}_{\omega =0}={\log }_2\frac{\mu ({\omega }_p)}{{\omega }_p}\text{\hspace{1em}(3.20)}$$

Bounded cutoff prewarping

Even though cutoff prewarping is an absolutely standard technique and is often used without any second thought, the need for a different choice of the prewarping point is actually not as exotic as it might seem. Consider e.g. the amplitude response of a 1-pole highpass filter prewarped by (3.12), shown in Fig. 3.18. One can notice a huge discrepancy between analog and digital amplitude responses occurring well into the audible frequency range $[0,16\text{kHz}]$. The error is getting particularly bad at cutoffs above 16kHz. In comparison, the responses of unprewarped filters in Fig. 3.19 even look kind of better, especially if only the audible frequency range is considered. This would be even more so, if higher sampling rates were involved, where the audible range error in Fig. 3.19 would become smaller, while the same error in Fig. 3.18 can still get as large, given a sufficiently high cutoff value.

Apparently, the large error within the audible range in Fig. 3.18 is due to the detuning error illustrated in Fig. 3.17. This error wasn’t as obvious in the case of the 1-pole lowpass filter, since this filter’s amplitude response is almost constant to the left of the cutoff point. On the other hand, highpass filter’s amplitude response is changing noticeably in the same area, which makes the detuning error is made much more promiment.

Does this suggest that we shouldn’t use cutoff prewarping with highpass filters? In principle this is engineer’s decision. However consider the unprewarped resonating 2-pole highpass filter’s amplitude response in Fig. 3.20. As with the responating lowpass, the resonance peak detuning is quite prominent here. Also the difference in the response value is magnified around the resonance point. All in all, we’d rather prewarp the cutoff (Fig. 3.21) and tolerate the detuning to the left of ${\omega }_c$.

Figure 3.18: Amplitude response of a prewarped bilinear-transformed 1-pole highpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

Figure 3.19: Amplitude response of an unprewarped bilinear-transformed 1-pole highpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

However notice, that as the cutoff peak is getting out of the audible range, we stop caring, where exactly it is positioned, since it can’t be heard anyway. So, why should we then tolerage the error in the audible range which continues to increase even faster? Instead, at this moment we could fix the prewarping point to the upper boundary of the audible range:

$${\omega }_p=\{\begin{array}{ll}{\omega }_c& \text{if }{\omega }_c\le {\omega }_{\max }\\ {\omega }_{\max }& \text{if }{\omega }_c\ge {\omega }_{\max }\end{array}$$

or simply

$${\omega }_p=\min \{{\omega }_c,{\omega }_{\max }\}\text{\hspace{1em}(3.21)}$$

where ${\omega }_{\max }$ is some point around 16kHz. At least then the detuning in the audible range won’t grow any further than it is at ${\omega }_c={\omega }_{\max }$ (Fig. 3.22). The picture gets even better at higher sampling rates (Fig. 3.23).

Figure 3.20: Amplitude response of an unprewarped bilinear-transformed resonating 2-pole highpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

Figure 3.21: Amplitude response of a prewarped bilinear-transformed resonating 2-pole highpass filter for a number of different cutoffs. Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

Figure 3.22: Effect of cutoff prewarping bounded at 16kHz. Sampling rate 44.1kHz.

Figure 3.23: Effect of cutoff prewarping bounded at 16kHz. Sampling rate 88.1kHz.

Substituting (3.21) into (3.18) we obtain

$${\tilde{\omega }}_c=\frac{\mu (\min \{{\omega }_c,{\omega }_{\max }\})}{\min \{{\omega }_c,{\omega }_{\max }\}}{\omega }_c=\{\begin{array}{ll}\mu ({\omega }_c)& \text{if }{\omega }_c\le {\omega }_{\max }\\ \frac{\mu ({\omega }_{\max })}{{\omega }_{\max }}{\omega }_c& \text{if }{\omega }_c\ge {\omega }_{\max }\end{array}\text{\hspace{1em}(3.22)}$$

The mapping defined by (3.22) is shown in Fig. 3.24. Note that thereby we become able to specify the cutoffs beyond Nyquist, and actually to specify arbitrarily large cutoffs, since the new mapping curve is crossing the former vertical asymptote at ${\omega }_c=\pi /2$.

Figure 3.24: Bounded cutoff prewarping. The thick dashed line shows the unbounded prewarping continuation. The thin oblique dashed line is the continuation of the straight part of the prewarping curve. The black dot marks the breakpoint of the prewarping curve.

The breakpoint at ${\omega }_c={\omega }_{\max }$ in Fig. 3.24 can be somewhat unexpected. In order to understand the mechanism behind its appearance suppose ${\omega }_c$ is varying with time. Using (3.18) we compute the time derivative of ${\tilde{\omega }}_c$:

$${\dot{\tilde{\omega }}}_c={\omega }_c\frac{d}{dt}\frac{\mu ({\omega }_p)}{{\omega }_p}+\frac{\mu ({\omega }_p)}{{\omega }_p}{\dot{\omega }}_c$$

As ${\omega }_c$ becomes larger than ${\omega }_{\max }$ the first term suddenly disappears and there is a jump in ${\dot{\tilde{\omega }}}_c$. In other words, the variation of the prewarping point makes its own contribution to ${\dot{\tilde{\omega }}}_c$. As soon as the variation stops, the respective contribution abruptly disappears.

The frequency detunings occurring in case of (3.22) can be found directly from Fig. 3.17, keeping in mind that (3.22) is simply another expression of (3.21).

Continuous-speed prewarping

The breakpoint occurring in Fig. 3.24 might be undesirable if the cutoff is being modulated, since there can be a sudden change of the perceived modulation speed as the cutoff traverses through the prewarping breakpoint. For that reason it might be desirable to smooth the breakpoint in one way or the other. The simplest approach would be to continue the curve as a tangent line after the breakpoint:

$${\tilde{\omega }}_c=\{\begin{array}{ll}\mu ({\omega }_c)& \text{if }{\omega }_c\le {\omega }_{\max }\\ \mu ({\omega }_{\max })+({\omega }_c-{\omega }_{\max }){\mu }^{\prime }({\omega }_{\max })& \text{if }{\omega }_c\ge {\omega }_{\max }\end{array}\text{\hspace{1em}(3.23)}$$

where

$${\mu }^{\prime }(\omega )=\frac{d}{d\omega }\left(\frac{2}{T}\tan \frac{\omega T}{2}\right)=\frac{1}{{\cos }^2\frac{\omega T}{2}}=1+{\tan }^2\frac{\omega T}{2}=1+{\left(\frac{T}{2}\mu (\omega )\right)}^2$$

is the derivative of $\mu (\omega )$. Note, however, that this will no longer keep ${\omega }_{\max }$ as the prewarping point and the situation would be something in between (3.12) and (3.22).

At ${\omega }_c\to \infty$ from (3.23) we have

$${\tilde{\omega }}_c=\mu ({\omega }_{\max })+({\omega }_c-{\omega }_{\max }){\mu }^{\prime }({\omega }_{\max })\sim {\mu }^{\prime }({\omega }_{\max }){\omega }_c\text{\hspace{1em}(3.24)}$$

Comparing the right-hand side of (3.24) to (3.18) we obtain the equation for the effective prewarping point at infinity:

$$\frac{\mu ({\omega }_p)}{{\omega }_p}={\mu }^{\prime }({\omega }_{\max })\text{\hspace{1em}(3.25)}$$

In principle ${\omega }_p$ can be found from (3.25), however we are not so much interested in how far off will be the prewarping point, as in the estimation of the associated increase in detuning. By (3.20)

$$\Delta P{|}_{\omega =0}={\log }_2\frac{\mu ({\omega }_p)}{{\omega }_p}={\log }_2{\mu }^{\prime }({\omega }_{\max })(\text{at }{\omega }_c=\infty )$$

which for ${\omega }_{\max }=16$kHz at 44.1kHz sampling rate gives ca. 2.5 octaves, while at 88.2kHz sampling rate it gives only about 0.5 octave. In comparison, at ${\omega }_c={\omega }_{\max }$ (and respectively ${\omega }_p={\omega }_{\max }$) we would have a smaller value

$$\Delta P{|}_{\omega =0}={\log }_2\frac{\mu ({\omega }_p)}{{\omega }_p}={\log }_2\frac{\mu ({\omega }_{\max })}{{\omega }_{\max }}$$

which for ${\omega }_{\max }=16$kHz at 44.1kHz sampling rate gives ca. 1 octave, while at 88.2kHz sampling rate it gives only about 2 semitones.

We have therefore found the effect of (3.23) on the detuning occurring at $\omega =0$ for ${\omega }_c\to \infty$. It would also be nice to estimate the same effect at $\omega ={\omega }_{\max }$. By (3.19)

$$\frac{\mu ({\omega }_d)}{{\omega }_a}=\frac{\mu ({\omega }_p)}{{\omega }_p}$$

which by (3.25) becomes

$$\frac{\mu ({\omega }_d)}{{\omega }_a}={\mu }^{\prime }({\omega }_{\max })$$

Letting ${\omega }_d={\omega }_{\max }$ we have

$${\omega }_a=\frac{\mu ({\omega }_{\max })}{{\mu }^{\prime }({\omega }_{\max })}=\frac{\frac{2}{T}\tan \frac{{\omega }_{\max }T}{2}}{{\cos }^{-2}\left(\frac{{\omega }_{\max }T}{2}\right)}=\frac{1}{T}\sin ({\omega }_{\max }T)$$

and the detuning itself is the logarithm of the ratio

$$\frac{{\omega }_a}{{\omega }_d}=\frac{{\omega }_a}{{\omega }_{\max }}=\frac{\sin ({\omega }_{\max }T)}{{\omega }_{\max }T}=\mathrm{sinc}({\omega }_{\max }T)$$

For ${\omega }_{\max }=16$kHz at 44.1kHz sampling rate this gives ca. $-1.5$ octaves, while at 88.1kHz it gives ca. $-4$ semitones.

Since (as illustrated by Fig. 3.17) the detuning is a monotonic function of $\omega$, at ${\omega }_c=\infty$ we are having the audible range detuning error in the range of ca. $[-1.5,2.5]$ octaves at 44.1kHz sampling rate and in the range of ca. $[-4,6]$ semitones at 88.2kHz sampling rate. Therefore it is more or less balanced out, although being a bit larger at higher frequencies.

Apparently, more elaborate ways of smoothing the breakpoint in Fig. 3.24 may be designed, but we won’t cover them here as the options are almost infinite.

Prewarping of systems of filters

According to (3.18), prewarping is nothing more than a multiplication of the cutoff gains of the integrators by $\mu ({\omega }_p)/{\omega }_p$, where ${\omega }_p$ is the prewarping point.

Suppose we are having a system consisting of several filters connected together. When prewarping filters in such system, it would be a good idea to choose a common prewarping point for all filters. In this case we multiply their cutoffs by one and the same coefficient. Thereby their amplitude and phase responses are shifted by one and the same amount (in the logarithmic frequency axis), and they all retain the frequency response relationships which existed between them prior to prewarping. Effectively this is the same as changing the “common cutoff” of the filter system, and the frequency response of the entire system is simply shifted horizontally by the same amount, fully retaining its shape.

On the other hand by prewarping them independently we shift the frequency response of each filter differently from the others and the amplitude and phase relationships between those are thereby destroyed. Therefore the amplitude and phase response shapes of the entire system of filters are not preserved.

As an illustration consider a parallel connection of a 1-pole lowpass and a 1-pole highpass filter, the lowpass cutoff being ${\omega }_c/2$, the highpass cutoff being $2{\omega }_c$ where ${\omega }_c$ is the formal cutoff of the system:

The transfer function of such system (written in the unit-cutoff form) is

$$H(s)=\frac{1}{1+2s}+\frac{s/2}{1+s/2}$$

The result of prewarping the lowpass and the highpass separately at their respective cutoffs ${\omega }_c/2$ and $2\omega$ is shown in Fig. 3.25. Actually for the ${\omega }_c=11.025$kHz and ${\omega }_c=16$kHz the highpass cutoff $2{\omega }_c$ needed to be clipped prior to prewaring, since it is equal or exceeds Nyquist and cannot be directly prewarped by (3.12). Compare to Fig. 3.26 where the prewarping of both filters has been done at the common point ${\omega }_c$.

Figure 3.25: Separate prewarping of system components (for a number of different cutoffs). Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

Figure 3.26: Common-point prewarping of system components (for a number of different cutoffs). Dashed curves represent the respective analog filter responses for the same cutoffs. Sampling rate 44.1kHz.

Other prewarping techniques

With 1-pole lowpass and highpass filters the only available control parameter is the filter cutoff. Thus the only option which we had for compensating the bilinear transform’s frequency detuning was correcting the cutoff value. Other filters may have more parameters available. Usually their parameters (e.g. resonance) will have a strong “vertical” effect on the amplitude response and thus are not very suitable for compensating the frequency detuning. However in some cases there will be further options of horizonally altering the filter’s amplitude (and phase) responses without causing noticeable changes in the vertical direction. In such cases we will have further options for more detailed compensation of the frequency detuning.

Note, however, that these compensations, being not expressible as cutoff multiplication, may destroy the frequency response of a system of filters, unless there is some other way to make them have identical effect on all of the filters in the system. We will discuss some examples of this later in the book.

3.9 Zero-delay feedback

There is a further problem with the trapezoidal integrator replacement in the TPT method. Replacing the integrators with trapezoidal ones introduces delayless feedback loops (that is, feedback loops not containing any delay elements) into the structure. E.g. consider the structure in Fig. 3.12. Carefully examining this structure, we find that it has a feedback loop which doesn’t contain any unit delay elements. This loop goes from the leftmost summator through the gain, through the upper path of the integrator to the filter’s output and back through the large feedback path to the leftmost summator.

Why is this delayless loop a problem? Let’s consider for example the naive lowpass filter structure in Fig. 3.5. Suppose we don’t have the respective program code representation and wish to obtain it from the block diagram. We could do it in the following way. Consider Fig. 3.27, which is the same as Fig. 3.5, except that it labels all signal points. At the beginning of the computation of a new sample the signals $A$ and $B$ are already known. $A=x[n]$ is the current input sample and $B$ is taken from the internal state memory of the $z^{-1}$ element. Therefore we can compute $C=A-B$. Then we can compute $D=({\omega }_{c}T)C$ and finally $E=D+B$. The value of $E$ is then stored into the internal memory of the $z^{-1}$ element (for the next sample computation) and is also sent to the output as the new $y[n]$ value. Easy, right?

Figure 3.27: Naive 1-pole lowpass filter and the respective signal computation order.

Now the same approach doesn’t work for the structure in Fig. 3.12. Because there is a delayless loop, we can’t find a starting point for the computation within that loop.

The classical way of solving this problem is exactly the same as what we had in the naive approach: introduce a $z^{-1}$ into the delayless feedback, turning it into a feedback containing a unit delay (Fig. 3.28). Now there are no delayless feedback paths and we can arrange the computation order in a way similar to Fig. 3.27. This however destroys the resulting frequency response, because the transfer function is now different. In fact the obtained result is not significantly better (if better at all) than the one from the naive approach. There are some serious artifacts in the frequency response closer to the Nyquist frequency, if the filter cutoff is sufficiently high.

Figure 3.28: Digital 1-pole lowpass filter with a trapezoidal integrator and an extra delay in the feedback.

Therefore we shouldn’t introduce any modifications into the structure and solve the zero-delay feedback problem instead. The term “zero-delay feedback” originates from the fact that we avoid introducing a one-sample delay into the feedback (like in Fig. 3.28) and instead keep the feedback delay equal to zero.

So, let’s solve the zero-delay feedback problem for the structure in Fig. 3.12. Notice that this structure simply consists of a negative feedback loop around a trapezoidal integrator, where the trapezoidal integrator structure is exactly the one from Fig. 3.11. We will now introduce the concept of the instantaneous response of this integrator structure.

So, consider the integrator structure in Fig. 3.11. Since there are no delayless loops in the integrator, it’s not difficult to obtain the following expression for $y[n]$:

$$y[n]=\frac{{\omega }_{c}T}{2}x[n]+u[n-1]\text{\hspace{1em}(3.26)}$$

Notice that, at the time $x[n]$ arrives at the integrator’s input, all values in the right-hand side of (3.26) are known (no unknown variables). Introducing notation

$$g=\frac{{\omega }_{c}T}{2}$$

$$s[n]=u[n-1]$$

we have

$$y[n]=gx[n]+s[n]$$

or, dropping the discrete time argument notation for simplicity,

$$y=gx+s$$

That is, at any given time moment $n$, the output of the integrator $y$ is a linear function of its input $x$, where the values of the parameters of this linear function are known. The $g$ parameter doesn’t depend on the internal state of the integrator, while the $s$ parameter does depend on the internal state of the integrator. We will refer to the linear function $f(x)=gx+s$ as the instantaneous response of the integrator at the respective implied time moment $n$. The coefficient $g$ can be referred to as the instantaneous response gain or simply instantaneous gain. The term $s$ can be referred to as the instantaneous response offset or simply instantaneous offset.

Let’s now redraw the filter structure in Fig. 3.12 as in Fig. 3.29. We have changed the notation from $x$ to $\xi$ in the $gx+s$ expression to avoid the confusion with the input signal $x[n]$ of the entire filter.

Figure 3.29: 1-pole TPT lowpass filter with the integrator in the instantaneous response form.

Now we can easily write and solve the zero-delay feedback equation. Indeed, suppose we already know the filter output $y[n]$. Then the output signal of the feedback summator is $x[n]-y[n]$ and the output of the integrator is respectively $g(x[n]-y[n])+s$. Thus

$$y[n]=g(x[n]-y[n])+s$$

or, dropping the time argument notation for simplicity,

$$y=g(x-y)+s\text{\hspace{1em}(3.27)}$$

The equation (3.27) is the zero-delay feedback equation for the filter in Fig. 3.29 (or, for that matter, in Fig. 3.12). Solving this equation, we obtain

$$y(1+g)=gx+s$$

and respectively

$$y=\frac{gx+s}{1+g}\text{\hspace{1em}(3.28)}$$

Having found $y$ (that is, having predicted the output $y[n]$), we can then proceed with computing the other signals in the structure in Fig. 3.12, beginning with the output of the leftmost summator.¹⁴

It’s worth mentioning that (3.28) can be used to obtain the instantaneous response of the entire filter from Fig. 3.12. Indeed, rewriting (3.28) as

$$y=\frac{g}{1+g}x+\frac{s}{1+g}$$

and introducing notations

$$G=\frac{g}{1+g}$$

$$S=\frac{s}{1+g}$$

we have

$$y=Gx+S\text{\hspace{1em}(3.29)}$$

So, the instantaneous response of the entire lowpass filter in Fig. 3.12 is again a linear function of the input. We could use the expression (3.29) e.g. to solve the zero-delay feedback problem for some larger feedback loop containing a 1-pole lowpass filter.

3.10 Implementations

1-pole lowpass

We are now going to convert the structure in Fig. 3.12 into a piece of code. Let’s introduce helper variables into Fig. 3.12 as shown in Fig. 3.30, where we have used the already known to us fact that, given the integrator’s instantaneous response $gx+s$, the value of $s$ equals the output of the $z^{-1}$ element.

Figure 3.30: 1-pole TPT lowpass filter with helper variables.

We already know $y$ from (3.29). Since $g={\omega }_{c}T/2$, we have

$$\begin{array}{rl}v=g(x-y)& =g\left(x-Gx-S\right)=g\left(x-\frac{g}{1+g}x-\frac{s}{1+g}\right)=\\ & =g\left(\frac{1}{1+g}x-\frac{s}{1+g}\right)=g\frac{x-s}{1+g}\end{array}\text{\hspace{1em}(3.30)}$$

Now (3.30) gives a direct expression for $v$ in terms of known signals, not using $y$. In order to avoid unnecessary computations, we can simply reobtain $y$ using the obvious from Fig. 3.30 fact that $y$ is just a sum of $v$ and $s$:

$$y=v+s\text{\hspace{1em}(3.31)}$$

We also need the $z^{-1}$ input (which we need to store in the $z^{-1}$’s memory) which is also obtained from Fig. 3.30 in an obvious way:

$$u=y+v\text{\hspace{1em}(3.32)}$$

The equations (3.30), (3.31) and (3.32) can be directly expressed in program code:

// perform one sample tick of the lowpass filter
// G = g/(1+g)
// the variable 's' contains the state of z^-1 block
v := (x-s)*G;
y := v + s;
s := y + v;

or instead expressed in a block diagram form (Fig. 3.31). Notice that the block diagram doesn’t contain any delayless loops anymore.

Figure 3.31: 1-pole TPT lowpass filter with resolved zero-delay feedback.

1-pole multimode

The highpass signal can be obtained from the structure in Fig. 3.31 in a trivial manner, since $y_{\text{HP}}=x-y_{\text{LP}}$, thereby turning Fig. 3.31 into a multimode 1-pole.

1-pole highpass

If we need only the highpass signal, we could do it in a smaller number of operations than in the multimode 1-pole. By noticing that

$$y_{\text{HP}}=x-y_{\text{LP}}=x-(v+s)=(x-s)-v=(x-s)-\frac{g}{1+g}(x-s)=\frac{1}{1+g}(x-s)$$

and that

$$u=y_{\text{LP}}+v=s+2v=s+\frac{2g}{1+g}$$

we can implement the highpass filter as follows:

// perform one sample tick of the highpass filter
// Ghp = 1/(1+g)
// the variable 's' contains the state of z^-1 block
xs := x - s;
y := xs*Ghp;
s := s + y*2g;

however this way we have traded an addition/subtraction pair for one multiplication, plus instead of one cutoff-dependent parameter $G$ we need to store and access Ghp and 2g. Therefore this is not necessarily a performance improvement.

1-pole allpass

The allpass signal can be obtained from the multimode 1-pole in a trivial manner, recalling that $y_{\text{AP}}=y_{\text{LP}}-y_{\text{HP}}$. However, if we need only the allpass signal, we could save a couple of operations. Noticing that

$$y_{\text{AP}}=y_{\text{LP}}-y_{\text{HP}}=v+s-(x-(v+s))=(s+2v)-(x-s)$$

and that $s+2v=u$ is the new state of the $z^{-1}$ block, we obtain

// perform one sample tick of the allpass filter
// 2Glp=2g/(1+g)
// the variable 's' contains the state of z^-1 block
xs := x - s;
s := s + xs*2Glp;
y := s - xs;

3.11 Direct forms

Consider again the equation (3.5), which describes the application of the bilinear transform to convert an analog transfer function to a digital one. There is a classical method of digital filter design which is based directly on this transformation, without using any integrator replacement techniques. In the author’s experience, for music DSP needs this method typically has a largely inferior quality, compared to the TPT. Nevertheless we will describe it here for completeness and for a couple of other reasons. Firstly, it would be nice to try to analyse and understand the reasons for the problems of this method. Secondly, this method could be useful once in a while. Particularly, its deficiencies mostly disappear in the time-invariant (unmodulated or sufficiently slowly modulated) case.

Having obtained a digital transfer function from (3.5), we could observe, that, since the original analog transfer function was a rational function of $s$, the resulting digital transfer function will necessarily be a rational function of $z$. E.g. using the familiar 1-pole lowpass transfer function

$$H_a(s)=\frac{{\omega }_c}{s+{\omega }_c}$$

we obtain

$$\begin{array}{rl}{H}_d(z)=H_a\left(\frac{2}{T}\cdot \frac{z-1}{z+1}\right)& =\frac{{\omega }_c}{\frac{2}{T}\cdot \frac{z-1}{z+1}+{\omega }_c}=\\ & =\frac{\frac{{\omega }_{c}T}{2}(z+1)}{(z-1)+\frac{{\omega }_{c}T}{2}(z+1)}=\frac{\frac{{\omega }_{c}T}{2}(z+1)}{\left(1+\frac{{\omega }_{c}T}{2}\right)z-\left(1-\frac{{\omega }_{c}T}{2}\right)}\end{array}$$

Now, there are standard discrete-time structures allowing an implementation of any given nonstrictly proper rational transfer function. It is easier to use these structures, if the transfer function is expressed as a rational function of $z^{-1}$ rather than the one of $z$. In our particular example, we can multiply the numerator and the denominator by $z^{-1}$, obtaining

$$H_d(z)=\frac{\frac{{\omega }_{c}T}{2}(1+z^{-1})}{\left(1+\frac{{\omega }_{c}T}{2}\right)-\left(1-\frac{{\omega }_{c}T}{2}\right)z^{-1}}$$

The further requirement is to have the constant term in the denominator equal to 1, which can be achieved by dividing everything by $1+{\omega }_{c}T/2$:

$$H_d(z)=\frac{\frac{\frac{{\omega }_{c}T}{2}}{1+\frac{{\omega }_{c}T}{2}}(1+z^{-1})}{1-\frac{1-\frac{{\omega }_{c}T}{2}}{1+\frac{{\omega }_{c}T}{2}}{z}^{-1}}\text{\hspace{1em}(3.33)}$$

Now suppose we have an arbitrary rational nonstrictly proper transfer function of $z$, expressed via $z^{-1}$ and having the constant term in the denominator equal to 1:

$$H(z)=\frac{\sum _{n=0}^N{b}_n{z}^{-n}}{1-\sum _{n=1}^N{a}_n{z}^{-n}}\text{\hspace{1em}(3.34)}$$

This transfer function can be implemented by the structure in Fig. 3.32 or by the structure in Fig. 3.33. One can verify (by computing the transfer functions of the respective structures) that they indeed implement the transfer function (3.34). There are also transposed versions of these structures, which the readers should be able to construct on their own.

Let’s use the direct form II to implement (3.33). Apparently, we have

$$N=1$$

$$b_0=b_1=\frac{\frac{{\omega }_{c}T}{2}}{1+\frac{{\omega }_{c}T}{2}}$$

$$a_1=\frac{1-\frac{{\omega }_{c}T}{2}}{1+\frac{{\omega }_{c}T}{2}}$$

and the direct form implementation itself is the one in Fig. 3.34 (we have merged the $b_0$ and $b_1$ coefficients into a single gain element).

In the time-invariant (unmodulated) case the performance of the direct form filter in Fig. 3.34 should be identical to the TPT filter in Fig. 3.12, since both implement the same bilinear-transformed analog transfer function (2.5). When the cutoff is modulated, however, the performance will be different.

Figure 3.32: Direct form I (DF1).

Figure 3.33: Direct form II (DF2), a.k.a. canonical form.

Figure 3.34: Direct form II 1-pole lowpass filter.

We have already discussed in Sections 2.7 and 2.16 that different topologies may have different time-varying behavior even if they share the same transfer function. Apparently, the difference in behavior between Fig. 3.34 and Fig. 3.12 is another example of that. Comparing the implementations in Figs. 3.34 and 3.12, we notice that the structure in Fig. 3.34 contains a gain element at the output, the value of this gain being approximately proportional to the cutoff (at low cutoffs). This will particularly produce unsmoothed jumps in the output in response to jumps in the cutoff value. In the structure in Fig. 3.12, on the other hand, the cutoff jumps will be smoothed by the integrator. Thus, the difference between the two structures is similar to the just discussed effect of the cutoff gain placement with the integrator.

We should conclude that, other things being equal, the structure in Fig. 3.34 is inferior to the one in Fig. 3.12 (or Fig. 3.31). In this respect consider that Fig. 3.12 is trying to explicitly emulate the analog integration behavior, preserving the topology of the original analog structure, while Fig. 3.34 is concerned solely with implementing a correct transfer function. Since Fig. 3.34 implements a classical approach to the bilinear transform application for digital filter design (which ignores the filter topology) we’ll refer to the trapezoidal integration replacement technique as the topology-preserving bilinear transform (or, shortly, TPBLT). Or, even shorter, we can refer to this technique as simply the topology-preserving transform (TPT), implicitly assuming that the bilinear transform is being used.¹⁵

In principle, sometimes there are possibilities to “manually fix” the structures such as in Fig. 3.34. E.g. the time-varying performance of the latter is drastically improved by moving the $b_0$ gain to the input. The problem however is that this kind of fixing quickly gets more complicated (if being possible at all) with larger filter structures. On the other hand, the TPT method explicitly aims at emulating the time-varying behavior of the analog prototype structure, which aspect is completely ignored by the classical transform approach. Besides, if the structure contains nonlinearities, preserving the topology becomes absolutely critical, because otherwise these nonlinearites can not be placed in the digital model.¹⁶ Also, the direct forms suffer from precision loss issues, the problem growing bigger with the order of the system. For that reason in practice the direct forms of orders higher than 2 are rarely used,¹⁷ but even 2nd-order direct forms could already noticeably suffer from precision losses.

3.12 Transient response

Looking at the 1-pole lowpass filter’s discrete-time transfer function (3.33) and noticing that ${\omega }_c=-p$ where $p$ is the analog pole, we could rewrite (3.33) as

$$H(z)=\frac{\frac{\frac{-pT}{2}}{1-\frac{pT}{2}}\left(1+z^{-1}\right)}{1-\frac{1+\frac{pT}{2}}{1-\frac{pT}{2}}{z}^{-1}}$$

Comparing this to (3.9) we notice that

$$\frac{1+\frac{pT}{2}}{1-\frac{pT}{2}}=\tilde{p}$$

where $\tilde{p}$ is the result of the application of the inverse bilinear transform formula (3.9) to $p$. Further noticing that

$$\frac{\frac{-pT}{2}}{1-\frac{pT}{2}}=\frac{1}{2}\cdot \left(\frac{1-\frac{pT}{2}}{1-\frac{pT}{2}}-\frac{1+\frac{pT}{2}}{1-\frac{pT}{2}}\right)=\frac{1-\tilde{p}}{2}\text{\hspace{1em}(3.35)}$$

we rewrite $H(z)$ as

$$H(z)=\frac{1-\tilde{p}}{2}\cdot \frac{1+z^{-1}}{1-\tilde{p}{z}^{-1}}=\frac{1-\tilde{p}}{2}\cdot \frac{z+1}{z-\tilde{p}}$$

Thus $\tilde{p}$ is the pole of $H(z)$, as we should have expected.

On the other hand, applying trapezoidal integration to the 1-pole lowpass differential equation in the pole form (2.14), we have

$$y[n]-y[n-1]=pT\cdot \left(\frac{y[n]+y[n-1]}{2}-\frac{x[n]+x[n-1]}{2}\right)$$

from where

$$\left(1-\frac{pT}{2}\right)y[n]=\left(1+\frac{pT}{2}\right)y[n-1]-\frac{pT}{2}\cdot (x[n]+x[n-1])$$

$$\begin{array}{rl}y[n]& =\frac{1+\frac{pT}{2}}{1-\frac{pT}{2}}y[n-1]+\frac{\frac{-pT}{2}}{1-\frac{pT}{2}}\cdot (x[n]+x[n-1])=\\ & =\tilde{p}y[n-1]+(1-\tilde{p})\frac{x[n]+x[n-1]}{2}\end{array}$$

where we have used (3.35).

Now consider a complex exponential $x[n]=X(z)z^n$. For such $x[n]$ we have

$$y[n]=\tilde{p}y[n-1]+(1-\tilde{p})\frac{1+z^{-1}}{2}X(z)z^n=\tilde{p}y[n-1]+\tilde{q}{z}^n\text{\hspace{1em}(3.36)}$$

where we introduced notation

$$\tilde{q}=(1-\tilde{p})\frac{1+z^{-1}}{2}X(z)$$

for conciseness. Recursively substituting (3.36) into itself at progessively decreasing values of $n$ we obtain

$$\begin{array}{rl}y[n]& =\tilde{p}y[n-1]+\tilde{q}{z}^n=\\ & =\tilde{p}\left(\tilde{p}y[n-2]+\tilde{q}{z}^{n-1}\right)+\tilde{q}{z}^n=\\ & ={\tilde{p}}^{2}y[n-2]+\left(\tilde{p}{z}^{-1}+1\right)\tilde{q}{z}^n=\\ & ={\tilde{p}}^2\left(\tilde{p}y[n-3]+\tilde{q}{z}^{n-2}\right)+\left(\tilde{p}{z}^{-1}+1\right)\tilde{q}{z}^n=\\ & ={\tilde{p}}^{3}y[n-3]+\left({\left(\tilde{p}{z}^{-1}\right)}^2+\tilde{p}{z}^{-1}+1\right)\tilde{q}{z}^n=\\ & \dots \\ & ={\tilde{p}}^{n}y[0]+\left({\left(\tilde{p}{z}^{-1}\right)}^{n-1}+{\left(\tilde{p}{z}^{-1}\right)}^{n-2}+\dots +\tilde{p}{z}^{-1}+1\right)\tilde{q}{z}^n=\\ & ={\tilde{p}}^{n}y[0]+\frac{{\left(\tilde{p}{z}^{-1}\right)}^n-1}{\tilde{p}{z}^{-1}-1}\tilde{q}{z}^n={\tilde{p}}^{n}y[0]+\frac{{\tilde{p}}^n-z^n}{\tilde{p}-z}\tilde{q}z=\\ & ={\tilde{p}}^{n}y[0]+\frac{z^n-{\tilde{p}}^n}{z-\tilde{p}}(1-\tilde{p})\frac{z+1}{2}X(z)=\\ & ={\tilde{p}}^{n}y[0]+(z^n-{\tilde{p}}^n)H(z)X(z)=\\ & =H(z)X(z)z^n+(y[0]-H(z)X(z))\cdot {\tilde{p}}^n=\\ & =y_s[n]+(y[0]-y_s[0])\cdot {\tilde{p}}^n=y_s[n]+y_t[n]\end{array}$$