clc
clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                                        %% Session 4 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=512;

fid = fopen('waziwaza.raw'); %fid negative mean the file couldn't be open
wave = fread(fid,'short');
soundsc(wave,16000); 
figure(1)
N=17100;
t=[1/16000:1/16000:N/16000];
sigma2=25*10^6;
noise = sqrt(sigma2)*randn(N,1); %Gaussian Centred White Noise with sigma2's variance
wave_noised = wave + noise; % noised wave
subplot 311
plot(t,wave);
xlim([0 1]);
subplot 312
plot(t,noise);
xlim([0 1]);
subplot 313
plot(t,wave_noised);
xlim([0 1]);
xlabel('time in second');
ylabel('GCWN');
title('Gaussian Centred White Noise with sigma2''s variance');


wavered= wave_noised([2100:2611]);
wave_windowed=wavered.*hann(M);

N=512;

for (k=2100:1:2611)
    S=10*log10(abs(fft(wave_windowed)).^2./(N));
end


load f2bark.mat;  

figure(2)

subplot(2,2,1)
plot([0:M/2-1]*16000/M,S(1:256));
xlabel(' Frequency (Hz)');
ylabel('Magnitude  (dB)');
title('S(k)');


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%% Periodogram with critic band %%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,2)

hold on
plot([0:M/2-1]*16000/M,S(1:256));
xlabel(' Frequency (Hz)');
xlim([0 8000]);
ylabel('Magnitude en dB');
title('S(k) with Bark scale');
for k=1:24
    line([f2bark(k,2),f2bark(k,2)],[0,100],'color','r');
end
legend('S(k)','bark');
hold off

                %%%%%%%%%%%%%%%%%%  
                %%%%% Step 3 %%%%%
                %%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%% Signal Component Nature Determination %%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% Tonale and untonale frequencies %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



ton=0;a=1;b=1; % a pour les tonales et b pour les atonales
for k=4:251
    if ( (S(k)>S(k-1)) && (S(k)>=S(k+1)) )        
       if ( k>2 && k<63 )
            for jj=[-2,2]
                if ((S(k)-S(k-jj))>=7)
                    ton=ton+1;
                end
            end
            if ton==2
                Tonale(a,1)=k-1;
                a=a+1;
                ton=0;
            else
                Tonale(b,2)=k-1;
                b=b+1;
                ton=0;
            end            
       end
        if ( k>= 63 && k<127 )
            for jj=[-3,-2,2,3]
                if ((S(k)-S(k-jj))>=7)
                    ton=ton+1;
                end  
            end
            if ton==4
                Tonale(k,1)=k-1;
                a=a+1;
                ton=0;
             else
                Tonale(b,2)=k-1;
                b=b+1;
                ton=0;
            end            
        end
        if ( k>=127 && k<250)
            for jj=[-6:-2,2:6]
                if  ((S(k)-S(k-jj))>=7)  
                    ton=ton+1;
                end
            end
             if ton==10
                a=a+1;
                ton=0;
              else
                Tonale(b,2)=k-1;
                b=b+1;
                ton=0;
             end
        end
    
    else
        Tonale(b,2)=k-1;
        b=b+1;
    end
end
clear ton a b jj k;


figure(2)

subplot(2,2,3:4)
hold on
plot([0:M/2-1]*16000/M,S(1:256));
xlabel(' Frequency (Hz)');
xlim([0 8000]);
ylabel('Magnitude en dB');
title('S(k) with tonale plotted');

for k=1:24
    line([f2bark(k,2),f2bark(k,2)],[0,100],'color','r');
end
for k=1:length(Tonale)
    if Tonale(k)~=0
    plot((Tonale(k))*16000/M, S(Tonale(k)+1),'*','color','green')
    end
end
legend('S(k)','bark');
hold off

figure(3)
hold on
plot([0:M/2-1]*16000/M,S(1:256));
xlabel(' Frequency (Hz)');
xlim([0 2000]);
ylabel('Magnitude en dB');
ylim([60 100]);
title('S(k) with tonale plotted');


for k=1:24
    line([f2bark(k,2),f2bark(k,2)],[0,100],'color','r');
end
for k=1:length(Tonale)
    if Tonale(k)~=0
    plot((Tonale(k))*16000/M, S(Tonale(k)+1),'*','color','green')
    end
end
legend('S(k)','bark');
hold off


                %%%%%%%%%%%%%%%%%%  
                %%%%% Step 4 %%%%%
                %%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%% Component Treatment %%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % %%% %%% %%% %% Tonal Treatment %% %%% %%% %%%  % 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Fech = 16000;
