Tik-61.246 Digital Signal Processing and Filtering (DISKO)
Solutions, exercise round 6.

1.

Normalizing and scaling of H(z) is done because most design algorithms and software (Matlab) require normalized values.

$$H(z) = K \frac{1+z^{-1}+z^{-2}+z^{-3}+z^{-4}}{z^{-4}}$$

(H(z) is an anti-causal FIR, $H(z) = K (z^4+z^3+z^2+z+1) \leftrightarrow y[n]=K(x[n+4]+x[n+3]+x[n+2]+x[n+1]+x[n])$.

H(z) is a lowpass filter (anti-causal moving average filter). Maximum value of lowpass filter is reached at $\omega = 0$. Because $z=e^{j\omega}$, we get z=1 (see figure 1).

$H(1) = K \frac{1+1+1+1+1}{1} = 1$

$\Rightarrow K = \frac15$

  

Figure 1: Note that scaling constant doesn't affect locations of poles and zeros

[Pole-zero diagram] [Frequency response of H(z) with K=0.2]

2.

Impulse-invariant transform (Mitra 7.2, p.425-431) and bilinear transform (Mitra 7.3, p. 432-435) for IIR filters. First, find an analog filter which satisfies the specifications. Then, convert it into a digital filter either by impulse-invariant method or bilinear transform.

In the impulse-invariant method the target is to get impulse response of digital filter h[n] to be the same as the sampled impulse response of analog filter h_a(nT). Because IIR filters have normally an impulse response of infinite length, this method brings distortion.

In the bilinear transform the whole left subspace of s-plane is mapped into unit circle of z-plane. This is done by substituting variables (see 2b).

The definition of the first order (N=1) analog Butterworth filter is (Sec. 5.3.2, Eq. 5.29)

$$\vert H_a(j\Omega)\vert^2=\frac{1}{1+\left(\frac{\Omega}{\Omega_c}\right)^2}$$

$$H_a(s)H_a(-s)=\frac{1}{1+\left(\frac{s}{j\Omega_c}\right)^2} =$$

$$\frac{1}{1-\left(\frac{s}{\Omega_c}\right)^2} = \overbrace{\f... ...verbrace{\frac{1}{1+\left(\frac{-s}{\Omega_c}\right)}}^{H(-s)}$$

where $s=j\Omega$

$$\Rightarrow H_a(s)=\frac{\Omega_c}{s+\Omega_c}$$

The pole is

$$s+\Omega_c = 0 \Rightarrow s= - \Omega_c$$

$\Omega$ refers frequency in analog domain ( $H(j\Omega)$) and $\omega$ frequency in digital domain ( $H(e^{j\omega})$).

a)

Impulse-invariant transform is acquired when (Eqs. 7.29 - 7.32)

$$H_a(s)=\frac{A}{s+s_k} \mapsto h_a(t)=Ae^{-s_k t} \mu(t) \mapsto$$

$$h[n] = h_a(nT) = Ae^{-s_k nT} \mu[n] \mapsto H(z)= \frac{A}{1-e^{-s_k T}z^{-1}}$$

The filter $\frac{\Omega_c}{s+\Omega_c}$ is already in that form ( $A_k = s_k = \Omega_c$), and thus we acquire (K is a constant to be set later so that maximum of $\vert H(e^{j\omega})\vert$ is 1)

$$H(z)=\frac{K}{1-e^{-\omega_c T}z^{-1}}.$$

Keeping in mind that $\omega_c=2\pi f_c$ and $T=\frac{1}{f_s}$, we get

$$H(z)=\frac{K}{1-e^{-2\pi f_c/f_s}z^{-1}}$$

By inserting the values of f_s=1 kHz (sampling frequency) and f_c=100 Hz (cut-off frequency), we get

$$H(z)=\frac{K}{1-e^{-\pi/5}z^{-1}}$$

The transfer function is then scaled by a factor K so that the maximum of the magnitude response will be one. We also know that the maximum is located at zero frequency, because the frequency response of a Butterworth filter is monotonic. Thus we get

$$\frac{K}{1-e^{-\pi/5}}=1 \Leftrightarrow K=1-e^{-\pi/5}$$

The amplitude (power) response of the filter is therefore (figure 2)

$$H(z)=\frac{1-e^{-\pi/5}}{1-e^{-\pi/5}z^{-1}}$$

  

Figure: Power spectrum $\vert H(e^{j\omega })\vert^2$ by Impulse-invariant method. Cut-off frequency (-3 dB) is at 100 Hz / (1000 Hz/ 2) = 0.2

[linear scale 0..1] [desibels]

b)

