file = 'sini1k.wav';
file = 'saha1k.wav';
file = 'kantti1k.wav';
file = 'moottoritie.wav';

siz = wavread(file,'size')

% Huomioita diskreettiaikaisesta signaalinkäsittelystä ja Fourier
% muunnoksesta
% Näytteenottotaajuuden (Fs) puolikkaaseen asti vain voi olla informaatiota
% tätä sanotaan Nyquist taajuudeksi, joka on Fs/2
%
% Aikatasosignaali: eli aikasarja tietyllä näytteenottotaajuudella
% Taajuustaso: Fourier muunnos muuntaa aikatasosignaalin taajuustasoon
% Aikatasosignaali merkitään x ja taajuustasoa X
% Numeerisessa laskennassa tehdään diskreetti Fourier muunnos (DFT, hidas: O(N^2))
% käytännössä muunnos tehdään FFT algoritmilla (Fast Fourier Transform,
% nopea: O(N log N) )
% FFT algoritmi (yleensä) saa syötteeksi aikasarjan ja muunnospisteiden
% lukumäärän
% Muunnospisteet (n) ovat 2^n
% eli esim 2^8 = 256 tuottaa 256 pisteisen muunnoksen, josta käytännössä
% voi hyödyntää 128 pistettä
% Muunnoksen resoluutio riippuu näytteenottotaajuudesta ja muunnospisteiden
% lukumäärästä (Fs/n). Resoluutioväliä sanotaan bin:iksi. Esim FFT bin =
% Fs=48000 Hz, n=256 = 48000/256 = 187.5 Hz eli yhteen muunnospisteeseen
% tulee kaikki väliltä 187.5 Hz oleva informaatio.
%
% Muunnoksen voi tehdä myös toisinpäin eli taajuustasosta aikatasoon
% tällöin muunnoksen nimi on käänteinen diskreetti Fourier muunnos (inverse
% DFT). Käytännössä iFFT algoritmi
% Huom! muunnoksessa otetaan kerrallaan signaalia muunnospisteiden
% lukumäärän verran, jos on pitkä aikasarja niin muunnoksen voi tehdä koko
% signalille monin eri tavoin
% esim:
% 1) ottamalla n pisteisen muunnoksen ja ottamalla tuloksen talteen. 
% Seuraava muunnos tehdään ottamalla signaalia alkaen edellisestä n
% pisteestä n pisteiden verran ja summaamalla edelliseen tulokseen ja
% seuraava vaihe samalla tavalla kunnes signaali on käyty loppuun.
% Lopuksi jakamalla muunnoksien lukumäärällä, jolloin tuloksena on n
% pisteinen keskiarvospektri
% 2) tekemällä n pisteisiä muunnoksia kuten edellä, mutta keräämällä
% tulokset matriisin, jolloin tuloksena on spektrogrammi

%% Esimerkki, josta näkee miten FFT algortimi tuottaa symmetrisen tuloksen
% FFT tuloksen ensimmäisestä pisteestä puoleenväliin asti on hyödylliset tulokset.
[data1,Fs,NBITS]=wavread(file);
n=4096;
X=fft(data1,n);
X=abs(X);
plot(X);

%% Esimerkki keskiarvospektristä
n=4096; %

X_avg=[n];
j=0;
for i=1:n:(siz-n)
    [data1,Fs,NBITS]=wavread(file,[i i+n]);
    X=fft(data1,n);
    X_avg=X_avg+abs(X); % Matlab abs osaa laskea absoluuttisen arvon kompleksiluvuista. sqrt(real^2 + imaginary^2);
    j=j+1;
end
X = X_avg./j;

% FFT bin = Fs/n
hz_labels=[];
hz=Fs/n;
for i=1:(n/2)
    hz_labels(i)=hz*i;
end
hz_labels = hz_labels';

% piirretään tulokset

% plot waveform
subplot(2,1,1);
temp=data1(1:200);
t=((0:length(temp)-1)/Fs)*1000;        % times of sampling instants
plot(t,temp);
xlabel('Time (ms)');
ylabel('Amplitude');
title('Waveform');

% plot spectrum 
subplot(2,1,2),
plot(hz_labels,20*log10(X(1:n/2))) % tehospektri 20*log10(X)
%set(gca,'XTick',[1:1000:Fs/2]); 		
xlabel('Hz'); 
ylabel('dB/Hz');
title('FFT') 