%% -------------------- PROJEKT: FREKVENCIJSKA KARAKTERISTIKA LJUDSKOG UHA ----------------------%%
%% -- CLANOVI TIMA: BRANIMIR DROPULJIC, NINA MAREVIC ------------------------------------------- %%

%% --- IDEJA: Ovaj program koristi se za realniju interpretaciju zvukova iz
%% okoline. Ljudsko uho ne cuje sve frekvencije jednako. Od raspona cujnih
%% frekvencija 20 - 20000 Hz, nize frekvencije razlikuje puno bolje od
%% visih. Frekvencije iznad 15 kHz gotovo i ne razlikuje, a jedva ih i
%% cuje. Racunalo nemoze interpretirati zvukove na nacin na koji to radi uho,
%% nego sve frekvencije, do fs/2 prikaze u jednakom rasponu. Za digitalnu obradu
%% govora vazno je da se zvuk prikaze na nacin na koji ga mi cujemo, da bi dobili
%% stvarnije i tocnije rezultaze. Zato je potrebno zvukove koje promatramo
%% predobraditi tj. prikazati ih onakvima, kakve ih uho stvarno i dozivljava.

%% Glavna funkcija koja se koristi u programu, je funkcija 'interpretacija_zvukova'.
%% Ona pokrece podfunkcije, i od ulaznog zvucnog signala radi realni
%% signal. Koristi podfunkciju 'filtri', koja napravi pojasne filtre koji su
%% odredjene sirine i odredjenog preklapanja (nisu idealno selektivni), te
%% su u odredjenoj skali ekvidistantni. Najcesce koristene skale su MEL i
%% BARK, koje postizu nelinearnu frekvencijsku os. Najnize frekvencije
%% imaju najvecu razlucivost, a najvise najmanju. Signal koji prolazi kroz
%% taj slog filtara predstavlja podatke o zvuku koji se dalje salje prema
%% mozgu pomocu neurona. Za prikaz signala na kompjuteru koristili smo
%% dvije mogu?nosti, slusanje signala, za koje je dovoljno ostaviti signale
%% takve kakvi jesu, i slikovno prikaz, za koje ih je potrebno obraditi.
%% Obrade se vrse podfunkcijama, 'log_anvelopa_en' i 'skraceni_kepstar'. Nakon
%% tog prikaza, signali se decimiraju i predstavljaju izlaze iz sklopa,
%% pogodne za daljnju obradu.

%% Na pocetku ove uvodne skripte - 'projekt' - koja dalje poziva glavnu
%% funkciju 'interpretacija_zvukova', odabire se zvucni signal koji se zeli
%% obraditi, oblik filtara, te skalu. Takodjer se odabire da li se zeli
%% prikazati slike filtara, te prikazati signali (zvukom i slikama) na
%% ulazu i izlazima iz kanala, i da li zelimo pohraniti izlaze iz kanala.
%% Ona nakon, sto pozove funkciju 'interpretacija_zvukova' i napravi sve sto
%% treba, prikazuje slike i zvukove te ih sprema ovisno o odabranom nacinu.

close all;
clear all;

% Odabir zvuka iz foldera 'zvukovi\'. 'zvuk_po_kanalima\' je folder u koji
% ce se kasnije spremati izlazi na kanalima.

pref='nina_sam';
zvukovi='zvukovi\';
zvuk_po_kanalima='zvuk_po_kanalima\';
zvuk=sprintf('%s%s',zvukovi,pref);

[y,FS,NBITS]=wavread(zvuk);
y=y(:);

% Odabir parametara:
%   M - broj pojasnih filtara [0:M] (uzima se oko 25 - 30)
%     - NAPOMENA!! - mora biti paran!
%   BW_f - koeficijent preklapanja filtara
%   Nimp - duzina imp odziva
%   CL - broj kepstralnih uzoraka kojih uzimamo
%   od_skala - odabrana skala

M=30;
BW_f=2;
Nimp=200;
CL=12;
od_skala=1;     % 1==MEL ,0==BARK
prikaz_slika_filtara=0;     %0==NE, 1==DA
prikaz_signala=0;           %0==NE, 1==DA
zapis_zvuka_po_kanalima=1;  %0==NE, 1==DA


% Poziv glavne funkcije, na ulaz joj dovodimo parametre filtara i ulazni
% signal. Na izlazu daje realizirane signale, spektre i impulsne odzive.

[skala,s_max,f,kor_dec,SP,SP_W,imp,imp_tr,imp_trw,ym,ym_en,ym_en_c,cep,ym_out,raz]=interpretacija_zvukova(y,FS,M,BW_f,Nimp,CL,od_skala);

