Course On DSP Design Using MATLAB

Marcelo Basilio Joaquim, Jos Carlos Pereira and Vilma Alves de Oliveira Department of Electrical Engineering School of Engineering of So Carlos University of So Paulo CP: 359, 13560-970 So Carlos SP - Brazil
Abstract - In this paper is reported a graduate Digital Signal Processing laboratory course at the School of Engineering of So Carlos - University of So Paulo. The course is based on MATLAB for windows programming language to aid in the teaching of DSP. Through some examples, the students are introduced in a fast and natural way to the DSP concepts and practical applications. Discrete-time signals are represented by a sequence of numbers where each element represents the signal value at a determined point of time. The Matlab and its Signal Processing Toolbox have incorporated many commands and functions useful to generate basic signals used in DSP as well as complex signals composed by the basic ones. At the beginning of the course is introduced the generation of waveforms like sine, square wave, unit sample real exponential and random signals with uniform and gaussian distributions. By the use of these signals concepts of linear and circular convolution and sampling are introduced. An important point as far as discrete time signals is concerned is the clean understanding of the sampling theorem and its implications. As an example, a 8 kHz sinusoidal signal sampled at 32 kHz is used for this purpose. The sampling frequency is reduced to 8.6 kHz, resulting a 600 Hz sinusoidal signal, outlining the tendency of high frequencies to become low as a result of aliasing effect. In order to better illustrate the aliasing effect, a damped sinusoidal signal (second order system response) is used. This signal is sampled at several frequencies (16, 8 and 4 Hz) to demonstrate the spectrum superposition as shown in figure 1. This procedure in the course has helped the students to better understand the sampling theorem and meaning of the convolution.
S ig n a l S p e c trum 1 0 .8 0 .6 0 .4 0 .2 0 0 10 S a m p ling a t 8 H z 1 0 .8 0 .6 0 .4 0 .2 0 0 20 1 0 .8 0 .6 0 .4 0 .2 0 0 20 1 0 .8 0 .6 0 .4 0 .2 0 0 20 S a m p ling a t 4 H z S a m p ling a t 1 6 H z

