Digital Signal Processing ESS040 Department of Electrical and Information Technique

Laboratory work Bengt Mandersson Lund 2012

2

Kapitel 1

Computer exercises Introduction to computer exercises The traditional way to motivate digital signal processing is the wish to process an analog signal digitally, i.e. using a computer or another form of digital hardware. Today digital signal processing is cheap and fast, can be accomplished on a small area, with great flexibility. The figure below shows the basic principle, where a real-world analog signal is sampled to a discrete-time signal consisting of a sequence of values or samples. The number of samples produced per second in this process is determined by the sampling frequency. The discrete-time signal can then be manipulated in different ways using digital signal processing, after which it is interpolated back to an analog signal, and transmitted to a proper analog system, e.g. a loudspeaker or an antenna. Sampling

Analogt Antivikningsfilter

Analog signal

Interpolation

S/H

A/D

Digital Signalbehandling

Tidsdiskreta signaler

D/A

Analogt Rekonstruktionsfilter

Analog signal

Already at an early stage of this course it is important to understand the principles and purposes of the different parts of this system. Therefore there follows a short summary: Discrete-time signals have a very fundamental property, namely they lack knowledge of at what rate they have been sampled or at what rate they will be reconstructed. Take the example of a CD disc containing a discrete-time signal. This signal does not know what time interval or sampling interval there should be between its samples, instead the CD player (the interpolator in the figure) provides the information that the signal is to be reconstructed with 44.1 kHz. This was also known in the studio where the disc was recorded (sampled). It is often said that a discrete-time signal has distance one between its samples, i.e. sampling period one without a unit. The sampling frequency is the inverse of the sampling period, and is hence also one without a unit. All frequencies of the discrete-time signal are given as a fraction of the sampling frequency, which is one. These frequencies are called normed frequencies and are between 0 and 1. It is only when the signal is reconstructed with a suitable sampling frequency the normed frequencies are transformed to “proper” frequencies, since it is then one decides what units are to be used for the time axis. Another very important property of discrete-time signals is the fact that the largest frequency 3

4

KAPITEL 1. COMPUTER EXERCISES

that can be represented is the frequency of a signal consisting of every second sample positive, and every second sample negative, with equal amplitude. This frequency has period 2 which is why the normed frequency is 0.5. One can try to record analog signals of faster variation, but these frequencies will not be represented correctly by the discrete signal. In fact, all frequencies of the analog signal which are above the normed frequency 0.5, corresponding to one half of the sampling frequency, will look like smaller frequencies in the discrete-time signal. This phenomenon is called aliasing, and causes the impossibility of remapping the signal back to the analog domain, since the aliased larger frequencies in the reconstruction process are considered to be in the interval 0 to 0.5. As an example we can consider two analog signals of 100 Hz and 350 Hz which are sampled at 500 Hz. The signal of 350 Hz is above half the sampling frequency and cannot be represented by the discrete-time signal, instead it is considered as a signal of frequency between 0 and 250 Hz (in this case 500350=150 Hz). If these two sampled signals are reconstructed we will get back the original 100 Hz tone, plus a tone of the aliased signal. The solution of this problem is naturally to only allow recording of signals which can be represented correctly of the discrete-time signal. Of this reason one often use an analog antialiasing filter, which lets through only frequencies smaller than half the sampling frequency, before a signal is transformed from the analog to the digital domain with aid of a sample-andhold circuit and an analog-to-digital converter. In this case the mapping between the analog and the digital domain is one-to-one, - the conditions of the so called sampling theorem are fulfilled. The fact that the conditions of the sampling theorem are fulfilled also means that we do not lose any information in the A/D converter, i.e. that the original analog signal can be reconstructed (interpolated) exactly from the discrete-time signal. This is accomplished by the digital-to-analog converter holding the value of each sample during a sampling interval at its output. The result is a piecewise constant signal, with stepwise changes at each new sampling interval. These stepwise changes introduce higher frequencies than half the sampling frequency. In fact all discrepancy from the original signal lies in frequencies above half the sampling frequency. Of this reason one again applies an analog filter which is called a reconstruction filter in order to eliminate frequencies above half the sampling frequency, which obviously were not sent into the system. Now the analog reconstructed signal looks exactly like the original analog signal, provided that we sampled twice as fast as the largest frequency of the input signal. If this had not been done, the anti-aliasing filter would have eliminated the large frequencies of the input signal, so the condition is fulfilled in any case then. Finally some words about the concept filter. A filter describes how a certain number of samples of a signal, up to the latest, shall be linearly combined in order to produce a certain output signal. The weights are put in a sequence forming a discrete-time sequence - this is the impulse response of the filter. The impulse response is a representation of a LTI system. Another description of the same system is as a difference equation. It is very important to realize that the concepts of filter and impulse response are not as strange as they first can seem. A filter is nothing but a tool in order to combine input data to obtain a desired output signal. Let us consider an example: Let us assume you each week save x[n] kronor on a bank account where n is the number of the week since you started. We call x[n] a signal - a discrete-time input signal to a system. Let us determine what output signal we desire; this will determine the system. As output signal from the system we choose for instance the sum of your savings the last month, i.e. the last four weeks. The difference equation for the system describes the output signal the current week as a sum of the four most recent savings where each of them is given weight 1, i.e: y[n] = 1x[n] + 1x[n − 1] + 1x[n − 2] + 1x[n − 3]. The filter now weights the last four savings and is therefore said to have length 4. The impulse response of the filter is the sequence of weights starting with the latest i.e. [ 1 1 1 1 ]. This type of filter has a bounded memory and weights only a limited amount of samples; the filter is therefore called

