LUCRAREA 2 SEMNALE ŞI SISTEME ÎN TIMP DISCRET · fft. Denumirea sa reprezintă prescurtarea de la...

2. SEMNALE ŞI SISTEME ÎN TIMP DISCRET

27

LUCRAREA 2 SEMNALE ŞI SISTEME ÎN TIMP DISCRET 2.1. Semnale în timp discret Un semnal în timp discret este definit prin valorile acestuia măsurate la momente discrete de timp. Semnalele în timp discret sunt reprezentate matematic prin secvenţe de numere notate:

21,][ NnNnx ≤≤ (2.1)

În MATLAB aceste secvenţe se pot defini ca vectori linie sau coloană, având elemente reale sau complexe. O primă limitare apare din faptul că aceşti vectori sunt de lungime finită în timp ce în problemele de prelucrarea numerică a semnalelor se poate lucra cu secvenţe de lungime infinită.

2.1.1. Definirea semnalelor în timp discret În studiul semnalelor şi sistemelor în timp discret se utilizează câteva secvenţe de bază ce vor fi prezentate în continuare, împreună cu modul lor de definire în MATLAB. • Impulsul unitate

Din punct de vedere matematic este definit astfel:

⎩⎨⎧

≠=

=0,00,1

][nn

nδ (2.2)

Utilizând proprietatea de deplasare în timp se poate scrie că

⎩⎨⎧

≠=

=−0

00 ,0

,1][

nnnn

nnδ (2.3)


28

Deoarece în MATLAB nu putem defini secvenţe de lungime infinită trebuie precizat domeniul de valori pentru n. Pentru a facilita definirea unor secvenţe de acest tip vom crea o funcţie MATLAB (vezi lucrarea 1, secţiunea 1.2.8.): function [y,n]=impuls(li,ls,k); %IMPULS - impulsul unitate (Dirac) in timp discret, %definit pe un suport temporal finit %-sintaxe: %y=impuls(li,ls,k) %[y,n]=impuls(li,ls,k) %-parametrii de iesire: %y=vector linie ce reprezinta delta(n-k) pe suportul %[li;ls] %n=vector linie ce reprezinta suportul [li;ls] %-parametrii de intrare: %li=limita inferioara a suportului temporal; %ls=limita superioara a suportului temporal; %k=indicele din delta(n-k) %-pentru afisare: stem(n,y) if nargin<3 error('Prea putine argumente de intrare') elseif nargin>3 error('Prea multe argumente de intrare') end if nargout>2 error('Prea multe argumente de iesire') end if li>=ls error('Suportul temporal este invalid') end if (k<li)|(k>ls) error('Indicele nu apartine suportului temporal') end L=ls-li+1; y=zeros(1,L); y(k-li+1)=1; n=li:ls;


29

Exemple: Utilizând funcţia MATLAB impuls creată anterior, să se definească şi să se reprezinte grafic secvenţele: 1. ][][1 nnx δ= 2. 2[ ] 0,5 [ 3]x n nδ= − , pentru 1010 ≤≤− n .

[y1,n]=impuls(-10,10,0) stem(n,y1),grid [y2,n]=impuls(-10,10,3) stem(n,0.5*y2),grid

• Treapta unitate

Din punct de vedere matematic este definit astfel:

⎩⎨⎧

<≥

=0,00,1

][nn

nu (2.4)

Utilizând proprietatea de deplasare în timp se poate scrie că

⎩⎨⎧

<≥

=−0

00 ,0

,1][

nnnn

nnu (2.5)

Se poate crea de asemenea o funcţie MATLAB pentru definirea secvenţelor de tip treaptă unitate:

function [y,n]=treapta(li,ls,k);

%TREAPTA - treapta unitate in timp discret, definita %pe un suport temporal finit %-sintaxe: %y=treapta(li,ls,k) %[y,n]=treapta(li,ls,k) %-parametrii de iesire: %y=vector linie ce reprezinta u(n-k) pe suportul %[li;ls] %n=vector linie ce reprezinta suportul [li;ls] %-parametrii de intrare: %li=limita inferioara a suportului temporal; %ls=limita superioara a suportului temporal; %k=indicele din u(n-k) %-pentru afisare: stem(n,y)

if nargin<3 error('Prea putine argumente de intrare') elseif nargin>3 error('Prea multe argumente de intrare') end


30

if nargout>2 error('Prea multe argumente de iesire') end

if li>=ls error('Suportul temporal este invalid') end

if (k<li)|(k>ls) error('Indicele nu apartine suportului temporal') end

L=ls-li+1; y=zeros(1,L); y(k-li+1:L)=1; n=li:ls;

Exemple: Să se definească şi să se reprezinte grafic secvenţele: 1. ][][1 nunx = 2. 2[ ] 0,7( [ 3] [ 3])x n u n u n= + − − , pentru 105 ≤≤− n .

[y1,n]=treapta(-5,10,0) stem(n,y1),grid [y21,n]=treapta(-5,10,-3) [y22,n]=treapta(-5,10,3) y2=0.7*(y21-y22) stem(n,y2),grid E1. Exerciţii: Să se definească şi să se reprezinte grafic următoarele secvenţe: 1. 1[ ] [ ] [ 1]x n n nδ δ= − − pentru 1010 ≤≤− n 2. 2[ ] [ ] 0,5 [ 1] 0,3 [ 2] 2 [ 1]x n n n n nδ δ δ δ= − − + − − + pentru 1010 ≤≤− n 3. 3[ ] [ ] 0,5 [ 4] 0,5 [ 4]x n u n u n u n= + − − + pentru 1010 ≤≤− n 4. 4[ ] [ 1] [ 5] [ 2] 2 [ 9]x n n u n n nδ δ δ= − + − + + − − pentru 2010 ≤≤− n

5. 11 1[ ]2 2

n n

x n ⎛ ⎞ ⎛ ⎞= − −⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠

pentru 100 ≤≤ n

6. 2[ ] ln cos sin15 15n nx n π π⎛ ⎞ ⎛ ⎞= −⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠ pentru 2020 ≤≤− n

7. 1[ ] ( 1) cos15

n nx n π⎛ ⎞= − ⎜ ⎟⎝ ⎠

pentru 100 ≤≤ n


31

2.1.2. Convoluţia liniară a semnalelor în timp discret Produsul de convoluţie liniară a două secvenţe numerice x1[n] şi x2[n] este definit astfel:

∑∞

−∞=

−⋅=∗=k

knxkxnxnxnx ][][][][][ 2121 (2.6)

unde secvenţele x1[n] şi x2[n] s-au presupus a avea lungime infinită. În MATLAB secvenţele numerice sunt privite ca nişte vectori şi ţinând cont de faptul că ele sunt finite şi că indicele primului element dintr-un vector nu poate fi zero (se efectuează o indexare cu 1) atunci definiţia de mai sus devine:

∑−

=

−⋅+=+1

021 ][]1[]1[

N

kknxkxnx (2.7)

unde N este maximul dintre lungimile celor două secvenţe. Sintaxa: conv(x1,x2) • returnează ca rezultat un vector de lungime egală cu lungimea vectorului x1

plus lungimea vectorului x2 minus 1, ce reprezintă produsul de convoluţie liniară al celor două secvenţe definite prin vectorii x1 şi x2.

Exemplu: Să se calculeze şi să se reprezinte grafic produsul de convoluţie liniară a secvenţelor ]5[][][1 −−= nununx ( 100 ≤≤ n ) şi nnx )9,0(][2 = ( 200 ≤≤ n ). x1=treapta(0,10,0)-treapta(0,10,5); n=0:20; x2=0.9.^n; x=conv(x1,x2); subplot(2,2,1),stem(0:10,x1),title(’x1’),grid subplot(2,2,2),stem(n,x2),title(’x2’),grid subplot(2,1,2),stem(0:length(x)-1,x),title(’x’),grid

0 5 100

0.2

0.4

0.6

0.8

1x1

0 5 10 15 200

0.2

0.4

0.6

0.8

1x2

0 5 10 15 20 25 300

1

2

3

4

5x


32

E2. Exerciţii: Să se calculeze şi să se reprezinte grafic (ca în exemplul precedent) produsul de convoluţie liniară a următoarelor perechi de secvenţe: 1. ]15[]5[][1 −−−= nununx nnx −= 10][2 pentru 200 ≤≤ n 2. ( )1[ ] sin 0.1x n nπ= ⋅ ⋅ 2[ ] [ 5]x n nδ= − pentru 200 ≤≤ n

3. 1[ ] cos(0,25 )x n nπ= ⋅ nnx )1(][2 −= pentru 200 ≤≤ n 2.1.3. Transformata Fourier discretă Transformata Fourier în timp discret (DTFT – Discrete Time Fourier Transform) a unei secvenţe ][nx este dată de relaţia:

∑ −=n

njj enxeX ωω ][)( (2.8)

unde ω este pulsaţia normată: 2S

FF

ω π= .

F este frecvenţa nenormată (exprimată în Hz) iar FS este frecvenţa de eşantionare.

Similar, frecvenţa normată este: S

FfF

= .

Funcţia )( ωjeX este periodică de perioadă π2 , deci este suficient să cunoaştem comportarea sa în intervalul ),[ ππ− (interval de bază). Datorită faptului că această funcţie este continuă, variabila ω putând lua o infinitate de valori, nu este posibilă o implementare pe o maşină de calcul.

Pentru a realiza totuşi o analiză în frecvenţă se utilizează transformata

Fourier discretă TFD (DFT – Discret Fourier Transform), obţinută prin discretizarea variabilei ω pe intervalul )2,0[ π în N puncte:

2kkN

ω π= , cu 1,,1,0 −= Nk … .

Astfel, transformata Fourier discretă a unei secvenţe ][nx este dată de

relaţia:

2

[ ] [ ]j kn

N

n

X k x n eπ

−=∑ cu 1,,1,0 −= Nk … (2.9)


33

În figură sunt prezentate reprezentările spectrului unui semnal discret în funcţie de pulsaţie sau frecvenţă normată şi corespondenţa cu frecvenţa analogică. De asemenea se observă corespondenţa între componentele spectrale de indice k calculate cu TFD şi spectrul reprezentat în pulsaţii normate.

În MATLAB, pentru calculul transformatei Fourier discrete se foloseşte funcţia fft. Denumirea sa reprezintă prescurtarea de la Fast Fourier Transform (transformata Fourier rapidă) şi indică faptul că este folosit pentru calcul un algoritm rapid. Sintaxe: y = fft(x) • dacă x este un vector se returnează un vector y de aceeaşi dimensiune cu

vectorul x ce conţine valorile transformatei Fourier discrete aplicată elementelor vectorului x; lungimea transformatei Fourier (numărul de puncte

-Ωmax 0 Ωmax

-Fmax 0 Fmax FS/2 FS

( )jaX Ω

2SΩ

max2SΩ > Ω Ω[rad/s]

F[Hz]

-π -ωmax 0 ωmax π 2π

-0.5 -fmax 0 fmax 0.5 1

( )jX e ω

ω

f

SFω Ω=

S

FfF

=

0 123 N/2 N-1

( )X k

k

2kkN

ω π=0,..., 1k N= −

2SΩ

−


34

N în care se calculează transformata) este egală în acest caz cu lungimea vectorului x.

• dacă x este o matrice se va returna matricea y de aceeaşi dimensiune cu matricea x; coloana i din matricea y va conţine valorile transformatei Fourier discrete aplicată elementelor coloanei i din matricea x.

y = fft(x,N)

• aceleaşi considerente ca în sintaxa precedentă cu deosebirea că în acest caz se specifică şi numărul de puncte N în care se calculează transformata.

Exemple: Să se calculeze transformata Fourier discretă a secvenţei:

]10[][][1 −−= nununx , pentru 0 20n≤ ≤ Să se reprezinte grafic partea reală, partea imaginară, modulul şi faza transformatei Fourier discrete calculate. x1=treapta(0,20,0)-treapta(0,20,10); X=fft(x1); figure(1),plot(X) // se obţine o reprezentare grafică eronată deoarece X e un vector cu valori complexe. figure(2),plot(abs(X)) // calculul transformatei Fourier discrete s-a efectuat într-un număr de puncte egal cu lungimea vectorului x1 (vezi sintaxa de la fft); pentru o mai bună reprezentare vom efectua calculul într-un număr mai mare de puncte (N = 512): X1=fft(x1,512); figure(3),plot(abs(X1)),grid // Această reprezentare corespunde însă intervalului de frecvenţă )2,0[ π (pe abscisă avem numărul de puncte al vectorului X1 pentru că în sintaxa funcţiei plot nu s-a specificat nimic altceva). De obicei se doreşte însă reprezentarea în intervalul de bază ),[ ππ− . Având în vedere faptul că funcţia este periodică de perioadă π2 atunci reprezentarea din intervalul )0,[ π− corespunde cu reprezentarea din intervalul )2,[ ππ . Prin urmare trebuie realizată o inversare a

0 100 200 300 400 500 6000

1

2

3

4

5

6

7

8

9

10


35

celor două jumătăţi ale vectorului X1. Acest lucru se realizează în MATLAB prin utilizarea comenzii fftshift (schimbă între ele cele două jumătăţi ale unui vector). În plus, pentru a avea pe abscisă reprezentarea în intervalul ),[ ππ− trebuie generat un vector cu pas liniar care să conţină în intervalul respectiv un număr de elemente egal cu lungimea transformatei Fourier discrete calculate:

w=-pi:2*pi/512:pi-2*pi/512; figure(4),plot(w,fftshift(abs(X1))),grid // În figura 4 s-a reprezentat spectrul de amplitudini iar în figura 5 se va reprezenta spectrul de fază calculat cu ajutorul funcţiei angle(). figure(5),plot(w,fftshift(angle(X1))),grid // În final putem reprezenta părţile reală şi imaginară ale TFD observăndu-se simetria, respectiv antisimetria acestora. figure(6) subplot(2,1,1),plot(w,fftshift(real(X1))),grid subplot(2,1,2),plot(w,fftshift(imag(X1))),grid

-4 -3 -2 -1 0 1 2 3 40

1

2

3

4

5

6

7

8

9

10

-4 -3 -2 -1 0 1 2 3 4-3

-2

-1

0

1

2

3

-4 -3 -2 -1 0 1 2 3 4-5

0

5

10

-4 -3 -2 -1 0 1 2 3 4-10

-5

0

5

10


36

Observaţii: • Reprezentarea spectrelor a fost făcută în funcţie de pulsaţia normată

[ , )ω π π∈ − . • Valorile calculate în X1=fft(x1,512) reprezintă spectrul calculat în

2kkN

ω π= , cu 1,,1,0 −= Nk … . Astfel, numărul de puncte N în care se

calculează fft determină rezoluţia spectrală. Cu cât N este mai mare cu atât aproximarea spectrului de frecvenţe continue )( ωjeX este mai bună.

• Dacă lungimea secvenţei ][nx este mai mică decăt N, secvenţa se completează cu zeroruri, spectrul calculat reprezentând convoluţia între spectrul secvenţei de lungime infinită şi spectrul funcţiei poartă de lungime egală cu a secvenţei date.

• Spectrul calculat cu funcţia fft poate fi reprezentat şi în funcţie de o indicele TFD: 1,,1,0 −= Nk … .

o frecvenţa normată: 2S

FfF

ωπ

= = , [ 0.5,0.5)f ∈ − .

f = -0.5:1/N:0.5-1/N; o frecvenţa nenormată SF f F= ⋅ [Hz]. [ / 2, / 2)Es EsF F F∈ − .

E3. Exerciţii: 1. a) Să se genereze secvenţa discretă [ ]s n obţinută prin eşantionarea cu frecvenţa de eşantionare 8kHzSF = a semnalului ( )0( ) sin 2πs t F t= de frecvenţă