% skala - odabrana skala
% s_max - njen maksimum
% f - frekvencijska skala
% kor_dec - odabrani korak decimacije
% SP - spektar odabrane skale
% SP_W - spektar odabrane skale nakon izgladjivanja impulsnih odziva
% imp - impulsni odziv
% imp_tr - odrezani impulsni odziv na 2*Nimp+1 uzoraka
% imp_trw - izgladjeni i skraceni impulsni odziv
% ym - matrica, signali na izlazu iz kanala
% ym_en - energija signala, nakon log_anvelopa_en funkcije
% ym_en_c - signali nakon primjene kratkog kepstra
% ym_out - decimirani ym_en_c - signal za dalje koristenje
% raz - razlika ym_en i ym_en_c


% Dalje slijede prikazi svih signala, imp odziva i filtara, u ovisnosti da
% li to zelimo. Takodjer naredbom soundsc slusanje izlaza iz pojedinih
% kanala, i snimanje '.wav' fileova u gore odabrani folder.

if prikaz_slika_filtara;
    
    figure(1)
    plot(f,skala);
    if od_skala;
        title('MEL SKALA');
    else
        title('BARK SKALA');
    end;
    xlabel('frekvencija [Hz]');
    axis([0 FS/2 0 s_max])
    grid;
    pause;

    figure(2)
    plot(skala,SP);
    if od_skala;
        title('MEL FILTRI');
    else
        title('BARK FILTRI');
    end;
    xlabel('frekvencija [Mel]');
    axis([skala(1,1) s_max 0 1])
    pause;
    
    plot(f,SP);
    if od_skala;
        title('MEL FILTRI');
    else
        title('BARK FILTRI');
    end;
    xlabel('frekvencija [Hz]');
    axis([0 FS/2 0 1])
    pause;
    
    figure(3)
    plot(f,abs(SP_W));
    if od_skala;
        title('MEL FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
    else
        title('BARK FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
    end;
    xlabel('frekvencija [Hz]');
    axis([0 FS/2 0 1])
    pause;

    figure(4)
    plot(skala,abs(SP_W));
    if od_skala;
        title('MEL FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
        xlabel('frekvencija [Mel]');
    else
        title('BARK FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
        xlabel('frekvencija [Bark]');
    end;
    axis([skala(1,1) s_max 0 1])
    pause;

    figure(5)
    plot(f,20*log10(abs(SP_W)));
    if od_skala;
        title('LOG-MEL FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
    else
        title('LOG-BARK FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
    end;
    xlabel('frekvencija [Hz]');
    ylabel('magnituda [db]')
    axis([skala(1,1) s_max -100 0])
    pause;
    
    plot(skala,20*log10(abs(SP_W)));
    if od_skala;
        title('LOG-MEL FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
        xlabel('frekvencija [Mel]');
    else
        title('LOG-BARK FILTRI NAKON IZGLADJENIH IMP. ODZIVA');
        xlabel('frekvencija [Bark]');
    end;
    ylabel('magnituda [db]')
    axis([skala(1,1) s_max -100 0])
    pause;

    figure(6)
    plot(skala,20*log10(sum(abs(SP_W)')/BW_f))
    ax=axis;
    axis([skala(1,1) s_max ax(3) ax(4)])
    title('Odstupanje sume karakteristika od ravne');
    ylabel('pogreska [dB]');
    pause;

    figure(7)
    plot([-Nimp:Nimp],imp_trw);
    title('IMPULSNI ODZIV IZGLADJENIH FILTARA')
    xlabel('2*Nimp+1')
    pause;

    figure(8)
    plot([-Nimp:Nimp],sum(imp_trw')-[zeros(1,Nimp) BW_f zeros(1,Nimp)],'r')
    title('Pogreska sume svih imp. odziva od jedinicnog impulsa');
    xlabel('2*Nimp+1')
    ylabel('pogreska');
    pause;

    figure(9)
    plot([-Nimp:Nimp],imp_trw(:,1));
    title('imp. odziv prvog filtra');
    pause;

    figure(10)
    plot([-Nimp:Nimp],imp_trw(:,M+1));
    title('imp. odziv zadnjeg filtra');
    pause;
    
end;


if prikaz_signala;

    figure(1+10*prikaz_slika_filtara)
    plot([1:length(y)]/FS,y);
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('SIGNAL KOJI OBRADJUJEMO')
    xlabel('vrijeme [sec]')
    soundsc(y,FS);
    pause;

    plot([1:length(y)]/FS,ym(:,5));
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('5. kanal signala')
    xlabel('vrijeme [sec]')
    soundsc(ym(:,5),FS);
    pause;

    plot([1:length(y)]/FS,ym(:,10));
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('10. kanal signala')
    xlabel('vrijeme [sec]')
    soundsc(ym(:,10),FS);
    pause;

    plot([1:length(y)]/FS,ym(:,15));
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('15. kanal signala')
    xlabel('vrijeme [sec]')
    soundsc(ym(:,15),FS);
    pause;

    plot([1:length(y)]/FS,ym(:,20));
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('20. kanal signala')
    xlabel('vrijeme [sec]')
    soundsc(ym(:,20),FS);
    pause;

    plot([1:length(y)]/FS,ym(:,25));
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('25. kanal signala')
    xlabel('vrijeme [sec]')
    soundsc(ym(:,25),FS);
    pause;

    plot([1:length(y)]/FS,sum(ym')')
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('Suma izlaza iz svih kanala')
    xlabel('vrijeme [sec]')
    soundsc(sum(ym')',FS);
    pause;

    plot([1:length(y)]/FS,[sum(ym')' y])
    ax=axis;
    axis([ax(1) length(y)/FS ax(3) ax(4)]);
    title('Suma izlaza iz svih kanala i ulazni signal (zeleni)')
    xlabel('vrijeme [sec]')
    pause;

    figure(2+10*prikaz_slika_filtara)
    plot([1:kor_dec:length(y)]/FS,ym_en(1:kor_dec:end,:));
    axis([0 length(y)/FS -100 -20]);
    title('Log-energija svakog pojasa');
    xlabel('vrijeme [sec]');
    ylabel('energija [dB]');
    pause;

    plot([1:kor_dec:length(y)]/FS,ym_en(1:kor_dec:end,15));
    axis([0 length(y)/FS -80 -30]);
    title('Log-energija 15. pojasa');
    xlabel('vrijeme [sec]');
    ylabel('energija [dB]');
    pause;

    pcolor([1:kor_dec:length(y)]/FS,[0:M],ym_en(1:kor_dec:end,:)');
    shading interp; 
    caxis([-80 -30]);
    colorbar;
    title('Log-energija svakog pojasa');
    xlabel('vrijeme [sec]');
    if od_skala;
        ylabel('indeks Mel-pojasa');
    else
        ylabel('indeks Bark-pojasa');
    end;
    pause;

    figure(3+10*prikaz_slika_filtara)
    plot([0:M],[ym_en(3000,:)' ym_en_c(3000,:)'])
    ax=axis;
    axis([0 M ax(3) ax(4)]);
    title('3000. uzorak, prije i nakon izgladjivanja kepstrom')
    pause;

    plot([1:length(y)]/FS,cep(:,1));
    axis([0 length(y)/FS -80 -30]);
    title('Nulti kepstralni koef cep[0]');
    xlabel('vrijeme [sec]');
    ylabel('[dB]');
    pause;

    pcolor([1:kor_dec:length(y)]/FS,[0:M],ym_en_c(1:kor_dec:end,:)');
    shading interp; 
    caxis([-80 -30]);
    colorbar;
    title('Log-energija svakog pojasa - nakon skracenog spektra');
    xlabel('vrijeme [sec]');
    if od_skala;
        ylabel('indeks Mel-pojasa');
    else
        ylabel('indeks Bark-pojasa');
    end;
    pause;

    figure(4+10*prikaz_slika_filtara)
    pcolor([1:kor_dec:length(y)]/FS,[0:M],raz(1:kor_dec:end,:)');
    shading interp;
    caxis([-10 10]);
    colorbar;
    title('Razlika Log-en spektra prije i nakon izgladjivanja kratkim kepstrom');
    xlabel('vrijeme [sec]');
    if od_skala;
        ylabel('indeks Mel-pojasa');
    else
        ylabel('indeks Bark-pojasa');
    end;
    pause;

end;

if zapis_zvuka_po_kanalima;
    
    for i=1:M+1,
        zvuk_kanal=sprintf('%s%s_%d',zvuk_po_kanalima,pref,i);
        wavwrite(ym(:,i)/max(abs(ym(:,i))),FS,NBITS,zvuk_kanal);
    end;
    zvuk_kanal=sprintf('%s%s_%s',zvuk_po_kanalima,pref,'sum');
    wavwrite(sum(ym')',FS,NBITS,zvuk_kanal);
    
end;