5

a FIR (Finite Impulse Response) filter. Let us instead assume you want your filter/system to compute your balance. A simple way to do that is to update the balance at each saving, i.e. the new balance is the old plus the saving amount. The difference equation for this system is y[n] = 1y[n − 1] + 1x[n]. It is important that not only the input signal is part of the computation, but also old output signals. This makes the computation of the impulse response more involved since it is, per definition, a rule for how to weight input signals in order to compute a certain output signal. If old output signals are part of the difference equation the system is called recursive. To solve the problem we see that the same equation can be used for old output signals, e.g. y[n − 1] = 1y[n − 2] + 1x[n − 1] and y[n − 2] = 1y[n − 3] + 1x[n − 2] etc. The important insight in this situation is that the output signal depends on all input signals (savings) up until now. The impulse response in this case is infinitely long [ 1 1 1 1 ... ... ] and is called an IIR (Infinite Impulse Response) filter. In this simple example we have only used weights equal to 1, but with different weights for input signals of different age the filter structure turns out to be very flexible. Each input signal can be divided into its frequency components, and in the rest of this course we will see how the filters can be used to enhance, attenuate, or delay different parts of these frequencies. In fact it is exactly the same thing to combine input signals from different time instances to an output signal as amplify, attenuate and delay the frequency components of a signal. Overview The detailed theoretical understanding of the different blocks in the system above will grow during the course. The computer exercises has three aims: 1) to illustrate applications of the theoretical concepts, 2) to describe how complicated systems can be built based on the basic components, and 3) to give examples of real-world signals from application areas where digital signal processing is used. Development in signal processing is often in the MATLAB environment. Therefore the computer exercises also aim at giving MATLAB routine and to teach its different signal processing tools and how to use them. The step from MATLAB development to a real-time algorithm in C, which can be run on a microprocessor or a digital signal processor, is treated in later courses at the department. 1 The heart and medicine - a system for recording/receiving/analysis 2 Digital audio tecnique (dedicated DSP hardware needed) 3 Design of a IIR-filter (dedicated DSP hardware needed) 4 Digital filtering of images.

6

1.1

KAPITEL 1. COMPUTER EXERCISES

Computer exercise 2: The heart and medicine

A system for recording/analysis/receiving: The example ECG-receiving The next system we will take a look at is shown in figure 1.1 (the system corresponds to the left half of figure ??) and describes how a measured analog signal is converted to digital data which is then processed for the purpose of e.g. information extraction.

Sampling Elektrod

Analogt Antivikningsfilter

Analog signal

Hjärtdata

S/H

A/D

Digital Signalbehandling

Tidsdiskret signal