0 500HzF = şi durată 40msMAXt = . Câte eşantioane are secvenţa discretă? b) Să se calculeze TFD a secvenţei în N=256 puncte cu funcţia fft. c) Să se reprezinte spectrul de amplitudini (modulul TFD) în funcţie de k=0:N-1 d) Să se reprezinte spectrul de amplitudini în pulsaţii normate [ , )ω π π∈ − e) Să se reprezinte spectrul de amplitudini în frecvenţe normate [ 0.5,0.5)f ∈ − f) Să se reprezinte spectrul de amplitudini în frecvenţe nenormate [Hz] g) Să se reprezinte cu subplot spectrul de amplitudini şi de fază în funcţie de frecvenţa normată. h) Să se reprezinte cu subplot partea reală şi partea imaginară în funcţie de frecvenţa normată. Fiecare spectru va fi reprezentat într-o figură separată (funcţia figure()). Graficul va fi completat cu funcţiile: grid, title şi xlabel. 2. Pentru programul de la exerciţiul anterior modificaţi pe rând următorii parametri şi explicaţi schimbările apărute în reprezentarea spectrului semnalului. a) Frecvenţa semnalului 0 200HzF = . b) Numărul de puncte al TFD N=1024 (pentru semnalul iniţial cu 0 500HzF = ).


