function demo2
% L170101

[y, fs, nbits] = wavread('sampsa.wav');
y = y - mean(y);
fprintf(1,'Seuraava äänitiedosto on tutkija Sampsa Laineen lahja maailmalle...\n');
soundsc(y, fs);
pause;

pituus_sek = length(y)/fs;
t = linspace(0, pituus_sek, length(y));
plot(t,y);
title('sampsa.wav: a-äänne');
xlabel('aika (sekunnit)');
ylabel('voimakkuus');
% Näyttää varsin sekavalta - otetaan pieni näyte kymmenysosasekunnin 
% ajalta:
fprintf(1,'Näyttää varsin sekavalta - otetaan pieni näyte...\n');
pause;

plot(t(1:4410),y(1:4410));
grid on;
title('sampsa.wav, t = 0..0.1 sek');
xlabel('aika (sekunnit)');
ylabel('voimakkuus');
% Signaali näyttää lähes periodiselta eli x(t) = x(t+T). 
% Mikä on sen "jaksonpituus" T?
fprintf(1,'Signaali näyttää lähes periodiselta eli x(t) = x(t+T).\n');
fprintf(1,'Mikä on sen "jaksonpituus" T?\n');
pause;

% Taajuustasoon fourier-muuntamalla ja sen spektri:
fprintf(1,'\n\nSignaaleja kannattaa tutkia niin aika- kuin taajuustasossa.\n');
fprintf(1,'Taajuustason tutkimiseen Fourier-muunnos on oiva työkalu.\n');
yF = fft(y);
hz = linspace(0,fs, length(yF));
plot(hz(1:8917), abs(yF(1:8917)));
grid on;
title('sampsa.wav, spektri f = 0..4000 Hz');
xlabel('taajuus (hertz eli 1/s)');
ylabel('amplitudi');
fprintf(1,'Jaksollisuus näkyy nyt piikkeinä spektrissä!\n');
fprintf(1,'Arvioi ensimmäisen piikin taajuus (perustaajuus)\n');
fprintf(1,'ja vertaa sitä yllä laskemaasi jaksonpituuteen T!\n');
pause;

% Poistetaan ensimmäiset piikit ylipäästösuodatuksella
fprintf(1,'\nSuodatetaan suurin osa pienistä piikeistä pois ylipäästösuodatuksella\n');
pause;
% KOMMENTTI: ensimmäiset peruskomponentit eivät vaikuttaneet juuri mitään,
%  onko energia sittenkin keskittynyt 2-3 kHz:iin?
Rs = 2*2200/fs;
Rp = 2*3000/fs;
[N,Wn] = ellipord(Rp, Rs, 1, 40);
[b,a] = ellip(N, 1, 40, Wn, 'high');
%freqz(b,a,1024,fs);
y2 = filter(b,a,y);

fprintf(1,'...ja kuunnellaan lopputulos\n');
%spektristä katoaa piikit
pause;

soundsc(y2, fs);
figure,clf,
plot(t(1200:1800),y(1200:1800),'b',t(1200:1800),y2(1200:1800),'r');
legend('Alkuperäinen','Suodatettu');

YY=fft(y2);
plot(hz,abs(YY))