Figur 1.1: System f¨or inspelning/analys/mottagning: Exempel EKG-m¨atning. This system is also very common and is used e.g. for recording of sounds for digital storage, to record other types of signals (for instance medical, electrical or mechanical) for analysis, and receiving transmitted data (if the microphone is replaced by an antenna). We call this system a recording/analysis/receiving unit. We repeat the functionalities of the different blocks: An analog signal is measured with aid of a microphone, electrode, antenna or transducer of any other kind. The signal has a frequency content up to a certain limit. The sampling theorem says that we must sample the signal (A/D convert) at a frequency which is twice as large as the largest signal frequency in order to get a correct discrete-time description of the analog signal. If we can not or do not want to sample at a such a high rate we can sample at a lower rate, provided the interesting information is in the lower frequency bands, but one must then eliminate frequencies above half the sampling frequency, in order to avoid aliasing and desctruction of the signal. This is accomplished with the anti-aliasing filter. Usually noise is received at measurements, and the noise may reside at high frequencies, which is why anti-aliasing filters are almost always used before A/D conversion. After this the signal is read during a sampling period in the S/H circuit so that the A/D converter gets time to read the value of each sample. In this exercise you will increase your understanding about a signal’s frequency representation and how the DFT tool can be used. The theory is illustrated with an ECG signal and your task is to analyze it in MATLAB and extract the meaning of the different frequencies. The analog (electrical) signal is sampled (recored) at the Department of Cardiology, Lund University Hospital, and stored digitally. Digital signal processing is performed off line by a computer (the analysis can also be made in real-time in order to use the system as an alarm).

An ECG signal The recorded ECG signal is sampled at 1 kHz i.e. all frequencies over 500 Hz have been eliminated. This is not a problem since the body is not expexted to produce faster variations than representable by approx. 250 Hz. Under Dator¨ ovningar at the course home page there is the following file to dowload to your MATLAB home directory: ekg1.mat.

1.1. COMPUTER EXERCISE 2: THE HEART AND MEDICINE

7

MATLAB exercises Exercise 1:

Download the ECG file to your MATLAB directory by writing

load ekg1 If you do not find your file in MATLAB you can change directory using dir and cd. Exercise 2:

Investigate the downloaded file by writing

size(ekg1) You find out that the file ekg1 consists of a discrete-time signal of 10000 samples that are stored as a row vector (1 row and 10000 columns). Using the sampling frequency 1000 Hz this means a length of 10 seconds. Investigate the file further by writing figure plot(ekg1) You see a ECG file of 10 seconds containing 10 heart beats, where the peaks correspond to the electrical impulse that contracts the heart muscle and squeezes the blood out to the body, and the following small peaks correspond return to the original configuration. You also see a slow movement in the base line, i.e. the amplitude level between the beats. This variation is caused by breathing, and also by ground level variation due to the measurement process at the body. The y axis unit is microvolts. Exercise 3: As mentioned the signal contains no information a bout how quickly it has been registered. MATLAB uses sampling enumeration as x axis for plots, as in the last exercise. It is always more intersting to have a time units on the x axis, i.e. to map the sample enumeration to time units. This can be accomplished using the sampling frequency (must be known) FT=1000; N=10000; n=0:N-1; t=n/FT; where F T is the sampling frequency, N is the signal length, n is the sample index counted from zero, and t is a new time indexation for the 10000 samples that is given by the total number of samples divided by the number of samples per second, i.e. a scale of seconds. Plot the signal again but now supply the the x axis unit plot(t,ekg1) xlabel(’tid (s)’) ylabel(’amplitud (mikroV)’) It can be seen that the pulse seem to be approx. 1/second i.e. 60 beats/minute. Exercise 4: Use the DFT tool to see what frequencies the signal contains; the heart beat frequency should be the strongest. In MATLAB the DFT is computed using the command FFT. This command needs to know the number of evenly spaced (normed) frequencies in the interval 0-1 to compute the spectrum. If you for instance choose to compute FFT for M 1 2 0 points, FFT will compute the spectrum at the normed frequencies [ M ... MM−1 ]. M M If you want to know the physical frequencies this correspond to you must multiply them with the sampling frequency. We choose to compute the FFT at M = 10000 points. M=10000; ekgspek=fft(ekg1,M);

8

KAPITEL 1. COMPUTER EXERCISES