A bilinear transform is acquired when (Eqs. 7.46)

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

is inserted into system function (Eqs. 7.47)

$$H(z) = H_a(s) \vert _{s= \frac2T (\frac{1-z^{-1}}{1+z^{-1}})}$$

For an analog filter, the cutoff frequency must then be altered to

$$\Omega_c=\frac{2}{T}\tan\left(\frac{\omega_cT}{2}\right)$$

where $\omega_c = 2\pi \, 100$ rad/s is the desired cutoff angular frequency. The warped angular frequency $\Omega_c = 2\pi \, 103.4$ rad/s is used in analog filter (see Figs. 7.7, 7.8 in Mitra). Along with the scaling coefficient K, the transfer function of the filter will be

$$H(z) = H_a(s) \vert _{s= \frac2T (\frac{1-z^{-1}}{1+z^{-1}})}... ...\Omega_c}^{=K'} (1 + z^{-1})}{2(1-z^{-1})+T\Omega_c(1+z^{-1})}$$

All constants (related to scaling of filter to unity) can be combined, and constants grouped in function of z

$$H(z)=K'\frac{1+z^{-1}}{(2+\Omega_c T)+(\Omega_c T-2)z^{-1}}$$

Because $(2+\Omega_cT)$ is constant, we can write

$$H(z)= K'' \frac{1+z^{-1}}{1+\frac{\Omega_cT-2}{\Omega_cT+2}z^{-1}}$$

Now the transfer function is in standart form. By inserting $\Omega_c=\frac{2}{T}\tan\frac{\omega_cT}{2}=\frac{2}{T}\tan\frac{\pi}{10}$, we get

$$H(z)= K''\frac{1+z^{-1}}{1+\frac{\tan(\pi/10)-1}{\tan(\pi/10)+1}z^{-1}}.$$

As earlier, the maximum of the magnitude response is at zero frequency

$$\Rightarrow K''\frac{1+1}{1+\frac{\tan(\pi/10)-1}{\tan(\pi/10... ...ac{1}{2}\left(1+\frac{\tan(\pi/10)-1}{\tan(\pi/10)+1}\right).$$

Thus, the frequency response will be (figure 3)

$$H(z)=\frac{1}{2}\left(1+\frac{\tan(\pi/10)-1}{\tan(\pi/10)+1}... ...\frac{1+z^{-1}}{1+\frac{\tan(\pi/10)-1}{\tan(\pi/10)+1}z^{-1}}$$

  

Figure: Power spectrum $\vert H(e^{j\omega })\vert^2$ by Bilinear transform. Cut-off frequency (-3 dB) is at 100 Hz / (1000 Hz/ 2) = 0.2

[linear scale 0..1] [desibels]

3.

Analog lowpass filter.

a)

We are examining a continuous-time lowpass filter (figure 4)

$$H_a(s)=\frac{1}{s+1}$$

$$\begin{array}{ll} {\rm zero :} & s=\infty\\ {\rm pole :} & s+1=0 \Rightarrow s=-1 \end{array}$$

  

Figure 4: The original LP filter.

b)

In the first case, we alter the cut-off frequency of a lowpass filter (see the figure 5).

$$H'(s)=H_c(s/\Omega_c)=\frac{1}{\frac{s}{\Omega_c} +1}=\frac{\Omega_c}{s+\Omega_c}$$

$$\begin{array}{ll} {\rm zero :} & s=\infty\\ {\rm pole :} & s=-\Omega_c \end{array}$$

  

Figure: LP $\rightarrow$ LP, altering cutoff frequency

c)

In the second case, we transform a lowpass filter to a highpass filter (see the figure 6).

$$H''(s)=H_c(\Omega_c/s)=\frac{1}{\frac{\Omega_c}{s} + 1}= \frac{s}{s+\Omega_c}$$

$$\begin{array}{ll} {\rm zero :} & s=0\\ {\rm pole :} & s=-\Omega_c \end{array}$$

  

Figure: LP $\rightarrow$ HP

d)

Bilinear transform (T=2)

$$s= \frac{1-z^{-1}}{1+z^{-1}} \Leftrightarrow z=\frac{1+s}{1-s}$$

is a mapping from Laplace plane s to z-plane (figure 7) so that

$$\Omega=0 ~\mapsto~ z=1$$

$$\Omega=\infty ~\mapsto~ z=-1.$$

  

Figure 7: Bilinear mapping

The zero of $H_a(s)=\frac{1}{s+1}$ is $s=\infty$ and pole s=-1.

The digital filters acquired by applying bilinear transform :

i)

zero:

