clear all; close all; clc;  

%% NACTENI PROMLUV

addpath Funkce
addpath Promluvy

for i=1:2
    eval(['[speech',num2str(i),'] = wavread(''',num2str(i),'_a.wav'');']); 
end

%% MFCC

Tw = 25;                % analysis frame duration (ms)
Ts = 10;                % analysis frame shift (ms)
alpha = 0.97;           % preemphasis coefficient
M = 13;                 % number of filterbank channels 
C = 12;                 % number of cepstral coefficients
L = 22;                 % cepstral sine lifter parameter
LF = 0;                 % lower frequency limit (Hz)
HF = 8000;              % upper frequency limit (Hz)
fs=16000;               % frequency

Nw = round( 1E-3*Tw*fs );    % frame duration (samples)
Ns = round( 1E-3*Ts*fs );    % frame shift (samples)

nfft = 2^nextpow2( Nw );     % length of FFT analysis 
K = nfft/2+1;                % length of the unique part of the FFT 

for i=1:2
    eval(['[ MFCCs' num2str(i) ', FBEs' num2str(i) ', frames' num2str(i) ' ] = mfcc( speech' num2str(i) ', fs, Tw, Ts, alpha, @hamming, [LF HF], M, C+1, L );']);
end

%% DYNAMICKE BORCENI CASU

MFPR = MFCCs2;
MFS = 2;

n2=length(MFPR);
 
for G=1:2

    if G == MFS
       continue
    end
  
    eval(['n1=length(MFCCs' num2str(G) ');']);

    eval(['SM' num2str(G) ' = simmx(abs(MFCCs' num2str(G) '),abs(MFPR));']);

    eval(['[p,q,CCC] = dp(1-SM' num2str(G) ');']);
    
    if ((n1>n2) || (n1<n2))
        
        ozq=zeros(1,length(q));
        ozp=zeros(1,length(p));
        countp=0;
        ddd=0;dddp=0;cdc=0;cdp=0;
        countq=0;

        for i=1:(length(p)-1)
            if p(i) ~= p(i+1)
                continue        
            elseif p(i) == p(i+1)
                ozp(i)=1;
                countp=countp+1;
                dddp=dddp+1;
                if dddp==1
                    cdp=i;
                elseif dddp>1
                    cdp=[cdp,i];
                end
            end
        end

        if countp~=0
            for i=1:(C+1)
                for j=1:(n1+countp)
                    eval(['new' num2str(G) '(i,j) = FBEs' num2str(G) '(i,p(j));']);
                end
            end
        end

        for i=1:(length(q)-1)
            if q(i) ~= q(i+1)
                continue        
            elseif q(i) == q(i+1)
                ozq(i)=1;
                countq=countq+1;
                ddd=ddd+1;
                if ddd==1
                    cd=i;
                elseif ddd>1
                    cd=[cd,i];
                end
            end
        end

        eval(['new' num2str(G) '(:,cd)=[];']);

    elseif n2>n1
         
        ozq=zeros(1,length(q));
        ozp=zeros(1,length(p));
        countp=0;
        ddd=0;dddp=0;cdc=0;cdp=0;
        countq=0;

        for i=1:(length(q)-1)
            if q(i) ~= q(i+1)
                continue        
            elseif q(i) == q(i+1)
                ozq(i)=1;
                countq=countq+1;
                ddd=ddd+1;
                if ddd==1
                    cd=i;
                elseif ddd>1
                    cd=[cd,i];
                end
            end
        end

        for i=1:(C+1)
            for j=1:((n1+countq)-9)
                eval(['new' num2str(G) '(i,j) = FBEs' num2str(G) '(i,q(j));'])
            end
        end
        
        for i=1:(length(p)-1)
            if p(i) ~= p(i+1)
                continue        
            elseif p(i) == p(i+1)
                ozp(i)=1;
                countp=countp+1;
                dddp=dddp+1;
                if dddp==1
                    cdp=i;
                elseif dddp>1
                    cdp=[cdp,i];
                end
            end
        end

        if countp ~= 0
            eval(['new' num2str(i) '(:,cdp)=[];'])
        end            
    end    
end

%% VYPOCET NOVYCH ENERGII
new2 = FBEs2;
POD = new1./new2;
PRUM = mean(POD,2);

NO = PRUM/max(PRUM);
INVNO = 1./NO;

%% TRIFBANK

hz2mel = @(hz)(1127*log(1+hz/700));     % Hertz to mel warping function
mel2hz = @(mel)(700*exp(mel/1127)-700); % mel to Hertz warping function

[ H1, freq ] = trifbank( M, K, [0 8000], fs, hz2mel, mel2hz );

%% ROZDELENI FILTRU

count=0;
for i=1:M
    for y=1:length(freq)
         if H1(i,y) ~= 0
             eval(['trif' num2str(i) '(1,y) = freq(1,y);']);
             eval(['trif' num2str(i) '(2,y) = H1(i,y);']);
         end
    end
    eval(['trif' num2str(i) '(:,all(~trif' num2str(i) ',1)) = [];']);
end

%% LOW SHELVING 1

G = mag2db(abs(INVNO(1)));
if G ~= 0
    fc = round(mean(trif1,2));
    fc = fc(1);
    Q = 150;
    [b, a] = shelving(G, fc, fs, Q, 'Base_Shelf');
    peak = filter(b,a,speech1);
end

%% PEAK2

G = mag2db(abs(INVNO(2)));
if G ~= 0
    fc = round(mean(trif2,2));
    Q = 25;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK3

G = mag2db(abs(INVNO(3)));
if G ~= 0
    fc = round(mean(trif3,2));
    Q = 30;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK4

G = mag2db(abs(INVNO(4)));
if G ~= 0
    fc = round(mean(trif4,2));
    Q = 25;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK5

G=mag2db(abs(INVNO(5)));
if G ~= 0
    fc = round(mean(trif5,2));
    Q = 30;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK6

G = mag2db(abs(INVNO(6)));
if G ~= 0
    fc = round(mean(trif6,2));
    Q = 28;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK7

G = mag2db(abs(INVNO(7)));
if G ~= 0
    fc = round(mean(trif7,2));
    Q = 15;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK8

G = mag2db(abs(INVNO(8)));
if G ~= 0
    fc = round(mean(trif8,2));
    Q = 30;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end


%% PEAK9

G = mag2db(abs(INVNO(9)));
if G ~= 0
    fc = round(mean(trif9,2));
    Q = 20;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK10

G=mag2db(abs(INVNO(10)));
if G ~= 0
    fc = round(mean(trif10,2));
    Q = 18;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK11

G = mag2db(abs(INVNO(11)));
if G ~= 0
    fc = round(mean(trif11,2));
    Q = 15;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% PEAK12

G = mag2db(abs(INVNO(12)));
if G ~= 0
    fc = round(mean(trif12,2));
    Q = 6;
    [b, a] = peaking(G, fc(1), Q, fs);
    peak = filter(b,a,peak);
end

%% HIGH SHELVING 13

G = mag2db(abs(INVNO(13)));
if G ~= 0
    fc = round(mean(trif13,2));
    fc = fc(1);
    Q = 100;
    [b, a] = shelving(G, fc, fs, Q, 'Treble_Shelf');
    speechhigh = filter(b,a,peak);
end

%% EQUALIZACE

Ei=(sum(speech1.^2)/length(speech1));
E0=(sum(speechhigh.^2)/length(speechhigh));

yyy=speechhigh.*(Ei/E0);


%% ZMENA VYSKY HLASU

shiftOffset = 0.7; 
DAFx1 = FunctionDirectPitchShift(yyy, fs, shiftOffset);

