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.
|
