% MATLAB-TIPS TIL OPPGAVE 5.41,
% BL.A. FUNKSJONENE POLYVAL og FILTER,
% og generell grunnleggende matlab,
% i "Digital Signal Processing" av Ashok Ambardar
%
% (laget av Henrik G. Sundt, hsundt@notam02.no, 9.oktober 2009)

% Se ellers kapittel 11 i boken for diverse nyttige Matlab-funksjoner.



% Elementer og funksjoner fra dette dokumentet kan brukes i løsningen av
% oppgaven, med innsatt andre verdier.


% Det er lagt inn "pause;" etter hvert plott, slik at hele
% avsnitt eller hele filen kan kjøres på en gang.
% Man må derfor trykke enter etter hvert plott.

% Hvis man tar bort ";" på slutten av linjen, vil Matlab skrive ut resultatet.
% Det kan være nyttig for å sjekke f.eks. at man har fått riktig type
% resultat (vektor eller enkeltverdi).




% EKSEMPEL 1
% -- GENERERING AV ET SIGNAL, OG EVALUERING AV FREKVENSINNHOLDET


% Generering av et signal x:

n=0:300;           % lager en vektor med verdier fra 0-300
x=0.5 .^n;         % .^ betyr at 0.5^n skal utføres på hvert element i vektoren n

% Plotting av x.
plot(n,x); pause;

% Lager et signal xn = x med litt støy:

xn=x+rand(size(x))/10;
plot(n,xn); pause;

% Plotter begge signaler i samme vindu:
subplot(2,1,1), plot(n,x);
subplot(2,1,2), plot(n,xn); pause;

% eller i samme plott med ulik farge:
subplot(1,1,1), plot(n,[x;xn]); pause;



% Vi skal regne i Matlab med Z-transformen til signalet x, som er:
% X(z) = z/(z-0.5);

% Oppsett av variabler for evaluering av verdier i Z-planet
% langs enhetssirkelen fra 0 til pi:

F=0:0.005:0.5;  % en vektor med steg 0.005
w=2*pi*F;
Z=exp(-j*w);    % Vektoren Z får verdier som går langs enhetssirkelen fra 0 til -pi

% plotting av absoluttverdi og fase til Z:
subplot(2,1,1), plot(F,abs(Z));   % absoluttverdien Z (alltid 1)
subplot(2,1,2), plot(F,angle(Z)); % fasen til Z
pause;

% Beregning av X, som blir evaluert langs enhetssirkelen i
% verdiene Z, ved direkte uttrykk av X(Z):

X = Z./(Z-0.5);  % OBS: Dette er en annen og enklere formel enn den som ble
                 % gjennomgått på regneverkstedet. Gir også et annet plott.

% Plotting av frekvensinnhold til signalet:
% Dette er det samme som å plotte frekvensresponsen til et system,
% hvis vi har impulsresponsen til systemet (= signalet som kommer ut,
% hvis man sender en impuls inn).
clf;
plot(F,abs(X));
pause;
 

% Alternativ metode:
% Ved hjelp av Matlab-funksjonen polyval(V,Z) kan vi evaluere Z ved
% å lage en vektor som representerer et polynom. Se "help polyval"
% i Matlab.

% Man finner først en vektor for polynomet:
% 4.-gradspolynomet A*z^4 + B*z^3 + C*z^2 + D*z + E blir:
%
% V = [A B C D E]
%
% 1.grads-polynomet D*z + E blir bare:
%
% V = [D E]
%
% Når man bruker funksjonen polyval(V,Z) blir dette polynomet beregnet for
% hvert element i Z, og resultater er en vektor.

% Vi lager tilsvarende funksjon som X over, men vha polynomer:
% Først gjør vi litt om på skrivemåten av X, for å tilrettelegge
% for vektorform:
% X = N1/D1 = (1*Z+0)./(1*Z-0.5)
% og lager følgende polynomer

N1=[1 0];     % Z+0
D1=[1 -0.5];  % Z-0.5

% evaluer med Matlab-funksjonen polyval, med Z som over:

XP=polyval(N1,Z)./polyval(D1,Z);

% plotter X og XP i samme vindu og vi ser at de er like:
subplot(2,1,1), plot(F,abs(X));
subplot(2,1,2), plot(F,abs(XP)); pause;

% I dette tilfellet er det hipp som happ om man bruker et direkte
% uttrykk for X, eller om man bruker polyval, men
% polyval er spesielt nyttig når man har polynomiske uttrykk som varierer
% under kjøring, som i eksempel 2 under.






% EKSEMPEL 2
% -- ET N-PUNKTS MOVING AVERAGE FILTER


% Et moving average filter h av lengde N.
N = 9;
h=ones(1,N)/N;

% Filtrering av signal med støy, med moving-average filteret.
% h er her et FIR-filter:
y=filter(h,1,xn);  % xn fra eksempel 1


% Plotting originalt signal, signal med støy, og filtrert signal.
clf;
plot(n,[x; xn; y]);   % x og xn fra eksempel 1
pause;

% Vi beregner frekvensresponsen Y til filteret:
Y=polyval(h,Z);

% For N lik 9 kan du evt prøve direkte uttrykk, som gir samme svar:
% Y = (Z.^8 + Z.^7 + Z.^6 + Z.^5 + Z.^4 + Z.^3 + Z.^2 + Z + 1)/9;
% Det er her fordelen ved polyval blir tydelig, f.eks. ved
% endring av N.

% Plotter amplitude og fase til filteret:
subplot(2,1,1), plot(F,abs(Y));
subplot(2,1,2), plot(F,angle(Y)); pause;


% Det andre parameteret i filter, som her bare er 1, er den tilbakekoblede delen (under brøkstreken i H).
%
% Man kan altså filtrere med f.eks. H=z/(z+0.5) ved å bruke
% y=filter([1 0],[1 0.5],x); 