Introduction
The main difficulty in DSP course is the large amount of mathematical models and equations. The repetitive algebraic and the complete understanding of the physical concepts embedded in the equations require the use of a suitable computational tool. The introduction of MATLAB in this course is due to its facility to build up mathematical functions and also due to its powerful graphical user interface in order to display the results. Important concepts on DSP are smoothly introduced by examples and graphical output. The following topics: Discrete-time signal generation; convolution; sampling theorem and aliasing effects on the signals; discrete linear systems and z- transform (impulse and frequency response, stability by pole analysis and series convergence, inversion of the z transform by partial fraction expansion and residue calculus; discrete Fourier transform and spectral analysis (windowing spectral resolution); FIR and IIR filter design; linear estimation and introduction to Wiener filters. In most examples the emphasis is given to practical application of the related theory. The exercises are built in such a way that the students can fix the important concepts on each experiment, and at the end of the course a final project is suggested evolving all concepts fixed before. The course has been organized in order to allow the students to design systems rather than just to solve exercises. Under this point of view the iterative procedure provided by MATLAB together with its high flexibility propitiates an ideal environment in which the students design DSP systems, verify results and change parameters, if necessary, to find out the best performance. This procedure inherently makes the students to incorporate the design quality in their tasks.

Discrete Signals, Convolution And Sampling

Figure 1. Aliasing Effect Of The Sampling.

Discrete Systems And Z-Transforms

Cumulative Sum of |h(n)| - (a) 4 1.5 1 3 2.5 2 1.5 1 -1 Imaginary part 0.5 Poles & Zeros Diagram - (b)

In this section, an introduction to linear discrete-time systems (LTI) is presented. Stability, causality, z-transform and its inverse using partial fraction expansion and residue theorem are studied. As we know, a LTI can be described in the time domain by a linear difference equation with constant coefficients like:

3.5

0 -0.5

a
k =0

y (n k ) = bk x (n k )
k =0

0.5

(1)

-1.5 0 20 40 60 80 100 -1.5 -1 -0.5 0 Real part 0.5 1

Or alternatively, in the frequency domain, the system can be represented by its system function H(z):

Figure 2. Cumulative Sum And Pole Locations. The inverse of the z-transform of a linear system can be obtained by expanding H(z) in partial fractions, i. e.,

H (z ) =

b0 + b1z + bM z a 0 + a1 z 1 + a N z N

(2)

H ( z) =

The impulse response, which is the inverse of the ztransform is another way to characterize the system. In the Matlab, this system is easily represented by two coefficients vectors, b=[ b0 b1 ... bM] and a=[a0 a1 ... aN]. The unit sample response can be observed using the unit sample function synthesized before and the filter function of the Matlab or directly from the impz function. Using sinusoidal signals is possible to verify frequency selective characteristics of a LTI. An important point is the verification of the system stability. A necessary condition for the system stability is that all its poles must be inside the unit circle of the z-plane. The functions root(a) and zplane(b,a) are used for this purpose. Another way to verify the system stability is to observe the convergence of the cumulative sum of h(n) absolute values. The frequency response of the system is determined by the freqz function. As an example considers the following system:

NM Rk + k jz j k k =1 1 p k z j =0 N

(4)

where Rk is the residue of pk and kj are the direct coefficients when N<M. Through this expansion the inverse of the z-transform can be obtained by inspection. The Matlab function [R P K] = residuez(b,a) is used to perform this calculus. A system where b = [1 2 1] and a = [ 1 -1.5 0.5] results in R = [8 -9]; P = [1 0.5] and k = [2]. Then, with these values, the unit sample response is:

h(n) = 8(1) n 9(0.5) n + 2 (n)

(5)

As can be observed, by its potentiality and friendly use, the software allows to characterize easily the LTI systems. Other examples are shown like: unstable systems, systems with poles with multiplicity higher than one, facilitating the theoretical DSP teaching.

y (n) + 0.2 y (n 1) + 0.8 y (n 2) = 0.3x (n) + 0.6 x (n 1) + 0.2 x (n 2)

where the poles are: -0.1 0.8888j In Matlab the following program part can be implemented to find the poles and their mapping in the zplane. The results are shown in figure 2. a=[1 0.2 0.8]; b=[0.3 0.6 0.2]; roots(a) zplane(b,a); cumsum(abs(h)); % Cumulative sum freqz(b,a,N,fs);

DFT And Spectral Analysis

The Discrete Fourier Transform (DFT) is one of the most important tools in DSP. It is used to represent a discretetime signal in the frequency domain. The DFT pair defined as:

X (k ) =

x ( n) e
n=0

N 1

jwkn

x ( n) = X ( k ) e
k =0

N 1

(6)
jwkn

For high N values, the equation (6) is computed fast and efficiently when FFT algorithms are used. Note that the

N length truncation of the sequence impose a periodic signal with period N, so special care must be taken when analyzing the results. In order to visualize the effects of the signal truncation an example is shown in figure 3. First, the DFT of a 3 cycles sinusoidal signal is computed (figure 3.a), then is calculated the DFT for the same signal with a non integer number of cycles (3.3 periods) as shown in figure 3.b. The signal truncation at a non integer number of periods produces spectral leakage.
3 periods of a Cossine wave (a) 0.5 0.45 0.4 0.35 0.3 0.25 0.2

parameter modification leads to a complete re-design. By other hand the computer aided design allows to overcome these difficulties as well as to visualize the results. Digital filters can be classified in two classes: IIR (Infinite impulse response) filter and FIR (Finite impulse response) filter. Normally, IIR filter are chosen due to its sharper frequency response, but, when linear phase characteristic is desired, only FIR filters can be used. The filter designs start from an analog counterpart followed by a transformation to the digital domain. IIR Filters In the section 3, three distinct way to represent a LTI has been shown: System function H(z), unit sample response and linear constant- coefficient difference equation. As far as the IIR filter design is concerned, these three representations lead to three different approaches of design, i e, impulse variance, bilinear transformation and derivative approximation. In the impulse invariance approach, the desired continuous impulse response is sampled as:

3.3 periods of a cossine wave (b) 0.45 0.4 0.5 0.35 0.3 0.25 0.3 0.2 0.15 0.2 0.4 0.6

TDF of High Resolution (c)

0.15 0.1 0.05 0 0 10 20 30 40

0.1 0.1 0.05 0 0 10 20 30 40 0 0 50 100 150 200 250 300

Figure 3. DFT Of A Sinusoidal Signal. Figure 3.c shows a higher spectral resolution for a better visualization of the signal spectrum, that can be performed by zero-padding in the time domain signal.
TDF Signal two tones (a) 0.6 0.18 0.16 0.5 0.14 0.4 0.12 0.1 0.3 0.08 0.2 0.06 0.04 0.1 0.02 0 0 10 20 30 40 0 0 10 20 30 40 Hanning Window (b)

h(n) = hd (nT )

(7)

High Resolution (c) 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0 50 100 150 200 250 300

where T is the sampling interval. The function system H(z) and, consequently, the filter coefficients are given by [1], [2], [5]:

H ( z) =

AK sk T 1 z k =1 1 e

(8)

Figure 4. DFT Using Windows The necessity to window the signal can be observed using a signal composed by two sinusoidal components. In the figure 4.a the signal is masked and it is difficult to identify clearly the spectral components of the signal. Figure 4.b shows the spectrum when the signal is windowed by a Hanning window, showing both spectral components. Increasing the spectral resolution (through a zero-padding in the time domain) the components are clearly identified as can see in the figure 4.c.

where sk are the poles in the s-plane and Ak are the coefficients of the partial fraction expansion. In the Matlab this design can be performed by the use of the impinvar function directly. The bilinear transformation is the most used method for IIR filter design. In this case the mapping of s into z is performed using the following transformation [3], [4]:

2 1 z 1 s= T 1 + z 1

and

H ( z ) = H ( s) s= 2 1 z 1
T 1+ z 1

(9)

The relationship between digital frequency and analog frequency for this transformation is:

Digital Filters
One of the most important class of LTI is the selective frequency filters. Their design by manual calculation evolve a large number of steps and calculus, and any filter

w = arctg T 2

(10)

It can be noted that this relationship applies a non linear compression on digital frequency. In order to

minimize this effect a pre-warping in s domain is required, i. e., a non linear expansion is applied for posterior compression.

Figure 5. Frequency Response Of The Filter FIR Filters The acronym FIR stands for finite impulse response filters, as a consequence, the system function (equation 2) has a0= 1 and ai = 0, i = 0, 1, ..., N. The most important property of these filters is that linear phase response can be obtained, which is impossible with IIR filters. By other hand, its great disadvantage is the large number of coefficients required for a given characteristic related to IIR filters. Some functions of Matlab used in the course for this topic are: fir1, fir2 and remez. The system filter function is:

' =

2 arctg T 2 T

(11)

As can be observed, the design of these filters evolve some phases which require a great deal of calculus like: specification of the filters characteristics, auxiliary analog filter design and H(s) to H(z) transformation. The functions buttord and butter of the Matlab can easily to design a filter with Butterworth characteristics. Furthermore, the Chebyshev and Elliptic characteristics are also incorporated in the software. As an example is shown a part of a program for a Butterworth digital filter design using bilinear transformation. disp(' --------------------------------------------'); disp('| Design of a Digital Low-pass Filter |'); disp(' --------------------------------------------'); % Data of the Analog Filter fap=input('Pass-band frequency: '); Rp=input('Maximun atenuation on pass-band: '); faa=input('Stop band frequency: '); Ra=input('Minimun atenuation on stop-band: '); fs=input('Sampling frequency: '); % Seletion of the filter order [N,wn] = buttord(2*fap/fs,2*faa/fs,Rp,Ra) % Filter Coefficients [bd,ad] = butter(N,wn) pause; % Frequency Response freqz(bd,ad,128,fs); title('Filtro Digital - Bilinear'); Example: Using fap=1.5 kHz, Rp=1 dB, faa=4kHz, Ra=40 dB and fs= 10kHz, results in the following filter coefficients: bd = [0.0834 0.2501 0.2501 0.0834] ad = [1.0000 -0.7355 0.4773 -0.0748]
Magnitude Response (dB Filtro Digital - Bilinear 0 -50 -100 -150

H ( z) =

h( n) z
n =0

N 1

(12)

It can be noticed that the filters coefficients is the unit sample response and the linear phase is implied when h(n) = h(M-n) [1]. [2]. The most used method to design FIR filters is windowing the data. The design begins by choosing an ideal filters characteristic hi(n). The response h(n) of the FIR filter is obtained by multiplying hi(n) by a finite length window w(n), i. e.,

h(n) = hi (n) w(n)

(13)

Figure 6. shows the frequency response for a low-pass filter by using rectangular and Hamming windows. It can be observed from this figure the ripple in the passband when a rectangular window is used. This effect comes from the high magnitude of the sides lobes in the spectrum of this window. If a Hamming is used the filter response is smoother but presents a larger transition band.
FIR Filter - Retangular and Hamming window 1.2 1 0.8

1000

2000 3000 Frequency (Hertz)

4000

5000

0.6 0.4

Phase (degrees)

0 -100 -200 -300

0.2 0

1000

2000 3000 Frequency (Hertz)

4000

5000

0.1

0.2

0.3

0.4

0.5

Figure 6. FIR Filters Also, other techniques are exemplified to the students like frequency sampling method and design using the Parks-McClellan algorithm.

comparison between the signal periodogram (DFT) and the results obtained.
P o w e r D e n s i ty S p e c trum a n d P e r i o d o g r a m 30

Linear Prediction And Wiener Filter

At this phase of the course some related topics of linear estimation and adaptive filters are presented. Particularly are studied the parametric estimation of the power spectral density for a stochastic signal (autoregressive model - AR) and an adaptive filtering of a embedded noise signal (Wiener filter). Linear Prediction Linear prediction is a parametric modeling technique used to describe a signal through time series. This technique is very efficient used in signal compression, voice analysis and synthesis and power spectrum estimation. Supposing x(n) to be a M-order autoregressive process, the linear prediction consists to estimate a non observed sample x(n) based on the past M observed samples {x(n-1), ..., x(n-M)}. The linear estimate of x(n), denoted by $ x(n) can be written as:

20

10

-10

-20

-30

10

20

30

40

Figure 7. DFT And Parametric Spectrum Wiener Filter The Wiener filter is an adaptive digital non recursive FIR filter designed to produce an optimal output in mean square sense. One of the filtering methods that are implemented in the course is shown bellow. Consider a stochastic signal x(n) composed by two components: a pure signal f(n) and an additive white noise v(n). Suppose also, that these signals are stationary in the wide sense.

$ x (n) = a1 x (n 1) a M x (n M )

(14)

where a(k) are constant coefficients and M is the order of the predictor. The coefficients a(k) are calculated based on the least square criterion, i. e., minimizing the mean square error $ e(n) = x(n) - x(n) . After some mathematical manipulation, the coefficients a(k) can be obtained from the following expression,

x ( n) = f (n) + v ( n)

(17)

The problem to be solved is to obtain the best estimate of f(n) when only x(n) and the joint statistics of x(n) and v(n) are known. A way to find out the filter coefficients is to use an adaptive algorithm [6] shown as follows: (1) Set an initial value for the filter coefficients

a ( j) R x ( k j) = Rx ( j), k = 1,2, M (15)

j =1

[h] = 0

(18)

where Rx(k) is the autocorrelation function of x(n) and the signal spectrum is given by:

(2) Compute the gradient vector [D(n)] which is the error derivative "minimum mean square (MMS)" related to the filter coefficients.

X ( z) =

1 E ( z) A( z )

(16)

[ D(n)] =

de(n) = 2 E {e(n) x (n k )] dh( k , n)

(19)

In order to solve the system equations (15) the Levinson-Durbin and Cholesky algorithms are presented and implemented. As an illustration of this model is determined, in the Matlab, the coefficients of the predictor filter and the power spectral density of an AR signal synthesized in the Matlab. In the figure 7. there is a

(3) Update the coefficients in a reversed order of the gradient vector.

[h(n + 1) = [h(n)] u[ D(n)]

(20)

where u is a constant which determine the correction step. (4) Compute the new error vector and repeat the process until to reach the optimum h(n). A good approximation of this algorithm is to replace the error in the step 2 by an instantaneous error called "least mean square (LMS)", [6] defined by:

solve related problems and beyond this, with the use of this program it was noted a higher progress and interest of the students in theoretical aspects of DSP and its practical applications.

References
1) Oppenheim, A. V. & Schafer, R. W., Discrete-time Signal Processing, Prentice Hall, 1989. 2) Proakis, J. G. & Manolakis, D. G., Digital Signal Processing: Principles, algorithms and applications, Macmillan Publishing Co., 1992. 3) The Mathworks Inc., Matlab Reference Guide, 1992. 4) The Mathworks Inc., Signal Processing Toolbox Users Guide, 1992. 5) Joaquim, M. B. & Pereira, J. C. Laboratrio de Processamento de Sinais, Classroom Notes, 1996. 6) Ifeachor, E. C. & Jervis, B. W., Digital Signal Processing: A Practical Approach, Addison-Wesley, 1993.

d (n) = 2e(n)[ x (n)]

(21)

In this approach, the coefficients are updated for any time as:

e(n) = f (n) [ x (n)]'[h(n)] [h(n + 1)] = [h(n)] + ue(n)[ x (n)]

(22)

The variable n is increased until a stable solution to be reached. Although this algorithm presents only an approximated solution of the optimum Wiener filter, it is simpler and can be realized in real time. This algorithm is implemented at the end of the course as an example of optimum filtering. A signal noised embedded with white noise is applied to the filter to observe the Wiener filter performance and the results are shown in figure 8.
Signal free of noise 1 0 -1

200

400

600

800

1000

Signal plus noise 1 0 -1

200

400

600

800

1000

Signal filtered 1 0 -1

200

400

600

800

1000

Figure 8. Wiener Filter Example

Conclusion
The introduction of DSP students to high level programming language is important in modern engineering curricula. As far as the students interested is concerned , the use of Matlab package in the experiments has been useful to assimilate the DSP concepts, to provide better capability to