%TEKNILLINEN KORKEAKOULU
%Tietoliikennelaboratorio
%S-72.060 Signaalit ja järjestelmät
%Seppo Saastamoinen 27.1.2005 ver. 1.2
clear; close all; %Nollataan muistit ja kuvat

Teekkari1 = 'Teemu Teekkari 12345A';
Teekkari2 = 'Liisa Teekkari 54321B';

%MAKRO: TYÖ 1
%SIGNAALIANALYYSI
%MAKRO: signanal.m

%Tämä makro laskee ja näyttää signaalin aika- ja taajuusalueessa
%seuravasti (4 kuvaa/tulostussivu).

%1. sivu - 1-puoleinen taajuusalue:
%a) aika-alueessa 2 kuvaa, jotka ovat samoja.
%Kaksi ylintä kuvaa zoomataan, jotta signaalin yksityiskohtia
%on helpompi tarkastella aika-alueessa.
%b) taajuusalueessa (kaksi alinta kuvaa) nähdään signaalin
%amplitudi (yksikkö V/Hz tai dBV/Hz) ja vaihekulma asteina 

%2. sivu - 2-puoleinen taajuusalue:
%a) ja b) kuten edellä, nyt käytössä on 2-puoleinen
%spektri 1-puoleisen sijaan.


%------------- ALKUASETUKSET ---------------------------------

%Matlabissa signaalit niin aika- kuin taajuusalueessa ovat
%näytteinä (samples). Yleensä lähtökohtana on aika-alue,
%jolloin näytetaajuus (sample frequency, Fs tai fs) on tunnettava.
%Fs voidaan antaa itse tai arvo saadaan signaalin mukana.

%Aluksi asetetaan Fs-arvo itse, koska signaali tehdään työohjeen
%mukaan. Työohjeesta on apua.
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Fs = 44100;
%Lasketaan näyteväli Ts (time between samples)
Ts = 1/Fs;
%Määritellään signaalin kestoaika t (aikavektori).
%Aika alkaa nollasta, jossa on 1. näyte. Askellusväli
%vektorissa on Ts. Loppuaika on noin 1 s (1-Ts).
t = [0:Ts:1-Ts]; 

%------------- SIGNAALI AIKA-ALUEESSA -------------------------

%Aluksi tehdään itse signaaleja. Alla on muutamia käskyjä,
%joiden avulla esimerkiksi suorakaideaallon tekeminen
%onnistuu.

%Huom! Kommentoi %-merkillä ne käskyt, joita ei käytetä.

%Asetetaan aluksi jaksollisen signaalin taajuus f.
%Työohjeesta löytyy sopiva arvo.
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
f = 1000;

%a) Kosiniaalto (jännite välillä +1 .. -1 V):
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
xt = 1*cos(2*pi*f*t);

%b) Suorakaideaalto (jännite välillä +1 .. -1 V).
%Pulssisuhde eli D (duty cycle) on ylhäällä olevan
%pulssin keston suhde jaksonaikaan (arvo %:eina).
%Työohjeessa on mainittu D:n arvo!
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%xt = 1*square(2*pi*f*t, D);

%c) Kolmioaalto (jännite välillä +1 .. -1 V). Nyt arvo D
%on aallon nousevan reunan keston suhde jaksonaikaan
%absoluuttiarvona, esimerkiksi 0.5.
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%xt = 1*sawtooth(2*pi*f*t, D);

%d) Joku muu itse tehty signaali.
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%xt = .....

%Alla olevilla komennoilla voidaan signaali hakea tiedostosta.

%Äänisignaalin hakeminen. Muuta työohjeen/assistentin ohjeen
%mukaan hakemistopolku ja tiedostonimi oikeiksi!
%Komennon avulla vektorin xt sisällöksi tulee signaali ja
%näytenopeuden Fs arvo saadaan samalla.

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%[xt,Fs] = wavread('c:\Temp\Sinkut\nimi1.wav');
%N = max(size(xt)); Ts = 1/Fs; t = 0:Ts:(N-1)*Ts;

%Kuunnellaan signaalia! Komento soundsc skaalaa amplitudit niin,
%että ääni ei säröydy eikä kuulu liian hiljaa, ts. pakottaa
%amplitudit välille -1 ... +1).
%soundsc(xt,Fs);

%------------- SIGNAALI TAAJUUSALUEESSA -----------------------

%Lasketaan näytemäärä aika-alueessa
N = max(size(xt)); 

%Lasketaan (kompleksinen) Fourier-muunnos signaalille xt
%ja skaalataan arvot oikeiksi jakamalla näytemäärällä N.
Xf = fft(xt)/N;

%Lasketaan taajuusmuunnoksen itseisarvo ja kerrotaan 2:lla 1-puoleista
%taajuusaluetta varten.
Xf1puol = 2*abs(Xf); 

%Lasketaan taajuusmuunnoksen itseisarvo ja muutetaan se 2-puoleista
%taajuusaluetta varten komennolla fftshift.
Xf2puol = fftshift(abs(Xf));

%Valmistellaan vaihespektrin laskemista. Jotta kulma-arvot olisivat
%oikein, lisätään spektrin reaaliosaan virhe 0.0001. Matlab tekee
%hyvin pieniä virheitä spektrin laskennassa. Esimerkiksi
%taajuusalkion arvo 1e-18 + 1e-18i antaisi vaihekulmaksi 45 astetta,
%vaikka oikea arvo on 0 astetta, koska amplitudi on nolla.
Xf1puolp = Xf + 0.0001;
Xf2puolp = fftshift(Xf) + 0.0001;

%Lasketaan 2- ja 1-puoleisten spektrien vaihekulmat asteina
%angle-komennolla.
Vaihe1puol = (180/pi)*angle(Xf1puolp);
Vaihe2puol = (180/pi)*angle(Xf2puolp);

%Taajuusalue on myös diskreetti eli spektriarvoja lasketaan
%tietyn kaavan mukaan. Taajuuden askellusväli on
f0 = Fs/N;

%Tehdään tulostusta varten 1- ja 2-puoleiset taajuusvektorit.
f1puol = [0:f0:(N-1)*f0];
f2puol = [-N/2*f0:f0:N/2*f0-1];

%------------- SIGNAALEIDEN TULOSTUS ---------------------------------

%Tulostetaan 1. sivu näytölle figure-komennolla
figure;

%Tulostetaan 4 kuvaa allekain
subplot(4,1,1);
%Tulostetaan haettu signaali aika-alueessa
plot(t,xt);
%Tehdään x- ja y-akselit järjestyksessä
%xmin xmax ymin ymax
axis([-0.05 1.05*max(t) 1.2*min(xt) 1.2*max(xt)]); 
%Kuvan pystyasteikon nimi
ylabel('Signaali/V');
%Kuvan vaaka-akselin nimi
xlabel('Aika/s');
%Asetetaan asteikkoviivoitus päälle
grid on;
%Sallitaan kuvan zoomaus
zoom on; 

%Haetaan päiväys ja kellonaika tulostusta varten eli
h = fix(clock);

%Luodaan kuville yhteinen otsikkoteksti, joka voi olla
%monirivinen. Muuta tiedot työhön sopiviksi.
%Huom! Kolme pistettä ... tarkoittaa Matlabille, että
%suoritus jatkuu seuraavalla rivillä. Näin voit pitkänkin
%komennon tai tekstin jakaa riveille.
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
title({['Kurssi: S-72.060 K-2005 Laboratoriotyö 1, osa XYZ'];...
      ['Signaali: 1000 Hz kosiniaalto, 1-puoleinen taajuusalue'];...    
      ['Tekijät: ' Teekkari1 ' ja ' Teekkari2];...
      [int2str(h(3)) '.' int2str(h(2)) '.' int2str(h(1)) ' klo '...
          int2str(h(4)) ':' int2str(h(5)) ':' int2str(h(6))]});

subplot(4,1,2); plot(t,xt);
axis([-0.05 1.05*max(t) 1.2*min(xt) 1.2*max(xt)]);
ylabel('Signaali/V'); xlabel('Aika/s');
grid on; zoom on;

subplot(4,1,3);
%Tulostetaan signaalin amplitudispektri
%1-puoleisessa taajuusalueessa
plot(f1puol,Xf1puol);
axis([-100 Fs -0.001 1.2*max(Xf1puol)]);
ylabel('Amplitudi/V/Hz'); xlabel('Taajuus/Hz');
grid on; zoom on;

subplot(4,1,4);
%Tulostetaan signaalin vaihespektri
%1-puoleisessa taajuusalueessa
plot(f1puol,Vaihe1puol); 
axis([-100 Fs -200 200]);
ylabel('Vaihekulma/astetta'); xlabel('Taajuus/Hz');
grid on; zoom on;

shg;

%Tulostetaan 2. sivu näytölle figure-komennolla.
%Nyt käytössä on 2-puoleinen taajuusalue.
figure;

subplot(4,1,1); plot(t,xt);
axis([-0.05 1.05*max(t) 1.2*min(xt) 1.2*max(xt)]); 
ylabel('Signaali/V'); xlabel('Aika/s');
grid on; zoom on; 

%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
title({['Kurssi: S-72.060 K-2005 Laboratoriotyö 1, osa XYZ'];...
      ['Signaali: 1000 Hz kosiniaalto, 2-puoleinen taajuusalue'];...    
      ['Tekijät: ' Teekkari1 ' ja ' Teekkari2];...
      [int2str(h(3)) '.' int2str(h(2)) '.' int2str(h(1)) ' klo '...
          int2str(h(4)) ':' int2str(h(5)) ':' int2str(h(6))]});

subplot(4,1,2); plot(t,xt);
axis([-0.05 1.05*max(t) 1.2*min(xt) 1.2*max(xt)]);
ylabel('Signaali/V'); xlabel('Aika/s');
grid on; zoom on;

subplot(4,1,3);
%Tulostetaan signaalin amplitudispektri
%2-puoleisessa taajuusalueessa
plot(f2puol,Xf2puol);
axis([-100-Fs/2 100+Fs/2 1.2*min(Xf2puol) 1.2*max(Xf2puol)]);
ylabel('Amplitudi/V/Hz'); xlabel('Taajuus/Hz');
grid on; zoom on;

subplot(4,1,4);
%Tulostetaan signaalin vaihespektri
%2-puoleisessa taajuusalueessa
plot(f2puol,Vaihe2puol); 
axis([-100-Fs/2 100+Fs/2 -200 200]);
ylabel('Vaihekulma/astetta'); xlabel('Taajuus/Hz');
grid on; zoom on;