The vector ekgspek contain both the amplitude and the phase function and is therefore complex valued. We are here interested only in which frequencies are in the signal, which is why we plot only the amplitude function (abs(ekgspek)). figure plot(abs(ekgspek)) Since we did not supply any x axis data the spectrum was plotted only against the enumeration of its samples. This information is difficult to interpret, but we do know the normed frequencies they correspond to, and since we know the sampling frequency we can transform the normed frequencies to the physical frequencies. f=(0:M-1)/M*FT; plot(f,abs(ekgspek)) xlabel(’frequency (Hz)’) An important observation: Frequencies above half the sampling frequency (500 Hz) are only a mirrored image of the frequencies below 500 Hz. It is hence below 500 Hz we are to study the figure. We zoom in at the left part (0-500 Hz and 0-800000 in amplitude). axis([0 500 0 800000]) It seems this signal is over-sampled, i.e. the prerequisites for the sampling theorem were filfilled with a great margin. Not much energy is above 50 Hz in the spectrum. We zoom, this time in both variables axis([0 20 0 150000]) We see the frequencies constituting the ECG signal. “Constituting” means that the signal mathematically can be broken down to these frequencies, not that the signal has been created by addition of a number of sinusoids with these frequencies. This is of course not the case; the signal has been created by many cells simultaneously changing their surface potential which give rise to an electrical field in the chest, measured with an electrode. We have now designed a computer interface to the physicians, and we would like to interpret the spectral plot. We see that the frequency zero has the largest energy. The reason is the varying base line which albeit not varying with a large amplitude compared to the beats, but yet is spread all over the signal. This spectral peak we are not interested in. The next top is almost exactly at 1 Hz. This frequency with harmonics at 2,3,4 Hz etc represents the peak waveform of the heart beats. The fundamental, i.e. the 1 Hz tone, corresponds to the repitition frequency of the signal. This frequency is thus a measure of the heart beat frequency. The computer program should therefore look for the first peak, overlooking the possible peak at 0 Hz. Exercise 5: Finally we shall investigate the way a simple filter may affect the ECG signal. We start by plotting the signal but this time at the upper half of a figure. figure subplot(211) plot(t,ekg1) We now try to decrease the noise of this signal by specifying that the output signal at time n is a weighted average of the fitfeen last samples. The filter is then fifteen equal coefficients, set e.g. to 0.2, for the input signals ekg1[n] to ekg1[n − 14]. The filtering is done using convolution (in MATLAB conv) of the input signal and the filter h=0.2*ones(1,15); y=conv(ekg1,h); y=y(15:length(y));

1.1. COMPUTER EXERCISE 2: THE HEART AND MEDICINE

9

The length of the convolution is the input signal length plus the filter length minus one. The last row takes away the transient at the beginning (the filter has not data at all its coefficients at the beginning), and makes the length of the output signal equal the input signal’s length. subplot(212) plot(t,y) We see a signal of considerably smaller noise. The signal is so clean that we can see a small peak before each of the great peaks. The small peaks correspond to atrial contraction which fill the ventricles with blood, so they are full when the great contraction takes place at the great peak.

10

1.2

KAPITEL 1. COMPUTER EXERCISES

Laboratury work 3: Digital audio

Introduction In this laboratory work we will use the DSK 6713 DSP kit from Texas Instruments. The programs in the DSK is started from Matlab.

Preparation Preparation exercise 1 Given the impulse response below. Determine the Fourier Transform.

h(n) =

1 [1 0 0 0 1] 3

(1.1) (1.2)

H(ˆ ω) =

X

h(k)e−j ωˆ k

(1.3)

k

Preparation exercise 2 Sketch the frequency response above. Use Matlab.

Matlab exercises FIR-filter i tidsplanet Solve the system below in Matlab (z −D delays the signal D sample). x(n)

y(t)

y(n)

A/D mikrofon

Simulate this in Matlab. Exemple of MATLAB code: t=-1:.0001;1; x=vco(t,[300,800],10000); soundsc(x,10000); noll=zeros(1,500); x1=[x,noll];x2=[noll,0.9*x]; y=x1+x2; soundsc(y,10000);

x(t)

D/A

z−D

högtalare 0.9

%generera 2 s l˚ ang tidsvektor, FS=1/10000 %generera en svept sinussignal mellan 300 och 800 Hz %lyssna %generera 500 nollor %generera ekot %lyssna, h¨ogtalare eller h¨orlurar m˚ aste vara anslutna

Chose now a real signal (speech) and listen to the echo effect.

