%% B9924 - LTI-suotimen analyysi
%
% $$H(z) = \frac{1 + 0.6 z^{-1}}{1 - 0.9 z^{-1}}, \ |z| > 0.9$$
%

% Jukka Parviainen, parvi@tkk.fi, 21.2.2011

%% Suotimen kertoimet Matlabia (Octavea) varten
% H(z) = [1 + 0.6 z^(-1)] / [1 - 0.9 z^(-1)] = B(z)/A(z)
%
% Osoittajapolynomi B(z): eteenpäin laskenta, nollat, "FIR"
% B(z) = [1   +0.6 z^(-1)]    
B      = [1    0.6];
% Nimittäjäpolynomi A(z): takaisinkytketty silmukka, navat, "IIR"
% A(z) = [1   -0.9 z^(-1)]
A      = [1   -0.9];
tf2latex(B, A)   % kurssin oma funktio, kts. R05-kierros Nopassa

%% (a) Lasketaan navat ja nollat, ja hahmotellaan napanollakuvio
% Nollat ovat B(z):n juuret eli B(z)==0.
disp('Nollat:');
disp(roots(B));
% Navat ovat A(z):n juuret eli A(z)==0.
disp('Navat:');
disp(roots(A));
% Tai suoraan Matlabin funktiolla
[z,p,k] = tf2zp(B, A);
% Napanollakuvio: kompleksitaso, yksikköympyrä, 
% navat rasteja "x", nollat ympyröitä "o"
figure(11); clf;
zplane(B, A);
%% (b) Kirjoita taajuusvaste H(e^jw)
% Korvaus 
%
% $$z \leftrightarrow e^{j \omega}$$
%
% jolloin saadaan
%
% $$H(e^{j\omega}) = \frac{1 + 0.6 e^{-j \omega}}{1 - 0.9 e^{-j \omega}}$$
%% (c) Magnitudivaste
figure(2); clf;
[H, w] = freqz(B, A);
plot(w/pi, abs(H)); 
grid on;
xlabel('normalisoitu kulmataajuus \omega (\times \pi rad/näyte)');
title('Magnitudivaste');
% sama vielä dB-asteikolla
figure(3); clf;
plot(w/pi, 20*log10(abs(H)));
grid on;
xlabel('normalisoitu kulmataajuus \omega (\times \pi rad/näyte)');
ylabel('dB');
title('Magnitudivaste');
%% (c) Magnitudivaste 3D-esityksenä:
% z-muunnos ja tulkinta yksikköympyrällä H(e^jw).
% Jos käytät Matlabia, niin pääset pyörittämään 3D-visualisointia.
figure(4); clf;
[X, Y] = meshgrid(-1.45: 0.03 : 1.45);
C = X + j*Y;
Z  = abs((C + 0.6)./(C - 0.9));  % == [z+0.6]/[z-0.9]
mesh(X, Y, Z);
axis([-1.3 1.3 -1.3 1.3 -3 13]);  % zoomaus
xlabel('Re');
ylabel('Im');
hold on;
%uu = contour(X, Y, Z);
%set(uu, 'LineWidth', 3);
%contour3(X, Y, Z);
wd = 2*pi*[0 : 0.01 : 1];
rx = real(exp(j*wd));
ry = imag(exp(j*wd));
% H(e^jw) 
rz = (exp(j*wd) + 0.6)./(exp(j*wd) - 0.9); % == [z+0.6]/[z-0.9]
pp = plot3(rx, ry, rz);
set(pp,'Linewidth', 3, 'Color', [1 0 0]);
% yksikköympyrä
pq = plot3(rx, ry, zeros(size(rx)));
set(pq,'Linewidth', 2, 'LineStyle', '-.');
% nollat
pze = plot3(real(roots(B)), imag(roots(B)), 0, 'o');
set(pze,'MarkerSize', 12, 'MarkerEdgeColor', [1 0 0], 'LineWidth', 2);
% navat
ppo = plot3(real(roots(A)), imag(roots(A)), 0, 'x');
set(ppo,'MarkerSize', 12, 'MarkerEdgeColor', [1 0 0], 'LineWidth', 2);
% koordinaattiakselit
paxx = plot3([-1.3 1.3],[0 0],[0 0]);
set(paxx,'LineWidth',2,'Color',[0 1 0]);
paxy = plot3([0 0],[-1.3 1.3],[0 0]);
set(paxy,'LineWidth',2,'Color',[0 1 0]);
hidden off
rotate3d
%% (d) Vaihevaste <H(e^jw)
figure(5); clf;
plot(w/pi, angle(H)); % phase response
grid on; 
xlabel('normalisoitu kulmataajuus \omega (\times \pi rad/näyte)');
ylabel('rad');
title('Vaihevaste');
%% (e) Differenssiyhtälö
[tftxt, tftex, frtxt, frtex, dtxt, dtex] = tf2latex(B, A);
['Differenssiyhtälö: ' dtxt]
%% (f) Lohkokaavioesitys
%% (g) Impulssivaste
figure(6); clf;
[h, n] = impz(B, A); % impulse response
stem(n, h);
grid on; 
xlabel('aikaindeksi n');
title('Impulssivaste h[n]');
imp2latex(h, n)   % kurssin oma funktio, kts. R05-kierros Nopassa
%% (h) Onko IIR vai FIR?
% IIR: takaisinkytkentä, äärettömän pitkä h[n], H(z):n alakerrassa
%      polynomi, y[n]:n laskemiseksi y[n]:n edellisiä arvoja
% FIR: ei yhtäkään takaisinkytkentää, äärellisen pitkä h[n],
%      H(z) pelkkä yksi polynomi, aina stabiili, mahdollisuus
%      lineaarivaiheiseen suotimeen.
%% (i) Asteluku?
max([length(B)-1; length(A)-1])
max([length(roots(B)); length(roots(A))])
%% (j) Kausaalinen?
% Kausaalinen järjestelmä ei ennusta mitään eli se ei tarvitse
% tulevaisuuden arvoja y[n]:n laskemiseksi.
%% (k) Stabiili?
% Stabiili järjestelmä pitää ulostulon rajattuna jos sisääntulo
% on rajattu. Vrt. mikrofoni liian lähellä kaiutinta.