$$s=\infty \Rightarrow z=\frac{1+s}{1-s}=\frac{\frac{1}{s}+1}{\frac{1}{s}-1}=-1$$

pole (omega_c will not be the same as $\Omega_c$, use warping)

$$s=-\Omega_c \Rightarrow z=\frac{1-\frac{T}{2}\omega_c}{1+\frac{T}{2}\omega_c} \in \Re, ~-1<z<1$$

  

Figure 8: H(z), Bilinear transform of lowpass filter

ii)

zero:

$$s=0 \Rightarrow z=\frac{1+1\cdot 0}{1-1\cdot 0}=1$$

The pole is the same as in previous case.

  

Figure 9: H(z), Bilinear transform of highpass filter

4.

FIR, truncated Fourier series. See p. 447-448 in Mitra.

Show that the best finite-length approximation to the ideal infinite length impulse response in the mean square error sense is obtained by truncated Fourier series method.

Let the $H_d(e^{j\omega})$ denote the desired frequency response function. Since $H_d(e^{j\omega})$ is a periodic function of $\omega$with the period $2\pi$, it can be expressed as a Fourier series:

$$H_d(e^{j\omega}) =\sum^\infty_{n=-\infty}h_d[n]e^{-j\omega n},$$

where the Fourier coefficients h_d[n] are the corresponding impulse-response samples and are given by

$$h_d[n]=\frac{1}{2\pi}\int^\pi_{-\pi} H_d(e^{j\omega})e^{j\omega n}d\omega,\ \ -\infty\leq n\leq\infty.$$

For most practical solutions, h_d is of infinite length and noncausal. Therefore we try to find a finite-duration impulse response sequence h_t[n] of length 2M+1, whose DTFT $H_t(e^{j\omega})$approximates the $H_d(e^{j\omega})$ in some sense. One commonly used approximation criteria is to minimize the integral squared-error

$$\Phi=\frac{1}{2\pi}\int_{-\pi}^\pi\vert H_t(e^{j\omega})-H_d(e^{j\omega})\vert^2d\omega$$

where

$$H_t(e^{j\omega})=\sum^M_{n=-M}h_t[n]e^{-j\omega n}.$$

Using Parseval's relation

$$\sum_{n=-\infty}^\infty g[n]h^*[n]=\frac{1}{2\pi}\int^\pi_{-\pi}G(e^{j\omega})H^*(e^{j\omega})d\omega,$$

we can rewrite the cost function

$$\Phi=\sum_{-\infty}^{\infty}\vert h_t[n]-h_d[n]\vert^2=\sum^M... ...+\sum^{-M-1}_{n=-\infty}h_d^2[n]+\sum_{n=M+1}^{\infty}h_d^2[n].$$

It obvious from the latter form that the integral-squared error is a minimum when h_t[n]=h_d[n], for $-M\leq n\leq M$, that is the best finite-length approximation to the ideal infinite-length impulse response in the mean-square error sense is obtained by truncation.

The causality is achieved by delaying the h_t[n] by M samples.

What is the disadvantage of this method and what are the solutions to this problem?

Disadvantage is the oscillatory behaviour of $H_t(e^{j\omega})$ (Gibbs phenomenon). This is caused by simple truncation (window function with abrupt transitions in time domain) and the instability of the ideal filter. Possible solutions are the use of tapered windows (fixed or adjustable) and specification of FIRs with smooth transitions.

Some examples of window functions:

i)

Rectangular N=11, Figure 10

ii)

Rectangular N=65, Figure 11

iii)

Hamming N=65, Figure 12

There are three figures for each item. Top left figure is the window function in time domain w[n]. The causal version can be obtained by shifting.

Bottom left figure is the window function in frequency domain $W(e^{j\omega})$.

The third figure in right is the amplitude frequency of actual filter which is obtained via window function method. The desired lowpass filter has normalized cut-off frequency at 0.2.

Notice that

i)

Rectangular N=11 gives insufficient result

ii)

Rectangular N=65 gives sharp transition band but oscillates (Gibbs phenomen)

iii)

Hamming N=65 is flat both in passband and stopband but the transition band is not as tight as in ii)

  

Figure 10: Rectangular window N=11

[ $w_1[n], W_1(e^{j\omega })$] [ $H_{t1}(e^{j\omega })=H_d(e^{j\omega })*W_1(e^{j\omega })$]

  

Figure 11: Rectangular window N=65

[ $w_2[n], W_2(e^{j\omega })$] [ $H_{t2}(e^{j\omega })=H_d(e^{j\omega })*W_2(e^{j\omega })$]

  

