function esim1
% ESIM1 Demo LTI-suodatuksesta
%
% Esitellään sirinäinen signaali ja suodatetaan sitä niin,
% että jäljelle jää lankapuhelimen kaista alle 3300 Hz.

% Jukka Parviainen, parvi@hut.fi, 12.9.2002
% Tiedostot saatavilla 
% http://www.cis.hut.fi/Opinnot/T-61.246/Demo/esim1.shtml
%
% Alkutoimenpiteet: tarkista että esityskoneessasi on 
% kaiuittimet tai salin kyseessä ollen audioliitäntä.
% Lataa tiedostot työhakemistoon tai aseta Matlabissa
% addpath-komennolla polku tunnetuksi.
% Kutista Matlab-työikkuna suhteellisen pieneksi ruudun alareunaan,
% sijoita piirtoikkuna, esimerkiksi ensimmäisen vaiheen jälkeen,
% sopivaan paikkaan, että se näkyy.
%
% Tiedostot: esim1.m, esim1_input.wav
%
% HUOM! Demoluonteen takia suurin osa _oikeasta_ laskennasta
% tehdään aivan koodin alussa.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;      % sulkee kaikki ikkunat
clc;            % puhdistaa ruudun

fprintf(1,'****************************************\n');
fprintf(1,'****************************************\n');
fprintf(1,'**************  esim1.m  ***************\n');
fprintf(1,'****************************************\n');
fprintf(1,'**** Suomen Virallinen Suodatusdemo  ***\n');
fprintf(1,'****************************************\n');
fprintf(1,'\n\n');

fprintf(1,'** Odota hetki...\n\n');

% Luetaan äänisignaali sisään Matlabiin
[y, fs, nbits] = wavread('esim1.wav'); 

Wp = 2*3000/fs; % Päästökaistan loppu
Ws = 2*3300/fs; % Estokaistan alku
Rp = 0.2;
Rs = 40;
[N,Wn] = ellipord(Wp, Ws, Rp, Rs);
[b,a] = ellip(N, Rp, Rs, Wn);
[hh,ww]=freqz(b,a,1024,fs);   % suotimen amplitudi- ja vaihevaste

Y = fft(y,2^14);        % Lasketaan spektriä 2^14 ensimm. pisteestä,
                        % jotta ei vie liikaa aikaa...

ysuod = filter(b,a,y);  % Tehdään suodatus valmiiksi

Ysuod = fft(ysuod,2^14); % Lasketaan vastaavasti spektriä vasteesta

fprintf(1,'** Paina ENTER...\n\n');
pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

fprintf(1,'Kuunnellaan ESIM1:n sirinäistä ääntä.\n\n');

t = linspace(0, length(y)/fs, length(y));

soundsc(y, fs);  %Kuunnellaan

fprintf(1,'** Paina ENTER äänen kuuntelun jälkeen...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

w = linspace(0, fs, length(Y));
subplot(311);
plot(w, abs(Y));
axis([0 fs/2 0 1.1*max(abs(Y))]);
grid on;
title(['Signaalin spektri 0..' num2str(fs/2) ' hertziä']);
xlabel('taajuus (Hz)');

fprintf(1,'Spektrissä on esitetty signaali taajuuskomponenttien avulla.\n\n');

fprintf(1,'** Paina ENTER...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

fprintf(1,'Äänessä on paljon korkeataajuista sirinää.\n\n');

fprintf(1,'Tehdään alipäästösuodin, joka päästää läpi\n vain lankapuhelintaajuudet.\n\n');
fprintf(1,'Aikatason nopeat muutokset signaalissa\n vastaavat korkeita taajuuksia.\n');
fprintf(1,'Aikatason hitaat muutokset ovat matalia taajuuksia.\n\n');

fprintf(1,'** Paina ENTER...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Seuraava pelkkää pluffia
h = waitbar(0,'Lasketaan suotimen kertoimia...');
for i=1:20,
  pause(0.02);
  waitbar(i/20,h)
end
close(h)

subplot(312);
ku = plot(ww,abs(hh).^2);
axis([0 ww(end) 0 1.05]);
set(ku,'LineWidth',2);
grid on;
title('Elliptinen alipäästösuodin H(z)');
subplot(311); xlabel('');
subplot(312); xlabel('taajuus (Hz)');

fprintf(1,'Toteutettiin elliptinen alipäästösuodin (kesk.).\n');
fprintf(1,'Tämä suodin käsittelee signaalia niin, \n');
fprintf(1,'että matalat taajuudet jäävät jäljelle\n\n');

fprintf(1,'** Paina ENTER...\n\n');
pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

fprintf(1,'Suodatetaan suhinainen signaali...\n\n');

fprintf(1,'\n Konvoluutioteoreema: aikatason konvoluutio on taajuustason tulo\n');
fprintf(1,'\n');
fprintf(1,'                 --------------\n');
fprintf(1,'                 |            |\n');
fprintf(1,'  X(e^jw)  ----> |  H(e^jw)   | ---->  Y(e^jw) = H(e^jw) . X(e^jw)\n');
fprintf(1,'                 |            |\n');
fprintf(1,'                 --------------\n');
fprintf(1,'\n');
fprintf(1,'\n');

% Seuraava pelkkää pluffia
h = waitbar(0,'Suodatetaan...');
for i=1:40,
  pause(0.03);
  waitbar(i/40,h)
end
close(h)

w = linspace(0, fs, length(Ysuod));
subplot(313);
plot(w, abs(Ysuod));
grid on;
axis([0 fs/2 0 1.1*max(abs(Y))]);
title('Suodatetun signaalin spektri');
subplot(312), xlabel('');
subplot(313), xlabel('taajuus (Hz)');

fprintf(1,'Kuunnellaan seuraavaksi suodatettu ääni.\n\n');

fprintf(1,'** Paina ENTER...\n\n');
pause;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;

soundsc(ysuod,fs);

fprintf(1,'Kuulostaako suodatettu signaali "paremmalta"?\n\n');

fprintf(1,'** Paina ENTER äänen kuuntelun jälkeen...\n\n');
pause;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Mitä on tapahtunut aikatasossa?!

clc

fprintf(1,'Piirretään vielä alkuperäinen ja suodatettu\n');
fprintf(1,'signaali aikatasoon\n\n');

figure(2), clf;
alku = 44800;
loppu = 45000;
plot(t(alku:loppu), y(alku:loppu));  
hold on;
ku3 = plot(t(alku:loppu), ysuod(alku:loppu), 'r');
set(ku3,'LineWidth',2);
title('Alkuperäinen ja suodatettu signaali x(t)');
xlabel('aika (sekunnit)');
ylabel('voimakkuus');
legend({'Alkuperäinen','Suodatettu'});
grid on;
axis([t(alku) t(loppu) 1.1*min(y(alku:loppu)) 1.1*max(y(alku:loppu))]);


fprintf(1,'** LOPPU. Paina ENTER.\n\n');
pause;

% Tällä voi vielä kirjoittaa signaalin ulos tiedostoon:
%wavwrite(ysuod, fs, nbits, 'esim1_output.wav');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%