[SIS] Vezba 5 z2

% KURS SISTEMI I SIGNALI
%
% VEZBA 5: Furijeova transformacija
%
% OSNOVNE OSOBINE FURIJEOVE TRANSFORMACIJE
%
% Osnovni transformacioni parovi
%
clear all
clc

% % % % disp ('Osnovni transformacioni parovi:');
syms w
syms t real
syms a T positive
% % % % f=dirac(t)
% % % % F=fourier(f)
% % % % ifourier(F,t)
% % % % disp ('-------------------------------------');
% % % % f=heaviside(t)
% % % % F=fourier(f)
% % % % ifourier(F,t)
% % % % disp ('-------------------------------------');
% % % % f=exp(-a*t)*heaviside(t)
% % % % F=fourier(f)
% % % % ifourier(F,t)
% % % % disp ('-------------------------------------');
% % % % f=exp(j*a*t)
% % % % F=fourier(f)
% % % % ifourier(F,t)
% % % % disp ('-------------------------------------');
% % % % f=cos(a*t)
% % % % F=fourier(f)
% % % % ifourier(F,t)
% % % % disp ('-------------------------------------');
% % % % f=sin(a*t)
% % % % F=fourier(f)
% % % % ifourier(F,t)
% % % % disp ('-------------------------------------');
%
% Linearnost
%
% % % % disp (sprintf('\nLinearnost:'));
% % % % f1=heaviside(-t);
% % % % f2=sin(t);
% % % % f3=dirac(t);
% % % % F1=expand(fourier(2*f1-3*f2+0.5*f3))
% % % % F2=expand(2*fourier(f1)-3*fourier(f2)+0.5*fourier(f3))
% % % % linearnost=(F1==F2)
%
% pomeranje u vremenskom domenu
%
% % % % disp (sprintf('\nPomeranje u vremenskom domenu:'));
% % % % syms T positive
% % % % f1=heaviside(t)
% % % % f2=heaviside(t-T)
% % % % F1=fourier(f1)
% % % % F2=fourier(f2)
% % % % pomeranje_u_vremenu=(exp(-j*T*w)*F1==F2)
%
% pomeranje u frekvencijskom domenu
%
% % % % disp (sprintf('\nPomeranje u frekvencijskom domenu:'));
% % % % syms W positive
% % % % F1=pi*dirac(w)-j/w
% % % % F2=pi*dirac(w-W)-j/(w-W)
% % % % f1=ifourier(F1, t)
% % % % f2=ifourier(F2, t)
% % % % pomeranje_u_frekvencijskom_domenu=(exp(j*W*t)*f1==f2)
%
% skaliranje u vremenskom domenu
%
%                         disp (sprintf('\nSkaliranje u vremenskom domenu:'));
%                         f1=dirac(t)
%                         f2=dirac(a*t)
%                         F1=fourier(f1)
%                         F2=fourier(f2)
%                         skaliranje_u_vremenu=(1/a*subs(F1,w,w/a)==F2)
%
% Furijeova transformacija izvoda funkcije
%
%                         disp (sprintf('\nFurijeova transformacija izvoda funkcije:'));
%                         f1=cos(t)
%                         F1=fourier(f1)
%                         f2=diff(f1,t)
%                         F2=fourier(f2)
%                         simplify((j*w)*F1)
%                         transformacija_izvoda=(F2==simplify((j*w)*F1))
%
% mnozenje sa t^n
%
% % % % disp (sprintf('\nMnozenje sa t^n:'));
% % % % f1=heaviside(t)
% % % % f2 = t*f1
% % % % F1=fourier(f1)
% % % % F2=fourier(f2)
% % % % mnozenje_sa_t=(simplify(j*diff(F1,w))==simplify(F2))
%
% Ogledanje u vremenu
%
% % % % disp (sprintf('\nOgledanje u vremenu:'));
% % % % f1=heaviside(t)
% % % % f2=heaviside(-t)
% % % % F1=fourier(f1)
% % % % F2=fourier(f2)
% % % % ogledanje_u_vremenu=(F2==subs(F1,w,-w))
%
% PRIMENA FURIJEOVE TRANSFORMACIJE U ANALIZI SISTEMA
%


disp (sprintf ('\nPrimena Furijeove transformacije u analizi sistema:'));
% Nalazenje odziva sistema resavanjem RUI u vremenskom domenu
u = sym ('cos(2*3.14*2*t)+0.5*cos(2*3.14*3*t)+0.5*cos(2*3.14*200*t)')

% IZVODI POBUDE
u_izvod_I = diff(u)
u_izvod_II = diff(u_izvod_I)

dif_jna = sprintf('D2y + 99.56681*Dy + 9070.20982*y = 0.89125*(%s)', char (u_izvod_II))
resenje1 = simplify(dsolve(dif_jna,'y(0)=0, Dy(0)=0','t'))
% resenje1 = simplify(dsolve(dif_jna,'y(1)=2, Dy(1)=0','t'))


disp('=====================================================');

% Nalazenje odziva koriscenjem Furijeove transformacije
syms w t G_w % ovo je oznaka za G(w), sliku od g(t)
slika_prvog_izvoda_g = simplify(j*w*G_w)
slika_drugog_izvoda_g = simplify((j*w)^2*G_w)

% TRAZI Grinovu FUNKCIJU.
% ODZIV NA Dirakovu POBUDU...
% Pri odredjivanju Grinove funkcije, svako y(t) na lijevoj strani RUI
% Mijenjam sa g(t), dok svako u(t) na desnoj mijenjam sa delta(t)
% Ovde imam drugi izvod pobude desno, tako da ovde trebam drugi izvod DELTE
slika_pobude = fourier(dirac(2, t))     % DRUGI IZOVD Diraka

