Este important de remarcat faptul că actuala versiune a compilatorului mcc acceptă functia input cu parametri variabile din spatiul de lucru. Aceasta permite cresterea flexibilitătii programelor, deoarece, dacă nu se întâmpla acest lucru, era nevoie de o recompilare la fiecare schimbare a parametrilor.
Calculul răspunsului la intrare tre 13513m127n aptă a unui SLIT continual de ordinul 2
Deoarece nu se pot folosi puternicele functii orientate pe obiect din Control Toolbox, este necesar să aplicăm pe cât posibil functii din biblioteca matematică a MATLAB.
Astfel, asa cum am arătat într-un exemplu anterior, nu putem aplica functii precum step, impulse, lsim etc. dacă acestea au ca parametri obiecte de tip sistem.
În cazul în care putem folosi ca parametri variabile vectoriale în loc de obiecte, compilatorul va putea functiona normal.
Cazul 1: Dacă dorim să calculăm răspunsul la intrare treaptă al unui sistem de ordinul 2, continual, liniar si invariant în timp, vom putea să aplicăm expresia cunoscută, calculată analitic:
Programul MATLAB este:
function []=test2();
iar='y';
while iar=='y'
T=input('timp simulare [sec]= ');
pas=input('pas de simulare= ');
omega=input('omega= ');
csi=input('csi= ');
ta=input('temporizare: ');
t=0:pas:T;
fi=atan(sqrt(1-csi^2)/csi);
y=1-exp(-csi*omega.*t)/(sqrt(1-csi^2)).*sin(omega*sqrt(1-csi^2).*t+fi);
plot(t,y);grid;
title(sprintf('Raspuns intrare treapta pentru\n
omegan=%d, csi=%d',omega,csi));
pause(ta);
iar=input('Continuam? [y/n] ','s');
end;
Această
functie permite ca, fără să parăsim programul, să putem varia atât
parametrii sistemului (factor de amortizare 
sau pulsatia
proprie 
) cât si parametrii de simulare (pas de integrare, timp
de simulare). Se poate astfel vizualiza răspunsul la intrare treaptă pentru
orice sistem de ordinul 2 si pe orice interval de timp, afisându-se,
o dată cu răspunsul sistemului, si parametrii sistemului.
Temporizarea a fost introdusă pentru ca programul să aibă timp să ploteze graficul cerut, altfel plotarea ar fi fost întârziată de cererea permanentă a programului dacă se reia sau nu.
Cazul 2: Calculul răspunsului indicial al SLIT monovariabile continuale
Dacă dorim să extindem problema anterioară, vom considera cazul mai general al sistemelor monovariabile, nu neapărat de ordinul 2. Vom folosi metoda răspunsului general al SLIT continuale:
Programul va implementa calculul atât al răspunsului indicial, cât si al celui ponderal, pornind de la formula de mai sus. Răspunsul indicial se obtine din cel ponderal prin integrare, deci se va folosi functia cumsum.
function []=test3();
% raspuns general la intrare treapta
% si impuls al sistemelor monovariabile
iar='y';
while iar=='y'
tf=input('timpul de simulare: ');
pas=input('pasul de integrare: ');
t=0:pas:tf;
% se introduce sistemul sub forma H(s)=num(s)/den(s)
num=input('numaratorul: ');
den=input('numitorul: ');
[a,b,c,d]=tf2ss(num,den);
for k=1:length(t)
y(:,k)=c*expm(a*t(k))*b;
end;
z=pas*cumsum(y'); % se calculeaza integrala
figure(1);plot(t,z);grid;
title('Raspuns intrare treapta monovariabil');
figure(2);plot(t,y);grid;
title(sprintf('Raspuns impuls monovariabil'));
pause(5);
iar=input('Continuam? [y/n] ','s');
end;
Dacă vom compila programul de mai sus cu comanda
>> mcc -m -B sgl test3.m
vom obtine un executabil, test3.exe care va putea fi lansat în executie din mediul DOS (Command Prompt). Se obtin rezultate ca în figura de mai jos.
Cazul 3: Calculul răspunsului indicial si ponderal în cazul SLIT continuale multivariabile
În acest caz se preferă reprezentarea de stare a sistemului, deci este nevoie de o interfată prin care se cer valorile efective ale matricilor A, B, C, D.
Programul MATLAB este extrem de simplu raportat la volumul de calcule efectuate:
function []=test3();
% raspuns general indicial si ponderal multivariabil
iar='y';
while iar=='y'
tf=input('timpul de simulare: ');
pas=input('pasul de integrare: ');
t=0:pas:tf;
% se introduce sistemul sub forma H(s)=(A,B,C,D)
a=input('matricea A: ');
b=input('matricea B: ');
c=input('matricea C: ');
d=input('matricea D: ');
z=zeros(size(c*b));
for k=1:length(t)
y(:,k)=c*expm(a*t(k))*b*ones(size(b,2),1);
end;
z=pas*cumsum(y');
plot(t,z);grid;
title('Raspuns intrare treapta multivariabil');
pause(5);
iar=input('Continuam? [y/n] ','s');
end;
Mai jos se prezintă un exemplu de rulare a executabilului test3.exe pentru cazul unui sistem cu două intrări si cu două iesiri:
Calculul răspunsului general al SLIT continuale
Cazul 1: Sisteme monovariabile
În cazul sistemelor monovariabile este mai simplu să introducem sistemul al cărui răspuns la impuls dorim să-l calculăm sub formă de functie de transfer. Fată de răspunsul indicial, diferenta în programele MATLAB constă în nefolosirea functiei cumsum (suma cumulativă) care implementa integrala de convolutie, ci a valorii la momentul t a cantitătii de sub integrală.
Cazul 2: Sisteme multivariabile
În acest caz, vom folosi extraordinara putere de calcul a MATLAB pentru a scrie acest program general:
function []=test4();
% raspuns general multivariabil
iar='y';
while iar=='y'
tf=input('timpul de simulare: ');
pas=input('pasul de integrare: ');
t=0:pas:tf;
% se introduce sistemul sub forma H(s)=(A,B,C,D)
a=input('matricea A: ');
b=input('matricea B: ');
c=input('matricea C: ');
d=input('matricea D: ');
in=input('intrarea nr. :');
% se genereaza intrarea u
u=zeros(size(b,2),length(t));
% la intrarea respectiva se aplica 1+sin(5t)
u(in,:)=1+sin(5*t);
z=[];
% se calculeaza integrala de convolutie
for n=1:length(t)
y=zeros(size(c,1),1);
for k=1:n
y=y+pas*c*expm(a*t(k))*b*u(:,n-k+1);
end;
z(:,n)=y; % valoarea integralei de convolutie la fiecare moment de timp
end;
plot(t,z);grid;
title(sprintf('Raspuns multivariabil la intrarea %d',in));
pause(5);
iar=input('Continuam? [y/n] ','s');
end;
Programul implementează formula generală de calcul prezentată anterior.
Se observă
cum programul a fost compilat pentru o valoare de forma 
aplicată pe una dintre
intrările sistemului. Mai jos se prezintă iesirile pentru semnalul de mai
sus aplicat la intrarea nr.1 si apoi la cea de-a doua, pentru sistemul cu
reprezentarea de stare anterioară.
În cazul în care dorim să extindem flexibilitatea programului, adică să putem genera ce intrare dorim, o posibilitate ar fi următorul cod MATLAB:
s=input('expresia intrarii f(t)= ','s');
% se genereaza intrarea u
u=zeros(size(b,2),length(t));u(in,:)=eval(s);
Dacă vom compila acest program modificat, nu vom primi nici o avertizare, dar la rulare, în momentul în care se apelează functia eval vom obtine
expresia intrarii f(t)= 1+sin(10*t)
ERROR: Reference to unknown function or variable 't' while evaluating expression
EXITING
Se observă că trebuie să evităm folosirea acestei functii.
Trasarea caracteristicilor de frecventă
Semnificatie imediată au caracteristicile de frecventă pentru SLIT monovariabile. Pentru a calcula aceste caracteristici nu putem beneficia de functiile orientate pe obiect disponibile în MATLAB, asa că vom proceda la apelarea unor functii matematice.
Pentru a calcula expresiile de mai sus, vom apela la expresia functiei de transfer în functie de reprezentarea de stare.
Programul MATLAB care implementează relatiile de mai sus este:
function []=test5();
num=input('numaratorul: ');
den=input('numitorul: ');
[a,b,c,d]=tf2ss(num,den);
wmin=input('wmin (10^wmin): ');
wmax=input('wmax (10^wmax): ');
w=logspace(wmin,wmax,400);
for k=1:length(w)
mag(k)=abs(c*inv((j*w(k)*eye(size(a))-a))*b+d);
phase(k)=angle(c*inv((j*w(k)*eye(size(a))-a))*b+d);
end;
subplot(211);semilogx(w,20*log(mag));grid;
title('Caracteristica amplitudine-pulsatie');
subplot(212);semilogx(w,phase);grid;title('Caracteristica faza-pulsatie');
Ideea de utilizare a MATLAB este că acesta scurtează spectaculos timpul de realizare, depanare si testare a programelor. Cu numai câteva linii de program realizăm functiuni care ar necesita mii de linii de program în C, spre exemplu.
În cele ce urmează, dorim să realizăm un program care să combine câteva dintre functiile construite mai înainte. Spre exemplu, dorim să vizualizăm atât răspunsul indicial si ponderal al unui sistem monovariabil, cât si caracteristicile de frecventă.
Un program MATLAB care implementează cerintele de mai sus ar putea fi:
function []=test6();
% calculul raspusurilor in timp
% si al caracteristicilor de frecventa
num=input('numaratorul: ');
den=input('numitorul: ');
[a,b,c,d]=tf2ss(num,den);
rasp_timp(a,b,c,d);
figure(3);
carbode(a,b,c,d);
în care functiile rasp_timp.m respectiv carbode.m pot avea forma:
function []=rasp_timp(a,b,c,d);
% raspuns indicial si ponderal multivariabil
tf=input('timpul de simulare: ');
pas=input('pasul de integrare: ');
t=0:pas:tf;
z=zeros(size(c*b));
for k=1:length(t)
y(:,k)=c*expm(a*t(k))*b;
end;
z=pas*cumsum(y');
figure(1);plot(t,z);grid;
title('Raspuns indicial');
figure(2);plot(t,y);grid;
title('Raspuns ponderal');
function []=carbode(a,b,c,d);
% trasarea caracteristicilor bode
wmin=input('wmin (10^wmin): ');
wmax=input('wmax (10^wmax): ');
w=logspace(wmin,wmax,400);
for k=1:length(w)
mag(k)=abs(c*inv((j*w(k)*eye(size(a))-a))*b+d);
phase(k)=angle(c*inv((j*w(k)*eye(size(a))-a))*b+d);
end;
subplot(211);semilogx(w,20*log(mag));grid;title('Caracteristica amplitudine-pulsatie');
subplot(212);semilogx(w,phase);grid;title('Caracteristica faza-pulsatie');
Se observă din functiile de mai sus că ele au ca argumente de intrare cvadruplul (a, b, c, d) si nu au nici un argument de iesire. Acestea realizează operatiunile de vizualizare a răspunsurilor în timp si în frecventă în interiorul acestor functii (vizualizările pentru un exemplu sunt prezentate mai jos).
Pentru a putea creste flexibilitatea programării, putem concepe functiile rasp_timp respectiv carbode ca functii ce implementează numai mecanismul de calcul efectiv, lăsând în sarcina programului principal operatiunile de vizualizare.
În acest caz functiile vor trebui să returneze toti parametrii necesari, deci vor trebui să dispună si de parametri de iesire.
function []=test6();
% calculul raspusurilor in timp
% si al caracteristicilor de frecventa
num=input('numaratorul: ');
den=input('numitorul: ');
[a,b,c,d]=tf2ss(num,den);
% apelul functiei de timp
[t,y,z]=rasp_timp(a,b,c,d);
figure(1);plot(t,z);grid;
title('Raspuns indicial');
figure(2);plot(t,y);grid;
title('Raspuns ponderal');
% apelul functiei de frecventa
[w,mag,phase]=carbode(a,b,c,d);
figure(3);
subplot(211);semilogx(w,20*log(mag));grid;
title('Caracteristica amplitudine-pulsatie');
subplot(212);semilogx(w,phase);grid;title('Caracteristica faza-pulsatie');
function [t,y,z]=rasp_timp(a,b,c,d);
% raspuns indicial si ponderal multivariabil
tf=input('timpul de simulare: ');
pas=input('pasul de integrare: ');
t=0:pas:tf;
z=zeros(size(c*b));
for k=1:length(t)
y(:,k)=c*expm(a*t(k))*b;
end;
z=pas*cumsum(y');
function [w,mag,phase]=carbode(a,b,c,d);
% trasarea caracteristicilor bode
wmin=input('wmin (10^wmin): ');
wmax=input('wmax (10^wmax): ');
w=logspace(wmin,wmax,400);
for k=1:length(w)
mag(k)=abs(c*inv((j*w(k)*eye(size(a))-a))*b+d);
phase(k)=angle(c*inv((j*w(k)*eye(size(a))-a))*b+d);
end;
5.2. Programe de timp real folosind placa de sunet
Să presupunem că trebuie să verificăm că frecventa fundamentală a unui diapazon este de 440 Hz. Pentru aceasta este nevoie să folosim un microfon si o placă de sunet pentru a achizitiona date la nivel de sunet. Se foloseste apoi transformata Fourier rapidă (FFT) asupra datelor achizitionate pentru a afla spectrul de frecventă al semnalului sonor achizitionat. Pasii care trebuie parcursi sunt prezentati mai jos.
Pentru exemplul prezentat mai sus, vom achizitiona o secundă de semnal sonor pe unul dintre canalele plăcii de sunet. Placa de sunet poate achizitiona sunet pe două canale (sunet stereofonic) sau pe un singur canal (monfonic). Deoarece diapazonul oscilează pe o frecventă de 440 Hz, se poate configura placa de sunet la cea mai mică frecventă de esantionare, anume 8000 Hz. Chiar si la această cea mai mică frecventă de esantionare nu putem avea un fenomen de aliasing, doarece frecventa Nyquist corespunzătoare este de 4000 Hz, mult peste frecventa semnalului sonor achizitionat.
Se plasează diapazonul oscilând aproape de microfon si achizitia de semnal se validează (se trigherează, se porneste) folosind un trigger manual (vezi daqdoc4_1.m).
1. Crearea unui obiect dispozitiv - Crează un obiect de forma analog input AI pentru o placă de sunet. Adaptoarele instalate precum si identificatoarele hardware se pot afla folosind functia daqhwinfo.
AI = analoginput('winsound');
2. Se adaugă canalele dispozitivului AI, unul în cazul nostru.
chan = addchannel(AI,1);
3. Se configurează valorile corespunzătoare - se asignează valorile de bază, esentiale pentru placa de achizitia si se crează variabilele de achizitie precum si frecventa de esantionare Fs care se foloseste ulterior pentru analiza spectrală. Această valoare se poate culege direct de la placa de sunet deoarace aceasta functionează numai pe valori predefinite de frecvente care pot diferi de valorile pe care le setăm noi.
duration = 1; %1 second acquisition
set(AI,'SampleRate',8000)
ActualRate = get(AI,'SampleRate');
set(AI,'SamplesPerTrigger',duration*ActualRate)
set(AI,'TriggerType','Manual')
blocksize = get(AI,'SamplesPerTrigger');
Fs = ActualRate;
4. Se achizitionează datele - Comanda Start AI, declansează manual o achizitie (manual trigger) si se extrag datele din motorul de achizitie. Înainte de extragerea datelor, este necesar ca sunetul emis de diapazon să fie achizitionat de placa de sunet, de unde vor fi ulterior preluate.
start(AI)
trigger(AI)
data = getdata(AI);
5. stergerea din memorie - atunci când nu mai avem nevoie de AI, îl vom scoate din memorie si din spatiul de lucru (workspace) al MATLAB.
delete(AI)
clear AI
Pentru acest exemplu, analiza constă în găsirea componentelor de frecventă ale diapazonului si afisarea rezultatelor. Pentru aceasta a fost creată functia daqdocfft. Această functie calculează FFT ale esantioanelor de sunet pe baza datelor initiale SampleRate si SamplesPerTrigger.
[f,mag] = daqdocfft(data,Fs,blocksize);
daqdocfft returnează frecventa si amplitudinea datelor care pot fi afisate ulterior.
function [f,mag] = daqdocfft(data,Fs,blocksize)
% [F,MAG]=DAQDOCFFT(X,FS,BLOCKSIZE) calculates the FFT of X
% using sampling frequency FS and the SamplesPerTrigger
% provided in BLOCKSIZE
xfft = abs(fft(data));
mag = 20*log10(xfft);
mag = mag(1:blocksize/2);
f = (0:length(mag)-1)*Fs/blocksize;
f = f(:);
Rezultatele pot fi afisate cu secventa de mai jos.
plot(f,mag)
grid on
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
title('Frequency Components of Tuning Fork')
Un mod simplu de a calcula frecventa fundamentlă este
[ymax,maxindex]= max(mag);
maxindex
maxindex = 441
Această sectiune ilustrează modul în care putem să generăm un sunet folosind iesirea plăcii de sunet a calculatorului. Pasii de bază care trebuie parcursi în general sunt subliniati în continuare.
Se instalează si se conectează placa de sunet, împreună cu un dispozitiv de iesire, cum ar fi boxele în cazul nostru.
Se configurează sesiunea de iesire de date. Aceasta înseamnă crearea unui obiect dispozitiv, adăugarea de canale, setarea parametrilor de bază si utilizarea functiilor specifice pentru îndeplinirea obiectivului.
În exemplul de mai jos, un semnal sinusoidal se generează în MATLAB si apoi este trimis către convertorul Numeric/Analogic al plăcii de sunet care este conectată pe iesire la boxe.
(Pentru rularea exemplului se poate rula daqdoc6_1.m de la linia de comandă MATLAB.
1. Crearea unui obiect specific plăcii de sunet. - se crează un obiect de forma analog output AO pentru placa de sunet. Adaptoarele precum si componentele hardware se pot afla cu functia daqhwinfo.
AO = analogoutput('winsound');
2. Se adaugă canalele de iesire pentru AO.
chan = addchannel(AO,1);
3. Se configurează valorile parametrilor de bază ai interfetei de iesire de sunet. Se defineste o durată de 4 secunde, se setează parametrii plăcii, se generează datele si se transmit la dispozitivul de iesire.
duration = 4;
set(AO,'SampleRate',8000)
set(AO,'TriggerType','Manual')
ActualRate = get(AO,'SampleRate');
len = ActualRate*duration;
data = sin(linspace(0,2*pi,len))';
putdata(AO,data)
4. Se transmit datele la iesire. Se startează operatia de iesire, se triggerează procesul, se asteaptă până când sunetul a fost transmis si dispozitivul semnalizează off.
start(AO)
trigger(AO)
while strcmp(AO.Running,'On')
end
5. Se curătă memoria si spatiul de lucru MATLAB.
delete(AO)
clear AO
Din nefericire, după cum am mentionat mai sus, aceste fisiere nu se pot compila deocamdată cu ajutorul mcc deoarece compilatorul nu este capabil încă să compileze functii orientate pe obiect.
|