1.2. LABORATURY WORK 3: DIGITAL AUDIO

1.2.1

11

Digital Audioteknik in real time

Task 1 Write first init digsig in Matlab. The echo processor delays the signal D samples and multiply the signal with the damping factor a, |a| < 1. Type echo menu.

Change the Delay and Gain. Talk in the microphone and listen. Try to estimate the value of the delay when it is difficult to talk. Turn gainM to zero Answer: ....... ms delay. Write the differense equation for y(n). Determine the poles and zeros and make a pole-zero plot. Determine H(z). Has the system linear phase? Explain. Stop the Matlab programm with exit. Task 2 Add some more delay lines shown in figure below.

Type echo menu and choose three delay lines. Adjust the parameters (D) and (a) using Delay and Gain, respectively. Talk in the microphone and listen. Determine the impuls response. Determine the impuls response if the number of delay lines go to infinity. Stop the program (the program in the DSP runs until a new program is downloaded).

12

KAPITEL 1. COMPUTER EXERCISES

Optional task Make an realization of the system with infinite numbers of delay lines.

Run echo menu and choose IIR-filter. Task 3 Now, try a little more advanced system of delay lines with time varying parameters

Call echo chorus menu to start the program. Try to find suitable values of the parameters Offset, Swing and Freq. Discuss your results with the supervisor

1.2. LABORATURY WORK 3: DIGITAL AUDIO

13

Program in the DSP Description of echo menu