RUI_leva_strana  = collect(slika_drugog_izvoda_g + 99.56681*slika_prvog_izvoda_g + 9070.20982*G_w,'G_w')
RUI_desna_strana = 0.89125*slika_pobude

% Odredjivanje prenosne karakteristike
prenosna_karakteristika = solve (sprintf ('%s = %s', char(RUI_leva_strana), char(RUI_desna_strana)), 'G_w')

% Amplitudska karakteristika izrazena u decibelima
% % Y = abs(X) returns the absolute value of each element in array X.
% % If X is complex, abs(X) returns the complex magnitude.
amplitudska_karakteristika = abs (prenosna_karakteristika);     % VRACA MODUO KOMP. BROJA
amplitudska_karakteristika_db = 20*log10(amplitudska_karakteristika);

% Fazna karakteristika izrazena u stepenima
%
fazna_karakteristika = 180/pi*atan (imag(prenosna_karakteristika)/real(prenosna_karakteristika));

% Crtanje amplitudske i fazne karakteristike
subplot (2, 1, 1), ezplot(amplitudska_karakteristika_db, [0, 1000]);
axis([0 1000 -90 10])
set(gca,'XScale','log');            % POSTAVLJA LOGARITAMSKU SKALU

title ('Amplitudska karakteristika sistema');
xlabel('\omega [rad/s]')
ylabel('20*log_1_0(|G(\omega)|) [dB]')  % Neki MATLAB-ov REGEX

subplot (2, 1, 2), ezplot(fazna_karakteristika, [0, 1000]);
axis([0 1000 -90 90])
set(gca,'XScale','log');
title ('Fazna karakteristika sistema');
xlabel('\omega [rad/s]')
ylabel('arg(|G(\omega)|) [stepeni]')


disp('=====================================================');
disp('=====================================================');
disp('=====================================================');

%%% --------------------------------------------------------------------
%%% --------------------------------------------------------------------
%%% --------------------------------------------------------------------
%%% RESENJE RACUNANJEM AMPLITUDSKE I FAZNE KARAKTERISTIKE
%%% PODRAZUMEVA NULTE POCETNE USLOVE, ZATO SE POCETAK ODZIVA U t-domenu
%%% I ODZIVA PREKO AMPL. I FAZNE KKE RAZLIKUJE... Takodje se razlikuju
%%% KADA NE KAZEM USLOVE U Nuli, nego kazem npr y(1) = nesto
%%% --------------------------------------------------------------------
%%% --------------------------------------------------------------------
%%% --------------------------------------------------------------------


% Racunanje odziva sistema na ulazni signal koriscenjem amplitudske i fazne
% karakteristike
% Vrednost amplitudske i fazne karakteristike na ucestanosti PRVE komponente iz
% ulaznog signala

% % %     u = sym ('cos(2*3.14*2*t)+0.5*cos(2*3.14*3*t)+0.5*cos(2*3.14*200*t)')
% % %     u = sym ('cos(2*3.14*2*t)+0.5*cos(2*3.14*3*t)+0.5*cos(2*3.14*200*t)')
% % %     u = sym ('cos(2*3.14*2*t)+0.5*cos(2*3.14*3*t)+0.5*cos(2*3.14*200*t)')

w = 2*3.14*2;
amplituda_1 = eval (amplitudska_karakteristika)
faza_1      = eval (fazna_karakteristika)

% Vrednost amplitudske i fazne karakteristike na ucestanosti DRUGE komponente iz
% ulaznog signala
w = 2*3.14*3;
amplituda_2 = eval (amplitudska_karakteristika)
faza_2      = eval (fazna_karakteristika)

% Vrednost amplitudske i fazne karakteristike na ucestanosti TRECE komponente iz
% ulaznog signala
w = 2*3.14*200;
amplituda_3 = eval (amplitudska_karakteristika)
faza_3      = eval (fazna_karakteristika)


% Izlazni signal sistema racunat na osnovu amplitudske i fazne
% karakteristike
resenje2 = amplituda_1*cos(2*3.14*2*t + faza_1) + ...
           amplituda_2*0.5*cos(2*3.14*3*t + faza_2) + ...
           amplituda_3*0.5*cos(2*3.14*200*t + faza_3)

figure
% Crtanje ulaznog signala
u = sym ('cos(2*3.14*2*t)+0.5*cos(2*3.14*3*t)+0.5*cos(2*3.14*200*t)');
subplot(3,1,1)
ezplot(u,[0, 2.4]);
axis([0 2.4 -2 2])

xlabel('t (sec)')
ylabel('u(t)')
grid on
title('Ulazni signal')

% Crtanje resenja dobijenog resavanjem RUI u vremenskom domenu
subplot(3,1,2)
ezplot(resenje1, [0, 2.4]);
axis([0 2.4 -2 2])
xlabel('t (sec)')
ylabel('y(t)')
grid on
title('Resenje dobijeno u t-domenu')

% Crtanje resenja dobijenog na osnovu amplitudske i fazne karakteristike
subplot(3,1,3)
ezplot(resenje2, [0, 2.4]);
axis([0 2.4 -2 2])
xlabel('t (sec)')
ylabel('y(t)')
grid on
title('Resenje dobijeno koriscenjem amplitudske i fazne karakteristike')


simplify(resenje1)
simplify(resenje2)