Figure 12: Hamming window N=65

[ $w_3[n], W_3(e^{j\omega })$] [ $H_{t3}(e^{j\omega })=H_d(e^{j\omega })*W_3(e^{j\omega })$]

5.

Systems in parallel: H(z) = H₁(z) + H₂(z)

0.7 #1#2#3#4#5 @font#1#2pt #3#4#5

Systems in cascade: H(z) = H₁(z) H₂(z)

0.7 #1#2#3#4#5 @font#1#2pt #3#4#5

The flow diagram of the system being investigated is

0.7 #1#2#3#4#5 @font#1#2pt #3#4#5

If the system being analyzed is complex, ``temporary points'' like w(n) may be used to help in the formation of the transfer function.

Using the flow diagram, we get the following equations:

$$\left\{ \begin{array}{ll} Y(z) & = F(z)[X(z)E(z) + W(z)] \\ W(z) & = G(z)Y(z) \\ \end{array} \right.$$

$$\begin{array}{ll} \Leftrightarrow & Y(z) = F(z)[X(z)E(z) + G(... ...Y(z) = X(z) \frac{F(z)E(z)}{1 - F(z)G(z)} \\ [4mm] \end{array}$$

$$\Rightarrow H(z) = \frac{Y(z)}{X(z)} = \frac{E(z)F(z)}{1 - F(z)G(z)}$$

6.

Let us use three temporary signals W₁, W₂, and W₃in the following locations:

From the above figure we get the following expressions:

W₁ = KX+z^-1W₃ (1)

$$(z^{-1}-\alpha)W_1$$ W₂ = (2)

$$\alpha W_1-\beta z^{-1}W_1 = (\alpha-\beta z^{-1})W_1$$ W₃ = (3)

$$z^{-1}W_2+\beta W_1$$ Y = (4)

Substituting the equation (3) in (1) we get

$$W_1 = KX + z^{-1}(\alpha-\beta z^{-1})W_1$$ (5)

or

$$(1-\alpha z^{-1}+\beta z^{-2})W_1 = KX$$ (6)

Next, substituting (2) in (4) we get

$$Y = [z^{-1}(z^{-1}-\alpha)+\beta]W_1$$ (7)

Then, from (6) and (7) we finally arrive at

$$H(z)=\frac{Y(z)}{X(z)} = K \frac{\beta-\alpha z^{-1}+z^{-2}}{1-\alpha z^{-1} + \beta z^{-2}}$$ (8)

a)

Since the structure employs 4 unit delays to implement a second-order transfer function, it is noncanonic.

b)

and c)

$$\begin{eqnarray*}H(z)H(z^{-1}) &=& K^2\left(\frac{\beta-\alpha z^{-1} + z^{-2}}... ...ha z^{-1} + 1} {z^{-2}-\alpha z^{-1} +\beta}\right) \\ &=& K^2 \end{eqnarray*}$$

Therefore $\vert H(e^{j\omega})\vert=K$ for all values of $\omega$ and hence $\vert H(e^{j\omega})\vert=1$ if K=1.

d)

H(z) is an allpass transfer function with a constant magnitude at all values of $\omega$.

7.

A canonic direct form II realization of H(z):

Its transposed realization:

[h]

0.6 #1#2#3#4#5 @font#1#2pt #3#4#5

The error shaping (FIR) filter has transfer function

H(z) = 1 - 2z^-1 + z^-2 = (1 - z^-1)².

(a)

The frequency response of the filter is

$$H(e^{j\omega}) = 1 - 2e^{-j \omega}+ e^{-2j \omega} = e^{-j \... ...omega} - 2 + e^{-j \omega}) = e^{-j \omega} (2 cos(\omega) -2)$$

The power spectrum of the noise has the same shape as the filter's squared amplitude response (because the noise was assumed to be white).

(b)

The variance of the noise at the output is defined by formula

$$\sigma^2_{e,tot} = \sigma^2_e \sum_j k_j \sum_{n=0}^{\infty} {\vert h^{\star}_j[n]\vert}^2.$$

In this case, there is only one error source (j=1), and
$h^{\star}[n] = Z^{-1} \left\{ 1 - 2z^{-1} + z^{-2} \right\} = \delta[n] -2 \delta[n-1] + \delta[n-2]$.

$$\Rightarrow \sigma_{e,tot}^2 = \sigma_e^2 \sum_{n = 0}^2 {\ve... ... \sigma_e^2 \left( 1^2 + {(-2)}^2 + 1^2 \right) = 6 \sigma^2_e$$

(c)

The filter reduces quantization error on the band of interest.