for k=1:length(Tonale)
    if Tonale(k,1)~=0
        P1(k,1)=10*log10(10^(S(Tonale(k,1))/10)+10^(S(Tonale(k,1)+1)/10)+10^(S(Tonale(k,1)+2)/10));
    end
end
clear k;
%Sa est le seuil d'audition absolue
f=[0.001:0.001:20];
Sa=3.64*f.^(-0.8)-6.5*exp(-0.6*(f-3.3).^2)+10^(-3)*f.^4;

%suppression des tonales sous le seuil d'audition absolue
for k=1:length(Tonale)
    if Tonale(k,1)~=0
        if P1(k)>Sa(1,floor((Tonale(k,1)-1)*16000/M )) %floor renvoie la partie enti鑢e
            TonaleTraited(k,1)=Tonale(k);
            TonaleTraited(k,2)=P1(k);
        end
    end
end
clear k;


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % %%% %%% %%% %% Non-Tonal Treatment %% %%% %%% %%%  % 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a=0; %a pour les cases du tableau
b=0; %b pour les barks 

a=0;
for b=1:24
    for k=1:length(Tonale)
        if (bark((Tonale(k,2))*16000/M)==b)
            a=a+10^(S( Tonale(k,2) )/10); 
            TonaleTraited(b,4)=10*log10(a);
            TonaleTraited(b,3)=f2bark(b,3);
        end
    end
    a=0;
end
clear a b k;


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%% Tonal et Untonal plotted after treatment %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(4)

subplot 211

hold on
plot([0:M/2-1]*16000/M,S(1:256));
xlim([0 8000]);
xlabel(' Frequency (Hz)');
ylabel('Magnitude en dB');
title('S(k)with tonale and non-tonale after treatment');

for k=1:24
    line([f2bark(k,2),f2bark(k,2)],[0,100],'color','red');
end
for k=1:length(Tonale)
    if Tonale(k)~=0
    plot((TonaleTraited(k))*16000/M, TonaleTraited(k,2),'*','color','green')
    end
end
for k=2:length(TonaleTraited)
    plot(TonaleTraited(k,3),TonaleTraited(k,4) ,'*','color','blue');
end

hold off


subplot 212
xlim([0 2000]);
ylim([20 100]);
hold on
plot([0:M/2-1]*16000/M,S(1:256));
xlabel(' Frequency (Hz)');
ylabel('Magnitude en dB');

for k=1:24
    line([f2bark(k,2),f2bark(k,2)],[0,100],'color','red');
end
clear k;
for k=1:length(Tonale)
    if Tonale(k)~=0
    plot((TonaleTraited(k))*16000/M, TonaleTraited(k,2),'*','color','green')
    end
end
for k=2:length(TonaleTraited)
    plot(TonaleTraited(k,3),TonaleTraited(k,4) ,'*','color','blue');
end

hold off

                %%%%%%%%%%%%%%%%%%  
                %%%%% Step 5 %%%%%
                %%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%% Component's Influence on nearby components %%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Voir la fonction dans le fichier masque.m


                %%%%%%%%%%%%%%%%%%  
                %%%%% Step 6 %%%%%
                %%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%% Masked Curve Calculation %%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N=512;

for (f=1:256) 
    Masquage_Courbe(f,1)=MasqueTotal(f*16/512);
    
end

figure (5)

subplot 211

hold on
plot([0:M/2-1]*16000/M,S(1:256));
plot([1:256]*16000/N,Masquage_Courbe(1:256),'color','black');
xlabel(' Frequency (Hz)');
xlim([0 8000]);
ylabel('Magnitude en dB');
ylim([0 100]);
title('Masking curve of S(k)');

for k=1:24
    line([f2bark(k,2),f2bark(k,2)],[0,100],'color','r');
end
hold off

subplot 212


hold on
plot([0:M/2-1]*16000/M,S(1:256));
plot([1:256]*16000/N,Masquage_Courbe(1:256),'color','black');
xlabel(' Frequency (Hz)');
xlim([0 2000]);
ylabel('Magnitude en dB');
ylim([20 100]);

for k=1:24
    line([f2bark(k,2),f2bark(k,2)],[0,100],'color','r');
end
clear k;
for k=1:length(Tonale)
    if Tonale(k)~=0
    plot((TonaleTraited(k))*16000/M, TonaleTraited(k,2),'*','color','green');
    end
end
hold off