37

2.2. Sisteme în timp discret Un sistem în timp discret transformă secvenţa de intrare ][nx într-o altă secvenţă de ieşire ][ny . Notând •S , operatorul sistemului, acesta este descris matematic prin relaţia:

][][ nxny S= (2.10) Reprezentarea simbolică a unui sistem în timp discret este următoarea: Pentru un sistem în timp discret, liniar şi invariant în timp (SDLIT) se numeşte funcţie pondere sau răspuns la impulsul unitate, secvenţa care se obţine la ieşirea sistemului dacă la intrare s-a aplicat impulsul unitate ][nδ :

][][ nnh δS= (2.11)

Astfel, un alt mod de a reprezenta simbolic un sistem în timp discret este: Răspunsul ][ny al unui SDLIT la orice secvenţă ][nx , poate fi determinat prin convoluţie dacă se cunoaşte răspunsul la impulsul unitate ][nh :

∑ −=∗=k

knhkxnhnxny ][][][][][ (2.12)

SDLIT pot fi reprezentate prin ecuaţii cu diferenţe finite cu coeficienţi constanţi, care dau legătura între secvenţa de intrare şi cea de ieşire:

∑∑==

−=−M

kk

N

kk knxbknya

00][][ (2.13)

Pentru analiza în frecvenţă a SDLIT reamintim formula pentru

transformata Z a unei semnal în timp discret ][nx :

∑ −=n

nznxzX ][)( (2.14)

•S][nx ][ny

][nh][nx ][ny


38

Pentru un SDLIT se numeşte funcţie de sistem sau funcţie de transfer )(zH , raportul între transformatele Z ale secvenţei de ieşire (răspunsul ][ny ) şi secvenţei de intrare (excitaţia ][nx ), reprezentând de fapt transformata Z a funcţiei pondere ][nh :

][)()(

][][)( nhZ

zXzY

nxZnyZzH === (2.15)

Ţinând cont de ecuaţia cu diferenţe finite şi de proprietatea de întârziere a transformatei Z ( )(][ zXzknx k−=−Z ), funcţia de sistem )(zH se mai poate scrie astfel:

∑

∑

=

−

=

−

+= N

k

kk

M

k

kk

za

zbzH

1

0

1)( (2.16)

în care s-a presupus că 10 =a . Rădăcinile polinomului de la numărător se numesc zerourile funcţiei de transfer iar rădăcinile polinomului de la numitor se numesc polii funcţiei de transfer. Pentru ca un sistem cauzal să fie stabil polii trebuie să fie situaţi în interiorul cercului de rază unitate (modulul lor să fie subunitar). 2.2.1. Răspunsul la impuls al unui SDLIT Funcţia MATLAB impz permite determinarea şi afişarea răspunsului la impuls (funcţia pondere) ][nh a unui SDLIT dacă se cunosc coeficienţii ka şi kb din ecuaţia cu diferenţe finite sau din expresia funcţiei de transfer )(zH . Sintaxe: [h,t] = impz(b,a) • vectorul b conţine coeficienţii kb ( b = [ b0, b1, ..., bM ] ) iar vectorul a conţine

coeficienţii ka ( a = [ 1, a1, a2, ..., aN ] ); se vor returna un vector coloană h care va conţine valorile eşantioanelor răspunsului la impuls al sistemului,

][nh , şi un vector coloană t ce va conţine momentele de pe axa timp (abscisa), alese în mod implicit, în care au fost calculate valorile eşantioanelor; parametrul de ieşire t poate să lipsească din sintaxă în cazul în care ne interesează doar răspunsul la impuls h.


39

[h,t] = impz(b,a,n) • se precizează numărul de puncte n în care se doreşte a fi calculat răspunsul la

impuls h; vectorul coloană t va conţine valorile acestor puncte (acestea vor fi 0, 1, 2, ..., n-1); parametrul de ieşire t poate să lipsească din sintaxă în cazul în care ne interesează doar răspunsul la impuls h.

[h,t] = impz(b,a,n,Fs) • aceleaşi considerente ca în sintaxa precedentă cu deosebirea că cele n valori

din vectorul t vor fi distanţate cu pasul 1/Fs (0, 1/Fs, 2/Fs, ..., (n-1)/Fs). [h,t] = impz(b,a,[],Fs) • aceleaşi considerente ca în sintaxa precedentă cu deosebirea că se alege în

mod implicit numărul de puncte în care se calculează răspunsul la impuls h. impz(b,a,...) • reprezintă grafic răspunsul la impuls calculat; punctele de suspensie au fost

introduse pentru a sugera faptul că se poate folosi oricare combinaţie a parametrilor de intrare din sintaxele precedente.

Exemple: Să se determine şi să se reprezinte grafic funcţia pondere a unui SDLIT definit prin: 1. ]2[6,0]1[6,0][3,0]1[9,0][ −+−+=−− nxnxnxnyny

2. 21

1

81,0)16/cos(8,115,01)( −−

−

+−+

=zz

zzHπ

b=[0.3,0.6,0.6]; a=[1,-0.9]; // s-au definit vectorii b şi a ce conţin valorile coeficienţilor kb şi ka din ecuaţia cu diferenţe finite. [h,t]=impz(b,a); // s-a calculat răspunsul la impuls al sistemului. size(h) → ans = 93 1 size(t) → ans = 93 1 // se verifică faptul că h şi t sunt vectori coloană (vezi sintaxa) cu 93 de elemente fiecare. Se poate verifica, tastând în fereastra de comenzi t urmat de enter, că valorile vectorului t sunt 0, 1, …, 92. Pentru reprezentarea grafică se poate proceda în două moduri:


40

impz(b,a),grid stem(t,h),grid // se va obţine acelaşi rezultat grafic. În mod asemănător procedăm şi pentru cel de al doilea exemplu: b=[1,0.5]; a=[1,-1.8*cos(pi/16),0.81]; h=impz(b,a); // vectorul coloană h va contine valorile eşantioanelor răspunsului la impuls al sistemului definit prin funcţia de transfer )(zH din exemplu. impz(b,a),grid 2.2.2. Răspunsul unui SDLIT la un semnal de intrare Funcţia MATLAB filter permite determinarea răspunsului ][ny al unui SDLIT dacă se cunosc coeficienţii ka şi kb din ecuaţia cu diferenţe finite sau din expresia funcţiei de transfer )(zH şi semnalul de intrare în sistem ][nx . Sintaxe: y = filter(b,a,x) • dacă x este un vector atunci se returnează un vector y de aceeaşi dimensiune

cu vectorul x; vectorul x conţine valorile semnalului de intrare în filtru (excitaţia); vectorul b conţine coeficienţii kb ( b = [ b0, b1, ..., bM ] ) iar vectorul a conţine coeficienţii ka ( a = [ 1, a1, a2, ..., aN ] ); dacă primul element din vectorul a este diferit de 1 atunci funcţia filter normează

0 10 20 30 40 50 60 70 80 900

0.2

0.4

0.6

0.8

1

1.2

1.4

0 10 20 30 40 50 60 70 80 90-1

0

1

2

3

4

5


41

coeficienţii ka ai sistemului la valoarea primului element din vectorul a (astfel primul element din vectorul a devine 1); vectorul obţinut y reprezintă răspunsul sistemului definit de coeficienţii din vectorii b şi a, dacă la intrare este aplicată secvenţa definită de vectorul x.

• dacă x este o matrice atunci se returnează o matrice y de aceeaşi dimensiune cu matricea x; funcţia filter va opera în acest caz pe coloane: coloana k din matricea y reprezintă răspunsul sistemului definit de coeficienţii din vectorii b şi a, dacă la intrare este aplicată coloana k din matricea x.

Exemple: Să se determine şi să se reprezinte grafic răspunsul unui sistem definit prin: 1. ]2[6,0]1[6,0][3,0]1[9,0][ −+−+=−− nxnxnxnyny

2. 21

1

81,0)16/cos(8,115,01)( −−

−

+−+

=zz

zzHπ

la semnalul de intrare ]10[][][1 −−= nununx , pentru 400 ≤≤ n . b=[0.3,0.6,0.6]; a=[1,-0.9]; // vectorii b şi a conţin valorile coeficienţilor kb şi ka . x=treapta(0,40,0)-treapta(0,40,10); // s-a definit vectorul x corespunzător secvenţei de intrare. y=filter(b,a,x); // s-a calculat răspunsul sistemului la secvenţa de intrare definită prin vectorul x. n=0:40; subplot(3,1,1),stem(n,x),grid,title(’x[n]’) subplot(3,1,2),impz(b,a),grid,title(’h[n]’) subplot(3,1,3),stem(n,y),grid,title(’y[n]’) În mod asemănător procedăm şi pentru cel de al doilea exemplu: b=[1,0.5]; a=[1,-1.8*cos(pi/16),0.81];

0 5 10 15 20 25 30 35 400

0.5

1x[n]

0 10 20 30 40 50 60 70 80 900

0.5

1

1.5h[n]

0 5 10 15 20 25 30 35 400

5

10y[n]


42

y=filter(b,a,x); subplot(3,1,1),stem(n,x),grid,title(’x[n]’) subplot(3,1,2),impz(b,a),grid,title(’h[n]’) subplot(3,1,3),stem(n,y),grid,title(’y[n]’) 2.2.3. Răspunsul în frecvenţă al SDLIT Funcţia MATLAB freqz permite determinarea răspunsului în frecvenţă al unui SDLIT dacă se cunosc coeficienţii ka şi kb din ecuaţia cu diferenţe finite sau din expresia funcţiei de transfer )(zH . Dacă este apelată fără parametrii de ieşire, cum se va vedea în sintaxă, această funcţie reprezintă grafic caracteristicile de amplitudine-frecvenţă şi fază-frecvenţă ale SDLIT respectiv. Sintaxe: [H,W] = freqz(b,a,n) • vectorul b conţine coeficienţii kb ( b = [ b0, b1, ..., bM ] ), vectorul a conţine

coeficienţii ka ( a = [ 1, a1, a2, ..., aN ] ), iar n reprezintă numărul de puncte în care se calculează răspunsul în frecvenţă H; vectorul W va conţine valorile acestor n puncte (valorile vor fi cuprinse între 0 şi π); este recomandat să se aleagă n putere a lui 2 (pentru a permite un calcul eficient folosind un algoritm FFT rapid); dacă n nu se specifică se alege în mod implicit 512.

[H,F] = freqz(b,a,n,Fs) • această sintaxă permite specificarea unei valori pentru frecvenţa de

eşantionare Fs (în Hz); vectorul F va conţine valorile celor n puncte în care se calculează răspunsul în frecvenţă H (în acest caz valorile acestor puncte vor fi cuprinse între 0 şi Fs/2);

[H,W] = freqz(b,a,n,’whole’)

0 5 10 15 20 25 30 35 400

0.5

1x[n]

0 10 20 30 40 50 60 70 80 90-5

0

5h[n]

0 5 10 15 20 25 30 35 40-20

0

20

40y[n]


43

• aceleaşi considerente ca în cazul primei sintaxe cu deosebirea că valorile celor n puncte de calcul, conţinute în vectorul W, vor fi cuprinse între 0 şi 2π; dacă n nu se specifică se alege în mod implicit 512.

[H,F] = freqz(b,a,n,’whole’,Fs) • aceleaşi considerente ca în cazul celei de a doua sintaxe cu deosebirea că

valorile celor n puncte de calcul, conţinute în vectorul F, vor fi cuprinse între 0 şi Fs; dacă n nu se specifică se alege în mod implicit 512.

H = freqz(b,a,W) • răspunsul în frecvenţă H se calculează la frecvenţele specificate în vectorul W;

valorile acestor frecvenţe trebuie să fie cuprinse între 0 şi 2π; dacă W nu se specifică se aleg în mod implicit 512 valori de frecvenţă.

H = freqz(b,a,F,Fs) • răspunsul în frecvenţă H se calculează la frecvenţele specificate în vectorul F;

valorile acestor frecvenţe trebuie să fie cuprinse între 0 şi Fs (frecvenţa de eşantionare în Hz).

freqz(b,a,...) • reprezintă grafic caracteristicile amplitudine-frecvenţă şi fază-frecvenţă ale

răspunsului în frecvenţă calculat; punctele de suspensie au fost introduse pentru a sugera faptul că se poate folosi oricare combinaţie a parametrilor de intrare din sintaxele precedente.

Exemple: Să se determine răspunsul în frecvenţă al SDLIT definite prin: 1. ]2[3,0]1[6,0][3,0]1[9,0][ −+−+=−+ nxnxnxnyny

2. 2

2

268.01634.0634.0)( −

−

−−

=z

zzH

Să se reprezinte grafic caracteristicile amplitudine-frecvenţă şi fază-frecvenţă ale răspunsului în frecvenţă calculat. b=[0.3,0.6,0.3]; a=[1,0.9]; [H,W]=freqz(b,a); // s-a calculat răspunsul în frecvenţă H în 512 puncte de frecvenţă cuprinse în intervalul ],0[ π . Reprezentările grafice cerute se pot realiza în două moduri: figure(1) subplot(2,1,1),plot(W,abs(H)),grid subplot(2,1,2),plot(W,angle(H)),grid


44

figure(2)freqz(b,a) În mod asemănător procedăm şi pentru cel de al doilea exemplu: b=[0.634,0,-0.634]; a=[1,0,-0.268]; [H,W]=freqz(b,a); freqz(b,a) 2.2.4. Diagrama poli-zerouri pentru funcţia de sistem a unui SDLIT Funcţia MATLAB zplane permite afişarea diagramei poli-zerouri în cazul funcţiei de sistem a unui SDLIT dacă se cunosc valorile polilor şi zerourilor sau dacă se cunosc doar coeficienţii ka şi kb din ecuaţia cu diferenţe finite sau din expresia funcţiei de transfer )(zH . Sintaxe: zplane(z,p) • dacă z şi p sunt doi vectori coloană ce conţin valorile zerourilor şi respectiv

polilor funcţiei de transfer H z( ) atunci se va afişa diagrama poli-zerouri, marcând zerourile cu semnul ‘o’ iar polii cu semnul ‘x’; dacă există poli sau zerouri multiple, acestea vor avea înscris, lângă semnul respectiv şi ordinul de multiplicitate.

• dacă z şi p sunt două matrice afişarea diagramei poli-zerouri se va face pentru fiecare coloană în parte cu culori diferite.

zplane(b,a)

0 0.5 1 1.5 2 2.5 3 3.50

0.2

0.4

0.6

0.8

0 0.5 1 1.5 2 2.5 3 3.5-4

-3

-2

-1

0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-200

-150

-100

-50

0

Normalized frequency (Nyquist == 1)

Pha

se (d

egre

es)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-80

-60

-40

-20

0

Normalized frequency (Nyquist == 1)

Mag

nitu

de R

espo

nse

(dB

)


45

• dacă b şi a sunt doi vectori linie ce conţin valorile coeficienţilor kb şi ka atunci se va afişa diagrama poli-zerouri a funcţiei de transfer H z( ) (calculându-se rădăcinile polinoamelor de la numărătorul şi numitorul funcţiei de sistem).

Reprezentarea polilor şi zerourilor unei funcţii de sistem (diagrama poli-

zerouri) se face în planul Z, în raport cu cercul de rază unitate, având pe abscisă partea reală şi pe ordonată partea imaginară. Fiind privite deci ca numere complexe, valorile respective pot fi exprimate în formă polară sau formă carteziană. Orice număr complex z poate fi exprimat în formă polară astfel: z z e j= ⋅ ϕ (2.17) în care: - z = modulul numărului complex z ; - ϕ = argumentul (faza) numărului complex z ;

În formă carteziană numărul complex z se exprimă sub forma: z real z j imag z= + ⋅( ) ( ) (2.18)

Reprezentarea în planul Z a unui număr complez z

Atenţie: În cazul în care dispunem de valorile coeficienţilor kb şi ka şi dorim să determinăm valorile polilor şi zerourilor funcţiei de sistem respective se poate utiliza funcţia MATLAB roots, care calculează rădăcinile unui polinom dacă sunt precizaţi coeficienţii acestuia.

z

|z|

realz

imagz

Re

Im

1

1

-1

-1 0ϕ


46

Dacă se doreşte determinarea valorilor coeficienţilor kb şi ka , având valorile polilor şi zerourilor funcţiei de sistem, putem utiliza funcţia MATLAB poly, care calculează coeficienţii unui polinom dacă sunt precizate rădăcinile acestuia.

Verificaţi sintaxele celor două funcţii folosind comanda help. Exemple:

1. Funcţia de transfer a unui SDLIT are un zerou de valoare 3

21 πj

er = . Ştiind că

această funcţie are zerouri şi în ∗r , r1 , ∗r

1 şi are doi poli în 32

31 πj

eq = şi ∗q să

se reprezinte diagrama poli-zerouri şi să se scrie forma nefactorizată a funcţiei de transfer (să se găsească valorile coeficienţilor kb şi ka ). r=1/2*exp(j*pi/3); q=1/3*exp(j*2*pi/3); z=[r;conj(r);1/r;1/conj(r)]; p=[q;conj(q)]; // s-au definit vectorii coloană z şi p (vezi sintaxa) ce conţin valorile zerourilor şi respectiv polilor funcţiei de transfer. zplane(z,p) // s-a reprezentat diagrama poli-zerouri. b=poly(z) → b =

1.0000 -2.5000 5.2500 -2.5000 1.0000 a=poly(p) → a = 1.0000 0.3333 0.1111 // s-au calculat valorile coeficienţilor kb şi ka ai funcţiei de sistem. Forma nefactorizată a acestei funcţii va fi:

21

4321

3333,015,225,55,21)( −−

−−−−

+−+−+−

=zz

zzzzzH

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

-1.5

-1

-0.5

0

0.5

1

1.5

Real part

Imag

inar

y pa

rt


47

2. Să se reprezinte diagramele poli-zerouri pentru sistemele în timp discret definite prin:

a. funcţia de transfer: 21

1

81,0)16/cos(8,115,01)( −−

−

+−+

=zz

zzHπ

b. funcţia de transfer: 21

21

81,0)16/cos(8,11)5,01()( −−

−

+−+

=zz

zzHπ

c. ecuaţia cu diferenţe finite: ]2[6,0]1[6,0][3,0]4[1.0]3[1.0]1[1.0][ −+−+=−+−+−+ nxnxnxnynynyny b=[1,0.5]; a=[1,-1.8*cos(pi/16),0.81]; zplane(b,a) b=[1,1,0.25]; zplane(b,a) // se observă ordinul de multiplicitate ale zeroului. b=[0.3,0.6,0.6]; a=[1,0.1,0,0.1,0.1]; zplane(b,a)

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Real part

Imag

inar

y pa

rt

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Real part

Imag

inar

y pa

rt

2

-1 -0.5 0 0.5 1

-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Real part

Imag

inar

y pa

rt


48

E4. Exerciţii: Se dau următoarele SDLIT definite prin: 1. ]4[125,0]3[5,0]2[81,0]1[27,1][][ −+−−−+−−= nxnxnxnxnxny 2. ][]1[9,0][ nxnyny =−+ 3. =−+−+−+ ]3[3,0]2[52,0]1[13,0][ nynynyny

]3[16,0]2[48,0]1[48,0][16,0 −−−+−−= nxnxnxnx 4. ]2[3,0]1[6,0][3,0]2[9,0][ −+−+=−+ nxnxnxnyny

5. 4321

321

024,04,064,08,01075,0125,05,01)( −−−−

−−−

+−+−−+−

=zzzz

zzzzH

6. 31 44,077,011)( −− +−

=zz

zH

7. 74321 3,0125,05,081,027,11)( −−−−− −+−+−= zzzzzzH La intrarea acestor sisteme se pot aplica oricare din semnalele de intrare de mai jos: I. ][][1 nnx δ= pentru 400 ≤≤ n II. ][][2 nunx = pentru 400 ≤≤ n

III. 3,0 10

[ ]20 ,11 20

n nx n

n n≤ ≤⎧

= ⎨ − ≤ ≤⎩

IV. 4[ ] sin5

nx n π⎛ ⎞= ⎜ ⎟⎝ ⎠

pentru 200 ≤≤ n

a) Să se determine răspunsul acestor sisteme la semnale de intrare I–IV şi să se

