Digital Signal Processing Exercises With Solutions

For L 1a:e aL 1 G(f, L) ' 1 XL (f ) X(f ) G(f, L) ' X(f ) Arg [G(f, L)] ' Arctan e aL sin(2πf L) ' e aL sin(2πf L) Arg [G(f, L)] e aL ' 0 Arg [XL (f )] Arg [X(f )] Arg [G(f, L)] ' Arg [X(f )]Thus the Fourier transform is not distorted a lot when the signal is observed on a su cient time duration (true for L a4 , then e aL ' 0.02 1).! !Be careful ! !If we call wL (t) the window modeling the signal limi-tation to a L duration, we actually have T F [x(t)wL (t)] X(f ) WL (f ),where WL (f ) represents the window Fourier transform (rectangular windowhere). We just show, in this exercise, that, in that case, the Fourier transformdegradation (due to the convolution with WL (f )) is a multiplicative errorcalled G(f, L), meaning : T F [x(t)wL (t)] X(f ) W (f ) X(f )G(f, L).2.1.3 Exercise 3 : DFT of a function with continuous spectrum, signal sampling e ectX(f ) (see equation 2.1) has a unbounded support. Thus, theoreticallyx(t) cannot be sampled without loss of information. However, X(f ) tends to0 when f . Practically the signal x(t) can be sampled considering an"acceptable" loss of information. This "acceptable" loss of information in theexercise consists in considering that X(f ) is null for an attenuation higherthan 40 dB compared to its maximum value. Then a maximum frequencyFmax can be deducted and used to dimension the sampling frequency Fe .20log X(Fmax ) 20log X(0) 40dB 20log X(0) 20log102 Fmax '100a100a Fe 2ππ12

The sampled signal Fourier transform is then written as : XXe (f ) x(kTe )e j2πf kTe(2.2)k which gives, for the considered signal :Xe (f ) hikXe (a j2πf Te ) k 011 e (a j2πf )TeX(f ) is periodical, of period Fe : X(f pFe ) X(f ), p Z. Forf Fe22fFe 1 it is possible to write a limited development of theexponential function (we also have aTe Xe (f ) 100π 1) leading to :1 Fe X(f )a j2πf TeThis result can be explained by the impact of the rst order periodization,increasing when f tends to deFe2but staying low when f Fe2 .Factor Feis due to the de nition of the DFT (equation (2.2)) forgetting factor Te .Remark :do not forget to use an anti-aliasing lter before sampling x(t).2.1.4 Exercise 4 : DFT of a function with continuous spectrum, signal sampling and limitation of the signal duration e ectsThe Fourier transform of the sampled and N points limited signal canbe written as :XD (f ) N 1Xx(kTe )e j2πf kTek 0and computed for the considered signal :XD (f ) N 1 hXe (a j2πf Te )k 0ik 1 e (a j2πf )N Te1 e (a j2πf )TeHere the two previous approximations are combined.13(2.3)

2.1.5 Exercise 5 : DFT of a function with non continuousspectrum, spectrum sampling e ectLet's consider :x(t) Aej(2πf0 t φ) , t R(2.4)1. Compute the DFT X(f ) of x(t).2. x(t) is observed on a limited duration [0, L]. Compute the DFT XL (f )of this limited part of x(t) and compare it to X(f ).3. Determine the DFT, XD (f ), considering N points of the signal sampled at Te and compare it to XL (f ).4. Only a nite number of points are able to be computed to get the DFT : X(n) XD (n FNe ) for n 0, ., N 1. Determine X(n) andplot X(n) , for the fwo following cases : when f0 when f0 n0N Fe , where n0 n0 N Fe , where n0Z, Z and 0 1Conclusion ?2.1.6 Exercise 6 : DFT of a function with non continuousspectrum, spectrum sampling e ectIn order to better appreciate the impact of the last approximation (onlyN points are computed for the DFT) we will work with a signal whosespectrum is discontinuous.1. x(t) Aejφ ej2πf0 t X(f ) Aejφ δ (f f0 )2. We have two ways to calculate XL (f ) :ZXL (f ) Lx(t)ej2πf t dt ALejφ sinc (π(f f0 )L) e jπ(f f0 )L0or LXL (f ) T F x(t)ΠL t X(f ) {Lsinc(πf L)e jπf L }2 ALejφ sinc (π(f f0 )L) e jπ(f f0 )LΠL t L2 (2.5)represents the rectangular truncation window of length Lcentered around t L2.14

3. The Fourier transform considering N points of the signal sampled atTe can be written as :XD (f ) AejφN 1 hXe jπ(f f0 )Te )ik k 0 Aejφ e jπ(f f0 )(N 1)Te1 e j2π(f f0 )N Te1 e j2π(f f0 )Tesin (π(f f0 )N Te )sin (π(f f0 )Te )(2.6)4. Only a nite number of points (N ) are able to be computed to getthe DFT. These N points are computed on a period Fe , leading to a computation step FNe : X(n) XD (n FNe ) for n 0, ., N 1. Thefrequency variable f becomes n FNe , with n 0, ., N 1. In the casewhere f0 n0 FNe , we obtain :XD (n) Aejφ e jπ(n n0 )(N 1)Nsin (π(n n0 )) 0)sin π (n nNgiving XD (n) Aejφ N for n n0 and 0 elsewhere. Figure 2.4 givesan example for XD (n) when n0 4. Note that this represents thetheoretical Fourier transform for an exponential with frequency isf0 : a Dirac at f0 . However, the case where f0 n0 FNe is a very particular case. Most of the time f0 (n0 ε) FNe . An example of thecorresponding XD (n) is plotted on gure 2.5. With this plot it is notobvious to deduct that this Fourier transform corresponds to an exponential function. The solution to get a better spectrum visualizationis to make an interpolation. This interpolation is usually done usingthe Zero Padding technic.Figure 2.4 Exemple de représentation de XD (n) dans le cas où n0 4.15

Figure 2.5 Exemple de représentation de XD (n) dans le cas où n0 4.2.2 z-Transform2.2.1 Exercise 11.x(n) an u(n), a 1 X(z) X az 1 n n 011 az 1Domain of existence :limn pn an z n 1 f or z a 2.y(n) an u( n 1), a 1( ) 1 XX 1 n X 1 n n Y (z) az a z a 1 z 1 n n 1n 0Domain of existence :limn pn a n z n 1 f or z a Conclusion : the inverse z-transform is not unique. Its expression depends on the closed curve chosen to compute it. If this curve enclosesall the poles a causal solution is obtained, else not.3. Several methods are possible to solve this question. For example, apartial fraction decomposition can be used :H(z) 1a b ab 1 az 1 1 bz 116 11 az 1

for z a :h(n) 1 n 1 au( n 1) bn 1 u( n 1)a bNo poles are included in the considered closed curve. h(n) is noncausal. pour a z b :h(n) 1 n 1au(n) bn 1 u( n 1)a bThe considered closed curved encloses z a and excludes z b.h(n) has a causal part and a non causal part. for z b :h(n) 1 n 1a bn 1 u(n)a bThe considered closed curved encloses all the poles (z a etz b). The inverse z-transform h(n) is causal.2.2.2 Exercice 21.y(n) ay(n 1) bn u(n) Y (z) az 1 Y (z) Y (z) 11 az 1 11 bz 111 bz 1 If we want to obtain a causal system, then :y(n) 1 n 1a bn 1 u(n)a b2.y(n) ay(n 1) x(n) Y (z) az 1 Y (z) X(z) H(z) h(n) T Z 1 [H(z)] an u(n) (causal f ilter).17Y (z)1 X(z)1 az 1

2.3 Digital ltering2.3.1 Exercise 1H(z) 1(1 az 1 ) (1 bz 1 )(2.7)1. The lter de ned by H(z) is a 2 order lter (degree of the denominator).2. Y (z) 1 (a b)z 1 abz 2 X(z)T Z 1 y(n) x(n) (a b)y(n 1) aby(n 2)3. This lter is a RII lter because it includes a feedback loop : theoutput at time instant n depends on the input at time instant n butalso on the output past values.4. This lter will be stable if the poles of H(z) have a modulus lowerthan 1. For that we need a 1 and b 1, which is the case here.Thus, the lter is stable. n 115. h(n) a ba bn 1 u(n) leads to a causal lter.2.3.2 Exercice 21.Y (z) X(z) az 1 X(z) H(z) 2.δ(n) (z) XY (z) 1 az 1X(z)δ(n)z n δ(0)z 0 1n T Z [δ(n 1)] z 1 T Z [δ(n)] z 1So :H(z) Y (z) 1 az 1 h(n) T Z 1 [H(z)] δ(n) aδ(n 1)X(z)18

Note : it was possible to directly obtain h(n) by taking x(n) δ(n)in the recursive equation de ning the system.3.U (z) Xn 0z n 11 z 1Domain of existence :limn pn z n 1 f or z 1Step response :Y (z) H(z)U (z) 1 az 11az 1 y(n) u(n) au(n 1)1 z 11 z 1 1 z 1Note : it was possible to directly obtain the step response by takingx(n) u(n) in the recursive equation de ning the system.4. This lter is a nite impulse response (FIR) lter because there is'ntany feedback loop : the output at time instant n does not depend onits past values.5. A FIR lter is stable if its coe cients are nite, which is the casehere.6. This lter is causal because its impulse response is null for n 0.2.3.3 Exercise 3 : FIR low pass lter synthesisThe ideal frequency response HIP B (f ) is periodical, thus can be decomposed in Fourier series :HIP B (f ) XhIP B (k) e j2πf k(2.8)ek where the Fourier series coe cients hIP B (k) represent the elements (or pointsor coe cients) of the lter impulse response. Here :ZhIP B (k) eHIP B (fe)e j2πf k dfe 19Zfec fecee j2πf k dfe(2.9)

leading to, after calculations, hIP B (k) 2fec sinc(2π fec k). Practically, thenumber of coe cients de ning the lter must be limited to a given number, 2N 1, called lter order. This limitation is modeled by the use of atruncation (or weighted) window, w(k), of length 2N 1 :hP B (k) hIP B (k) w(k)It leads to an approximated frequency response :HP B (fe) HIP B (fe) W (fe)where W (fe) is the Fourier transform of w(k).Two truncation windows are proposed in this exercise : rectangular or Hamming. Using the rectangular window it comes :hP B (k) 2fec sinc(2π fec k) f or k N, ., N 0 elsewhere(2.10)Using the Hamming window we get : 2πkeehP B (k) 2fc sinc(2π fc k) 0.54 0.46 cosf or k N, ., N2N 1 0 elsewhere(2.11)Figures 2.6 and 2.7 plot the impulse responses and transfer functionsof the two obtained lowpass lters for a cut o frequency fc 100 Hz, asampling frequency Fe 800 Hz and an order N 31.Figure 2.8 plots an example of input and output signals using the lterwith rectangular truncation window. A delay between the output and theinput signal can be observed. It is due to the shifted impulse response (to geta causal lter). The output signal is distorted compared to the input one.Indeed some frequencies have been removed by the lter, as it is shown ingure 2.9.20

Figure 2.6 Low pass lter impulse responses for fc 100 Hz, Fe 800Hz and N 31.Figure 2.7 Low pass lter transfer functions for fc 100 Hz, Fe 800Hz and N 31.21

Figure2.8 Example of input and output signals using a FIR lowpasslter synthesized with the following parameters : fc 100 Hz, Fe 800 Hz,order 31, rectangular truncation window.Figure 2.9 DFT of input and output signals using a FIR lowpass ltersynthesized with the following parameters : fc 100 Hz, Fe 800 Hz,order 31, rectangular truncation window.22

2.3.4 Exercice 4 : second order lter cellPurely recursive cell1.TZy(n) x(n) a1 y(n 1) a2 y(n 2) Y (z) X(z) a1 z 1 Y (z) a2 z 2 Y (z)So :H(z) Y (z)1 1X(z)1 a1 z a2 z 22. This lter will be stable if the poles of H(z) have a modulus lowerthan 1 (enclosed in the circle with a radius of 1). Let's compute thepoles :H(z) z2z2 z 2 a1 z a2a(z z1 )(z z2 )We have to solve z 2 a1 z a2 0, leading to compute a21 4a2and consider the 3 possible cases :(a) 0 : two conjugated complex poles of modules a2 . In thatcase, it will be necessary to have a2 1. Note that we reallyhave here a second order cell.(b) 0 : a dual pole z0 a21 . In that case, we need a1 2.(c) 0 : two real poles (two rst order cells in fact) : z1 and z2 a1 .2 a1 2By pointing out that 1 z1 z2 1, weget the three following conditions : a2 a1 1, a2 a1 1 et 2 a1 2.By regrouping all the possible cases, we obtain, in the coe cientsplan (a1 , a2 ), the stability domain of the lter (a triangle) : see gure2.10.3. Filter frequency response as a function of a1 and a2 : 1 H (eω ) H(z)Hz z ej ωe2 leading to : H (eω ) 2 1 a21 a221, 2a1 (1 a2 ) cos(eω ) 2a2 cos(2eω)23

Figure 2.10 Stability triangle for second order IIR lters.where ωe is the normalized pulsation : ωe 2π Ffe .4. Search for the normalized resonance pulsation ωe0 : excepted for ωe 0and ωe 2π 0.5, the derivative of H (eω ) 2 is null when cos(eω0 ) 2) a1 (1 a. Thus, we obtain a resonance if4a22)5. cos(eω0 ) a1 (1 aand cos(2eω0 ) 24a2calculations : H(eω0 ) a1 (1 a2 )4a2a21 (1 a2 )216a2 1. 1 lead to, after some 2 a2p(1 a2 ) 4a2 a21For z1 rejθ and z2 re jθ (conjugated complex poles) : a2 r2and a1 2r cos(θ) and then : H(eω0 ) 1(1 r)(1 r) sin(θ)Indeed, H(ω0 ) is inversely proportional to the distance of the polesto the circle with radius 1. Closer is the pole to the circle (a2 close to 1)stronger will be the resonance. It is so possible to link the pole positionin the circle with radius 1 of the complex plane to the producedspectral e ect.24

6.H(z) 11z21z1 (1 z1 z 1 )(1 z2 z 1 )z1 z2 1 z1 z 1 z2 z1 1 z2 z 1T Z 1 h(n) z1n 1 z2n 1sin ((n 1)θ)u(n) rnu(n)z1 z2sin (θ)7. A way to implement this lter is proposed in gure 2.11.Figure 2.11 Purely recursive second order cell.General second order lter cell1.y(n) x(n) b1 x(n 1) b2 x(n 2) a1 y(n 1) a2 y(n 2)TZ Y (z) X(z) b1 z 1 X(z) b2 z 2 X(z) a1 z 1 Y (z) a2 z 2 Y (z)So :H(z) 1 b1 z 1 b2 z 21 a1 z 1 a2 z 22. H(z) 1 b1 z 1 b2 z 2 11 a1 z 1 a2 z 23. By implementing the FIR lter and the purely recursive second orderIIR cell in series ( gure 2.12), we can write :Y1 (z) X(z)11 a1z 1T Z 1 a2z 225 y1 (n) x(n) a1 y(n 1) a2 y(n 2)

T Z 1Y (z) Y1 (z) 1 b1 z 1 b2 z 2 y(n) y1 (n) b1 y1 (n 1) b2 y1 (n 2)This leads to a possible implementation : see gure 2.13. This implementation is called canonic structure : only one FIFO is needed byimplementing the IIR part of the lter followed by the FIR part.Figure 2.12 Cascade d'un ltre RIF et d'un ltre RII.Figureordre.2.13 Struture canonique de réalisation de la cellule du second4. For b2 1 we get : HF IR (z) 1 b1 z 1 z 2 . And thus :HF IR (eω ) e j ωe (b1 2 cos(eω )) Arg [HF IR (eω )] eωThe phase is linear and so the group propagation time (GPT) constant.It was also possible to note that hF IR [1 b1 1] is symmetrical, whichleads to a constant GPT for a FIR lter.2.3.5 Exercise 5 : guided IIR synthesis1. The IIR synthesis is a digital lter synthesis method based on ananalog synthesis. From the speci cations on H(f ), we rst have todeduct those on H(f ) to conduct the analog synthesis and get H(p).The speci cations in terms of attenuation are the same, but the fre1quency axis has to be pre-distorted : f πTe tan π fe (in order to26

anticipate the distortion introduced by the bilinear transformationgiving at the end H(z) from H(p)). This leads to the speci cationsgiven in gure 2.14.Figure 2.14 Analog speci cations2. Analog synthesis :(a) First step : we have to choose an analog model and x its parameters in order to obtain a transfer function H(f ) satisfying thespeci cations. We choose, in this exercise, to use a Butterworthlow pass model, whose transfer function is given by : H(ω) 2 1 1 2N ,ωωcwhere ωc represents the cuto pulsation ( 2πfc , fc denoting thecuto frequency). The attenuation for ω ωc 2πfc is xed to 3 dB. Parameter N has to be calculated in order to satisfy thewanted attenuation in the attenuated bandwidth :110 log101 ωaωc 2N 30dB 10 logwhere ωa 2π fea . This leads to choose N 3.(b) Second step : go from H(ω) 2 to H(p).hi H(ω) 2 H(p) 227p jω1103

so for N 3 : H(p) 2 1 H(p)H( p)1 p6(we will forget for the momentbypωc1ωc ,knowing that p will be replacedat the end).Among the 6 poles of H(p) 2 (which are the sixth roots of 1), 3belong to H(p), 3 belong to H( p). We will choose as the 3 polesfor H(p) : p0 1, p1 12 j 32 ,p2 12 j 32because theirreal parts are negative leading to a stable analog lter.It comes :H(p) so, by replacing p bypωcH(p) 1(p 1)(p2 p 1):ωc3(p ωc )(p2 pωc ωc2 )(c) Third step : After applying the bilinear transform on H(p), weobtain :H(z) H1 (z)H2 (z)withH1 (z) 0.43(1 z 1 ),1 0.29z 1H2 (z) 0.135 0.27z 1 0.135z 21 0.753z 1 0.4z 23. The obtained lter is composed of a stable rst order lter (its polehas a modulus 0, 29 1) and a second order lter, also stable becausethe point (a1 , a2 ) ( 0.753, 0.4) is included in the stability triangle(see previous exercise).4. Obtained lter is not a resonant lter becausea1 (1 a2 )4a2 6.588 1(see previous exercise).5. If we need to lter a signal x using the synthesize IIR lter, we canwrite the following code on Matlab :y1 f ilter([0.43 0.43], [1 0.29], x);y f ilter([0.135 0.27 0.135], [1 0.753 0.4], y1);28

Chapitre 3References Digital Signal Processing : Course, by Nathalie Thomas, ENSEEIHTdocumentation. Signal and Systems, by Simon Haykin and Barry Von Veen, Wiley Digital Signal Processing, by Alan V. Oppenheim, Ronald W. Schafer,Prentice-Hall. Documents on complex variable, Laplace transform, z tranform :http ml29

Digital Signal Processing Exercises with solutions Nathalie Thomas Master SATCOM 2018 2019. Chapitre 1 Exercises . signal sampling and limitation of the signal du- . z transform is the tool used to study the time invariant linear digital systems. 3. 1.2.1 Exercise 1 Let's consider a real a2]0;1[ and call u(n) the Heaviside function : .