Command: global cc, echo menu, cc is the Matlab adress to the DSP (global). Funktioner Download Download the parameters to the DSP. Open DSP Starts the CCS, open a Matlab link to the DSP, download and starts the program in the DSP. Close DSP Close the DSP (Close the Matlab link to the DSP. The program in the continues until a new prgram is downloaded. FIR Choose FIR-filter (1 to 3 delay lines). IIR Choose IIR-filter (1 delay line). Exit Exit the Matlab program. (The Matlab link and DSP are not affected).

14

KAPITEL 1. COMPUTER EXERCISES

Description of echo chorus menu

Command: global cc, echo chorus menu, cc ¨ar Matlab-address till DSPn. Funktioner Download L¨aser ner parametrar till DSPn. Open DSP Startar CCS, ¨oppnar Matlab-l¨ank till DSPn, laddar ner och startar program i DSPn. Close DSP St¨anger DSP (St¨anger Matlab-l¨anken till DSPn. Programmet i DSP forts¨atter tills nytt program laddas ner). Exit Avslutar Matlab-programmet (Matlab-l¨anken och programmet i DSPn p˚ avekas ej).

1.3. LABORATORY WORK 4: FILTER DESIGN OF IIR FILTERS

1.3

15

Laboratory work 4: Filter design of IIR filters

1.3.1

Introduction

In this laboratory work we will use the DSK 6713 DSP kit from Texas Instruments. The programs in the DSK is started from Matlab.

1.3.2

Preparation

Preparation exercise 1 Combine the pole-zero plot and the magnitude spectra |H(f )| below.

Imaginary Part

Imaginary Part

Imaginary Part

Imaginary Part

Amplitudfunktioner |H(f)| 1

Z−plan 1 0 −1 −2

0.5 0 Real Part

2

0

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

1 1 0 −1 −2

0.5 0 Real Part

2

0 1

1 0 −1 −2

0.5 0 Real Part

2

0 1

1 0 −1 −2

0.5 0 Real Part

2

0

Preparation exercise 2 Sketch the magnitude spectra |H(f )| for the given the pole-zero plots.

16

KAPITEL 1. COMPUTER EXERCISES

1.3. LABORATORY WORK 4: FILTER DESIGN OF IIR FILTERS

1.3.3

17

Laboratory tasks

IIR filter The program MKIIR is used for interactive IIR filter design.

Task 1 Place a pole or a pair of complex poles in the z plane and check the result. Type MKIIR to start the program. Press first Replace last. Move the pole inside the unit circle and check the result. Investigate the relation of the angle of the pole and the magnitude spectra and the impulse response. Release Replace last.

18

KAPITEL 1. COMPUTER EXERCISES

Task 2 Choose poles and zeros for fulfilling the demands below. Press Import Workspace and open the file iirspec1.mat to load the specifications into MKIIR.

|H(f )|

6¡¡¡¡¡¡¡¡¡¡ ¡¡¡¡¡¡¡¡¡¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ −3 ¡ ¡¡ ¡¡ ¡¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡¡ ¡¡ ¡¡ ¡¡ ¡¡ ¡¡ ¡ ¡ ¡ −40 ¡ ¡ ¡ ¡ 0

0.1

0.15

f

Choose poles and zeros to fulfill the requirements. Save your filter with Export Workspace and choose a file name (i.e. filter1.mat. Then press Import filter. Optional task Compare with standard filter solutions. Stop the program MKIIR, press Exit. Task 3 Now, you will exame notch filters. The system function for both the FIR filter and the IIR filter are given below HF IR (z) = (1 − ejω0 z −1 )(1 − e−jω0 z −1 ) = 1 − 2 cos(ω0 )z −1 + z −2 HIIR (z) =

(1 − ejω0 z −1 )(1 − e−jω0 z −1 ) 1 − 2 cos(ω0 )z −1 + z −2 = (1 − rejω0 z −1 )(1 − re−jω0 z −1 ) 1 − 2r cos(ω0 )z −1 + r2 z −2

with ω0 the frequency which shall be cancelled. Plot the magnitude spectra for the filters with varying pole positions using the program MKIIR . Alternatively,you can use Matlab (freqz, poly, help). Task 4 In this task, you shall cancell sinousoids added to music files. Sampling rate 8000 Hz, stereo. Call global S, jukebox ) and choose a song from the list. The variabel S will contents 100 000 samples from the song. Use Windows Soundrecorderto listen to the song. Avoid Matlab due to difficulties to stop the outputs in Matlab.

1.3. LABORATORY WORK 4: FILTER DESIGN OF IIR FILTERS

19

Now, use MKSPA to determine the frequencies of the sinousoids. Press PEAK. The frequencies are multiples of 10 Hz. Press Exit to stop MKSPA. Now, you shall determine a notch filter which cancell the sinousoids. Try both FIR and IIR filters. Type your choice of poles and zeroes in Matlab and put them into vectors. Then, determine the A and B polynomes. Example: z1=exp(-j*2*pi*F0/Fs); z2=conj(z1); Z=[z1,z2]; Z=[ z0 z1 . . . z5 ]; P=[ p0 p1 . . . p5 ]; B=poly(Z); A=poly(P ); You must use j for the imagin part Test the filter using MKIIR. Then, download the filter to the DSP. Type (global cc B A, filter menu). Press Open DSP and wait untill the DSP has started. Play your song and test your filter using Download. Test both the FIR and the IIR filter.

20

1.3.4

KAPITEL 1. COMPUTER EXERCISES

Programbeskrivningar

IIR-filterdesign, MKIIR Beskrivning MKIIR m¨ojligg¨or interaktiv design av IIR-filter, genom utplacering av poler och nollst¨allen i z-planet med hj¨alp av muspekaren. Amplitudfunktion, fasfunktion och trunkerat impulssvar kan unders¨okas. Om kravspecifikation ¨ar angiven kan standardfilter fr˚ an Matlab Signal Processing Toolbox importeras. Programmets funktioner beskrivs i korthet nedan.

Funktioner Import Workspace Inl¨asning fr˚ an fil av kravspecifikation, och i f¨orekommande fall, filter. Export Workspace Lagring till fil av kravspecifikation, och i f¨orekommande fall, filter. Set filter specs Inmatning av kravspecifikation; Passbands- och sp¨arrbandsfrekvens (normerad frekvens), passbandsrippel och sp¨arrbandsd¨ampning (dB). F¨or att helt ta bort kravspecifikation, ange samma passbands- och sp¨arrbandsfrekvens. Filter coeffs Inmatning eller utskrift av filterkoefficienter. Help Ej implementerad ¨annu. About Information om programversion. Reset Rensar alla poler och nollst¨allen, men bevarar eventuell kravspecifikation. Exit Avslutar programmet. Sparar arbetsytan till filen current.mat .

1.3. LABORATORY WORK 4: FILTER DESIGN OF IIR FILTERS

21

Add pole Utplacering av poler i z-planet. Peka p˚ a omr˚ adet innanf¨or enhetscirkeln, tryck v¨anster musknapp f¨or att placera pol. F¨or att avsluta funktionen, tryck h¨oger musknapp. Delete pole Borttagning av poler i z-planet. Peka p˚ a den pol som ska tas bort, tryck v¨anster musknapp. F¨or att avsluta funktionen, tryck h¨oger musknapp. Add zero Utplacering av nollst¨allen i z-planet. Peka p˚ a omr˚ adet p˚ a eller innanf¨or enhetscirkeln, tryck v¨anster musknapp. F¨or att avsluta funktionen, tryck h¨oger musknapp. Delete zero Borttagning av nollst¨ allen i z-planet. Peka p˚ a det nollst¨alle som ska tas bort, tryck v¨anster musknapp. F¨or att avsluta funktionen, tryck h¨oger musknapp. Move Flytta poler eller nollst¨allen i z-planet. Peka p˚ a pol eller nollst¨alle som ska flyttas, tryck v¨anster musknapp. Det valda objektet markeras med svart. Peka p˚ a ny position, tryck v¨anster musknapp. F¨or att avsluta funktionen, eller avbryta omplacering, tryck h¨oger musknapp. Import filter Import av standardfilter fr˚ an Matlab Signal Processing Toolbox. F¨oruts¨atter att kravspecifikation finns. Observera att Butterworthfiltret i Matlab Signal Processing Toolbox anpassas efter sp¨arrbandsd¨ampningen, och ej n¨odv¨andigtvis sk¨ar genom −3 dB-punkten. Magnitude Visar filtrets amplitudfunktion. Phase Visar filtrets fasfunktion. Impulse Visar filtrets impulssvar. F¨or IIR-filter visas det trunkerade impulssvaret. Replace last Det senast utlagda nollst¨allet tas bort d˚ a ett nytt l¨aggs ut. Fungerar p˚ a motsvarande s¨att f¨or poler. Reciproc Zeros Anger att reciproka nollst¨allen ska l¨aggas ut. Fungerar endast f¨or nollst¨allen. HighPass Anger att filtret ska vara av h¨ogpasstyp. Styr den konstanta f¨orst¨arkningen. F¨or h¨ogpassfilter g¨aller att |H(f )||f =0.5 = 1, och f¨or l˚ agpassfilter att |H(f )||f =0 = 1. Zoom Anger att Matlabs zoom ska vara aktiverad. Unwrap Phase Anger att fasen ska ritas upp utan att begr¨ansas till intervallet [−180◦ , 180◦ ]. Focus Zoomar in ¨overg˚ angszonen i amplitudfunktionen. F¨oruts¨atter att kravspecifikation finns. Radius Zoomar in z-planet inom den radie som anges till h¨oger. Adjust Gain Justerar den konstanta f¨orst¨arkningen. Anges i dB. ¨ Ovrigt Muspekarens aktuella position i z-planet visas som radie och normerad frekvens (vinkel) vid Radius och Frequency. Antal utlagda nollst¨allen och poler visas vid Zeros respektive Poles. Den konstanta f¨orst¨arkningen som kr¨avs f¨or att |H(f )||f =0 = 1 f¨or l˚ agpassfilter och |H(f )||f =0.5 = 1 f¨or h¨ogpassfilter, visas vid Gain. Tips! Poler kan endast placeras innanf¨or enhetscirkeln. F¨or att l¨agga en pol s˚ a n¨ara enhetscirkeln som programmet till˚ ater, sikta strax utanf¨or cirkeln.

22

KAPITEL 1. COMPUTER EXERCISES

Nollst¨allen kan placeras p˚ a och innanf¨or enhetscirkeln, om man bortser fr˚ an reciproka nollst¨allen. F¨or att placera ett nollst¨alle p˚ a enhetscirkeln, sikta strax utanf¨or, s˚ a l¨aggs det p˚ a enhetscirkeln. Om olika poll¨agen/nollst¨allesl¨agen ska unders¨okas, anv¨and funktionen Replace last, s˚ a ers¨atts senast angivna pol/nollst¨alle d˚ a ett nytt l¨aggs ut. Om det skulle intr¨affa att ett f¨onster inte st¨angs, g¨or d˚ a s˚ a h¨ar: Peka p˚ a det f¨onster som ska st¨angas, och tryck v¨anster musknapp. Skriv d¨arefter delete(gcf) vid promten.

1.3. LABORATORY WORK 4: FILTER DESIGN OF IIR FILTERS

23

Spektralanalys, MKSPA Beskrivning MKSPA anv¨ands f¨or att analysera frekvensinneh˚ allet i ljudfilerna. Programmet ¨ar uppbyggt kring MATLAB-funktionen psd, med till¨agg f¨or markering av toppar i spektrum.

Funktioner Window Lista f¨or val av f¨onsterfunktion. FFT length Lista f¨or val av FFT-l¨angd. Overlap Inmatning av hur m˚ anga sampel segmenten ska ¨overlappa varandra med. Fs(Hz) Inmatning av samplingsfrekvens, som ligger till grund f¨or graderingen av frekvensaxeln. Zoom Anger om Matlabs zoom ska aktiveras. Analyze Ber¨akning av spektrum, baserat p˚ a vald f¨onsterfunktion, FFT-l¨angd, och antal overlappande sampel. ¨ Peak Markerar toppar i spektrum. Genom att flytta muspekaren ¨over grafen, markeras den st¨orsta toppen i n¨arheten av muspekaren. Amplitud (dB), samt frekvens (Hz) f¨or aktuell topp visas vid textf¨alten Peak:. Den aktuella toppen markeras med en triangel. Exit Avslutar programmet.

24

KAPITEL 1. COMPUTER EXERCISES

Filtrering i DSPn, filter menu Beskrivning Programmet filter menu:s funktioner beskrivs i korthet nedan.

Kommando: global cc B A, filter menu, cc ¨ ar Matlab-address till DSPn. Funktioner Download L¨ aser ner parametrar till DSPn. Filtret best¨ ams av B(z) och A(z). Open DSP Startar CCS, ¨ oppnar Matlab-l¨ ank till DSPn, laddar ner och startar program i DSPn. Close DSP St¨ anger DSP (St¨ anger Matlab-l¨ anken till DSPn. Programmet i DSP forts¨ atter tills nytt program laddas ner). FIR V¨aljer FIR-filter. IIR V¨aljer IIR-filter. Exit Avslutar Matlab-programmet (Matlab-l¨ anken och programmet i DSPn p˚ avekas ej).

1.4. COMPUTER EXERCISE 5: FILTERING OF IMAGES

1.4

25

Computer exercise 5: Filtering of images

In this computer exercise wi will apply digital filtering into image processing.

1.4.1

Tasks

Task 1 Apply a sliding filter to some various waveforms. Test with a square wave and some various sinusoids. Plot the results. Start with the example below.

>> >> >> >> >> >> >> >>

N=256; fhat=12/N; %Normerad frekvens f¨ or fyrkantssignalen n=(1:N);xx=square(2*pi*fhat*n); %skapar en fyrkantssignal som insignal L=5; %filterl¨ angd h=1/L*ones(L,1); yy=conv(h,xx); % subplot(211), stem(yy); %rita upp vektorn yy subplot(212), stem(xx); %rita upp vektorn xx

1.4.2

Filtering of images

Task 2 The Matlab program below load an image and determine a sliding average for the pixels.

>> >> >> >>

xx=double(imread(’lena.tif’))./255; %L¨ aser in bilden, skalar pixelv¨ ardena N=9; %S¨ atter maskstorleken hh=ones(N)./(N*N); % Skapar en kvadratisk mask med sidan N yy=conv2(xx,hh); %Ber¨ aknar 2-D faltningen mellan bilden och filtret

Type colormap gray. Then, plot the image using imagesc. Add some noise and redo the filtering.

>> xx_brus=xx+rand(size(xx))*0.2; % skapa en brusig version av xx Change the filter order and try again.

1.4.3

High pass FIR filter

Task 3 Redo task 1 with the high pass filter below. 1 1 1 4 1 1 h(n) = [00100] − [11111] = [− − + − − ] 5 5 5 5 5 5

(1.4)

26

KAPITEL 1. COMPUTER EXERCISES

Task 4 High pass filtering of images.  0 0 h(k, l) =  0 1 0 0

  0 1 1 1 0 −  1 1 9 0 1 1

   1 −1 −1 −1 1 1  =  −1 8 −1  9 1 −1 −1 −1

Apply this mask on the image used above. Used either imagesc or image.

(1.5)