reprezinte grafic în domeniul timp (folosind subplot) semnalul de intrare, funcţia pondere a sistemului şi semnalul de ieşire.

b) Pentru fiecare caz analizat să se reprezinte grafic în domeniul frecvenţă (folosind subplot) semnalul de intrare, funcţia de transfer a sistemului şi semnalul de ieşire.

c) Să se reprezinte diagramele poli-zerouri asociate sistemelor 1–7. Tema de casă 1. Pentru semnalele primite ca tema la lucrarea 1 să se calculeze transformata Fourier discretă folosind funcţia fft. Obs. Semnalul 3 de tip dreptunghiular multinivel aleator se va înlocui cu

semnal multinivel periodic (obţinut prin repetarea aceleiaşi succesiuni de niveluri). Semnalul 8 de tip sinusoidal de frecvenţă variabilă se va înlocui cu un semnal sinusoidal de perioadă 2.5 ms şi amplitudine 3 limitat superior la valoarea maximă 2.


49

Lungimea N a TFD va fi aleasă o putere a lui 2 (128, 256, 512 sau 1024) mai mică decât lungimea secvenţei. Să se reprezinte grafic: 1) Spectrul de amplitudini

a) În funcţie de indicele k al TFD calculată. k=0…N-1. b) În funcţie de pulsaţii normate [-π, π]. c) În funcţie de frecvenţa normată [-0.5, 0.5]. d) În funcţie de frecvenţa nenormată [-Fs/2 – Fs/2].

2) Spectrul de amplitudini în dB în funcţie de frecvenţe normate. 3) Spectrul de faze în funcţie de frecvenţe normate. Să se reia reprezentarea spectrelor pentru lungimea N a TFD aleasă astfel ca să îndeplinească condiţia ,SN k T F k= ⋅ ⋅ ∈ (N să fie de ordinul sutelor, mai mic decât lungimea secvenţei). Explicaţi diferenţele care apar între spectre în cele două situaţii (deşi se reprezintă spectrul aceluiaşi semnal). 2.

a) Să se genereze următorul semnal discret: [ ] [ ] [ ]x n s n v n= + unde • v[n] este un semnal aleator cu distribuţie normală (gaussiană). • s[n] este un semnal sinusoidal de amplitudine 2 şi frecvenţa F0 aleasă

aleator în intervalul [FS/10, FS/3]. Frecvenţa de eşantionare FS = 8kHz. Să se reprezinte semnalul x[n] în domeniul timp şi spectrul semnalului în frecvenţe nenormate şi în pulsaţii normate. Identificaţi din spectru pulsaţia ω0 corespunzătoare frecvenţei F0 a componentei sinusoidale de la intrare. b) Se consideră sistemul discret cu funcţia de transfer 1 2

0 1 2( )H z b b z b z− −= + + având atenuare infinită în ω0 determinat la punctul a). Atenuare infinită înseamnă câştig 0 adică 0( ) 0jH e ω = deci 0

