Skip to content

Instantly share code, notes, and snippets.

@dugagjin
Last active January 9, 2016 22:08
Show Gist options
  • Save dugagjin/1d59942dce90b9d6cd0d to your computer and use it in GitHub Desktop.
Save dugagjin/1d59942dce90b9d6cd0d to your computer and use it in GitHub Desktop.
lineaire tijdsinvariante systemen
% Maak de FRF van een eerste ordesysteem, gegeven de transferfunctie: H(jw)
% = A/1+jw*Tau. Waarbij Tau = 0.5 en A = 1. Maak een Bodeplot, duid het -3 dB punt
% aan door middel van een rood cirkel.
clear all;clf;
A = 1;
Tau = 1/2;
w = linspace(0,100,10000);
s = j * w;
H = A ./ (Tau*s+1);
figure(1);
subplot(2,1,1);
semilogx(w, db(H),1/Tau,-3,'or'); % 1/Tau,-3 om de -3dB aan te duiden voor dB grafiek.
subplot(2,1,2)
semilogx(w, angle(H),1/Tau,atan(-1),'ro'); % 1/Tau,atan(-1) om -3dB punt te hebben voor hoek grafiek.
% Om verschil te zien van twee verschillende manieren van doen:
figure(2);
NUM = [A];
DEN = [Tau 1];
SYS = tf(NUM,DEN);
bode(SYS);
% Maak een FRF van een tweede ordesysteem, gegeven de transferfunctie:
% H(jw) = wn²/(-w²+2E*wn*j*w+wn²) waarbij wn = 10 en E = 0.1. Maak een bode
% plot, en duid de resonantiefrquentie aan door middel van een rode lijn
% (voor amplitude) en rood cirkel (voor fase).
clear all;clf;
wn = 10;
E = 0.1;
N = 1000;
w = linspace(0,100,N);
% amplitude in functie van rad/s
figure(1);
NUM = [wn.^2];
DEN = [-1 2*E*wn*1i wn.^2];
SYS = tf(NUM,DEN);
polen = pole(SYS);
H1 = wn.^2./(-w.^2+2*E*wn*1i*w+wn.^2);
subplot(2,2,1);
semilogx(w,db(H1),[real(polen(1)) real(polen(1))],[-100 100],'r');
% faseplot in functie van rad/s
% H2 voor resonantiefreq (rode cirkel)
subplot(2,2,2);
H2 = wn^2./(-real(polen(1))^2+2*E*wn*1i*real(polen(1))+wn^2);
semilogx(w,angle(H1),real(polen(1)),angle(H2),'ro');
% ligging van de polen (kruisje) en nullen (cirkel) plotten
subplot(2,2,3);
plot(-imag(polen(1)),real(polen(1)),'xb',-imag(polen(2)),real(polen(2)),'xb');
% limieten zetten anders zien we niet deftig wat er geplot word.
xlim([-15 5]);
ylim([-25 25]);
% Maak tenslotte ook de impulsresponsfunctie op basis van de FRF, en men plot deze. Laat dan
% de wn en de E varieren en observeer hoe de FRF en de polenligging verandert.
H1=real(ifft(H1,N));
H1 = H1(1:1000);
subplot(2,2,4);
plot(w,H1);
% Construeer nu een FRF op basis van de volgende polen en nullen:
% polen: p = -0.1+5j, -0.1-5j, 3
% nullen: z = -2
% Maak opnieuw een Bodeplot, de polen nullenligging en de
% impulsreponsfunctie. Varieer nu de ligging van de polen en observeer het
% resultaat. Wanneer heb je te maken met een onstabiel systeem? Hoe zie je
% dat de FRF(FrequencyReponsFunctie, aan de IRF (ImpulsResponsFunctie),
% en aan de polen-nullen ligging?
clear all;clf;
figure(1);
Z = -1; % nul op -2
K = 1; % versterking = 1 => dus geen versterking
p = [-0.1+5*j,-0.1-5*j,-3];
% bode plot
sys = zpk(Z, p, K); % ludiek manier om bode te maken
subplot(2,2,1);
bode(sys);
% ligging van de polen (kruisje)
subplot(2,2,3);
plot (real(Z), -imag(Z),'bo', real(p), -imag(p), 'xb');
% limieten zetten anders zien we niet deftig wat er geplot word.
xlim([-10 10]);
ylim([-10 10]);
% weet niet of juist vanaf hier maar normaal wel bro.
% Maak tenslotte ook de impulsresponsfunctie op basis van de FRF, en plot deze.
w = 0.1:0.1:100;
s = j*w;
H = (s-Z)./((s- p(1)).*(s- p(2)).*(s - p(3)));
N = 1000;
h = 2*real(ifft(H,2*N));
h = h(1:N);
subplot(2,2,2);
plot(h);
% Maak een tweede-ordesysteem aan zoals in oef 1 met wn = 10 en E = 0.02.
% Je mag daarvoor commando tf gebruiken. Behalve het maken van een Bode
% plot omdat je daar nog interessante dingens kan mee doen.
% > Maak bv. de FRF expliciet uitgerekend op de zelfgekozen frequenties
% met behulp van "freqresp". Kies w tussen 0 en 100 Hz. met een
% frequentieresolutie van 10^-3 Hz. Plot deze op het bode diagram.
% > Maak nu een zelfgekozen signaal aan dat je laat inwerken op het
% systeem met behulp van "lsim". Kies bv. een sinusoïdaal signaal en
% varieer de frequentie. Kun je zo het bodediagram reconstrueren?
% > Kun je zo ook de impulsreponsiefunctie h(t) bepalen?
% > Gebruik dan deze impulsresponsiefunctie om H(w) uit te rekenen, en plot
% dit over de bodeplot.
clear all;
wn = 10;
E = 0.02;
% Eerste deel:
FreqRes = 0.01;
N = 10000; % sampling gekozen
Ts = 1/N/FreqRes; % opgelet! niet 1/Fs, maybe omdat Fs = N/FreqRes???
tvec = [0:N-1]*Ts;
fvec = [0:N-1]/N/Ts;
% formule voor polen te verkijgen invullen
polen1 = -E * wn + 1j*wn*sqrt(1-E^2);
polen2 = -E * wn - 1j*wn*sqrt(1-E^2);
w = 2*pi*fvec; % vrijwillige keuze tussen 0 en 100. Maar w = 2*pi*fvec!
s = 1j * w; % omdat s = j*w (zie laplace)
H = wn^2 ./ (s-polen1) ./ (s-polen2);
subplot(2,2,1);
semilogx(w,db(H));
% tweede deel
polen = [polen1, polen2];
subplot(2,2,2);
plot(real(polen),imag(polen),'xr');
xlim([-15 15]);
ylim([-15 15]);
% vierde deel
subplot(2,2,3);
semilogx(w,angle(H));
% derde deel
h = 2*real(ifft(H,2*N));
h = h(1:N);
subplot(2,2,4);
plot(tvec, h);
% Laad de matlabfile "trueFRF" in, en maak een bodeplot. Dit is de FRF van
% een metalen plaat, getest in labo omstandigheden.
clf;
clear all;
load('trueFRF.mat');
subplot(2,2,1)
semilogx(fre,db(Hm(:,1)))
grid on;
subplot(2,2,2)
semilogx(fre,angle(Hm(:,1)))
grid on;
subplot(2,2,3)
semilogx(fre,db(Hm(:,2)),'r')
grid on;
subplot(2,2,4)
semilogx(fre,angle(Hm(:,2)),'r')
grid on;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment