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

Teekkari1 = 'Teemu Teekkari 12345A'; %LUE HUOM! ALLA!
Teekkari2 = 'Liisa Teekkari 54321B';

%HUOM! TÄHÄN TYÖHÖN LIITTYVÄÄ TIETOA ON EEFORUMILLA!

%MAKRO: TYÖ 3, osa C1
%IDEAALINEN KAISTANPÄÄSTÖSUODATIN
%MAKRO: ideal2005.m

%Tässä työssä mallinnetaan ideaalista kaistanpäästösuodatinta, jonka läpi ajetaan ääninäytteitä.
%Suodattimen kaistaleveyttä B ja keskitaajuutta fc voidaan muuttaa halutuiksi.

%Perusasetukset:
Fs = 11025; %Näytenopeus (Hz tai S/s eli samples/s) - samalla äänikortille sopiva arvo.
Ts = 1/Fs; %Näyteväli.

T = 2; %Näytteistysaika 2 s - samalla ääninäytteen kesto.
t = -0.5*T:Ts:0.5*T-Ts; %Aika -0,5*T sekunnista lähes 0,5*T sekuntiin.
N = T/Ts; %Näytemäärä eli näytteistysaika T jaettuna näytevälillä Ts.
f0 = Fs/N; %Spektrin taajuuksien askellusväli hertseinä.

f = -Fs/2:f0:Fs/2-f0; %Taajuusvektori spektriä varten.

%Aluksi määritellään kaistanpäästösuodattimen ala- ja ylärajataajuudet.
%Alla olevia arvoja voit myöhemmin muuttaa sopivimmiksi!

f1 = 1000; %Vastaa "spektrilaatikon" vasenta reunaa.
f2 = 2000; %Vastaa "spektrilaatikon" oikeaa reunaa.

%Näin ollen on suodattimen kaistaleveys ns. kantataajuisena 1-puoleisessa taajuusalueessa
B = (f2-f1)/2;

%Valitaan suodattimen vahvistukseksi suunnilleen 1, joten sinc-pulssin korkeus on A = 0,5 eli
A = 0.5; 

%Nyt voidaan lausua impulssivaste ideaaliselle kantataajuiselle
%alipäästösuodattimelle:
h1t = A*(sin(pi*2*B*t)./(pi*2*B*t)); %Eli h1(t) = A*sinc(2Bt).

%Lasketaan suodattimen keskitaajuus fc
fc = (f1+f2)/2;

%"Siirretään laatikko" taajuusalueessa fo:n ympärille,
%jolloin aika-alueen esitys on taajuussiirretylle suodattimen impulssivasteelle

ht = h1t.*cos(2*pi*fc*t);

%Kyseessä on siis lineaarinen modulaatio, jota hyödynnetään mm. modeemien yhteydessä.

%Esitetään Hf 2-puoleisessa taajuusalueessa:
Hf = abs(fftshift(fft(ht))/(0.5*N));
%Eli ensin tehdään DFT/FFT, joka skaalataan oikeaksi jakamalla näytemäärällä N.
%Sitten siirretään 2-puoleinen spektri oikeaan paikkaan komennolla fftshift.
%Lopuksi esitetään ampltudin itseisarvo taajuusalueessa komennolla abs.

%Ääninäytteiden haku: 
[xt] = wavread('c:\Temp\Sinkut\x1t.wav');

%Tehdään xt:lle Fourier-muunnos DFT/FFT:llä, ja käytetään taas 2-puoleista spektriä
Xf=abs(fftshift(fft(xt))/N);

%Konvoluutio eli konvoloidaan omaa ääntäsi ideaalisen kaistanpäästösuodattimen
%impulssivasteen kanssa eli suodatetaan äänestäsi kuuluviin vaikkapa taajuusalue
%1000 - 2000 Hz:

yt = conv(ht,xt);

%Koneella on melkoinen työ eli käännetään ht aikahetken 0 ympäri,
%vedetään h(lambda) N/2 näytettä vasemmalle ja siirretään näytteittäin
%h(lambda)- ja x(lambda)-arvot ns. päällekkäin,
%kerrotaan sitten päällekkäiset arvot keskenään, ja lopuksi
%lasketaan kertolaskujen arvot yhteen sekä sijoitetaan "integrointitulos" kohtaan t.

%yt:n pituus vektorialkioina on ht:n pituus + xt:n pituus -1.
%Miksikö? Ajattele, että ht:n ja xt:n pituudet ovat 2 ja 2.
%Tulee vain 3 tilannetta, jolloin oikeita päällekäisyyksiä on.
%Huomiodaan tämä jatkossa, koska virheilmoitus vektoreiden pituuksista tulee herkästi.

%Lasketaan lähtöspektri Yf
Yf=abs(fftshift(fft(yt))/N);

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

%Lopuksi tulostetaan kuvat, ja kuunnellaan suodatettua ääntä.
%6 kuvaa, 2 rinnakkain, vasemmalla aika-alue ja oikealla taajuusalue:

subplot(3,2,1); plot(t,xt); axis([-0.5*T 0.5*T 1.2*min(xt) 1.2*max(xt)]);
%Asteikot skaalattu automaattisiksi, jottei ns. näpräämiseen menisi aikaa.
xlabel('Aika/s'); ylabel('Amplitudi/V'); grid on; zoom on;

title({['Kurssi: S-72.060 K-2005 Laboratoriotyö 3, osa C1'];...
      ['1000 - 2000 Hz ideaalinen kaistanpäästösuodatin, Minun ääni'];...    
      ['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(3,2,2); plot(f,Xf); axis([-Fs/2 Fs/2 0 1.2*max(Xf)]);
xlabel('Taajuus/Hz'); ylabel('Amplitudi/V'); grid on; zoom on;

subplot(3,2,3); plot(t,ht); axis([-0.5*T 0.5*T 1.2*min(ht) 1.2*max(ht)]);
xlabel('Aika/s'); ylabel('Amplitudi'); grid on; zoom on;

subplot(3,2,4); plot(f,Hf); axis([-Fs/2 Fs/2 0 1.2*max(Hf)]);
xlabel('Taajuus/Hz'); ylabel('Amplitudi'); 
grid on; zoom on;

subplot(3,2,5); t1 = -0.5*T+Ts:Ts:1.5*T-Ts; %Lasketaan konvoluutiotulokselle aikavektori oikean mittaiseksi.
plot(t1,yt); axis([-0.5*T T 1.2*min(yt) 1.2*max(yt)]); xlabel('Aika/s'); ylabel('Amplitudi/V');
grid on; zoom on;

subplot(3,2,6); f1 = -Fs/2+f0/2:f0/2:Fs/2-f0/2; %Lasketaan Yf:n taajuusasteikko oikeaksi.
plot(f1,Yf); axis([-Fs/2 Fs/2 0 1.2*max(Yf)]); xlabel('Taajuus/Hz'); ylabel('Amplitudi/V');
grid on; zoom on;

%Kuunnellaan, miltä suodatettu ääni kuullostaa:
soundsc(yt,Fs);

shg;