1jz e ω= este zerou al lui H(z).

• Determinaţi zerourile astfel încât funcţia H(z) să aibă coeficienţi reali şi scrieţi expresiile teoretice ale lui H(z) şi ( )jH e ω .

• Reprezentaţi grafic diagrama poli-zerouri. • Reprezentaţi grafic răspunsul în frecvenţă al sistemului. • Reprezentaţi grafic răspunsul la impuls al sistemului. • Calculaţi ieşirea sistemului dacă la intrare se aplică semnalul x[n]. • Reprezentaţi spectrul semnalului de la ieşire în frecvenţe nenormate şi în

pulsaţii normate.


50

c) Se consideră sistemul discret cu funcţia de transfer 1 2

0 1 21 2

1 2

( )1b b z b zH z

a z a z

− −

− −

+ +=

+ +

având zerourile zk determinate la punctul b) şi polii de forma 1 ρ kp z= , 0.7<ρ<1 (polii au acelaşi argument ca zerorurile dar sunt de modul subunitar).

• Determinaţi expresia teoretică a coeficienţilor şi scrieţi expresiile H(z) şi ( )jH e ω .

• Calculaţi toate zerourile şi toţi polii funcţiei H(z) şi reprezentaţi grafic diagrama poli/zerouri. Se va alege 0.7<ρ<1.

• Reprezentaţi grafic răspunsul în frecvenţă al sistemului. Cum se modifică răspunsul în frecvenţă dacă polii se apropie de cercul de rază unitate (ρ→1)?

• Reprezentaţi grafic răspunsul la impuls al sistemului. • Calculaţi ieşirea sistemului dacă la intrare se aplică semnalul x[n]. • Reprezentaţi spectrul semnalului de la ieşire în frecvenţe nenormate şi în

pulsaţii normate. 3. Reluaţi problema 2 înlocuind semnalul v[n] cu un semnal vocal. (se va folosi exemplul din arhiva demo1.zip aflată pe site-ul laboratorului).

LUCRAREA 2 SEMNALE ŞI SISTEME ÎN TIMP DISCRET · fft. Denumirea sa reprezintă prescurtarea de la...

Documents

Transcript of LUCRAREA 2 SEMNALE ŞI SISTEME ÎN TIMP DISCRET · fft. Denumirea sa reprezintă prescurtarea de la...