%Katelyn Nye, EEN 502 - Fall 2007: Acoustics Project 3 
%1.
% Original Signal

[x,Fs,bits] = wavread('Speech');
t=[0:.5/26623:.5];
figure(1)
plot(t,x)
title('Time Waveform of Speech')
xlabel('Time (s)')
ylabel('Amplitude')

% Spectrogram plot of input signal obtained by Hann window & 
% converted to spectral representation via DFT.
[x,Fs,bits] = wavread('Speech');
N=2048; 
HannWindow=hanning(N);      
figure(2)
[s,F,t]=SPECTROGRAM(x,HannWindow,(N/2),N,Fs);
SPECTROGRAM(x,HannWindow,(N/2),N,Fs)
title('Spectrogram Plot of Speech')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 2.
% 25 critical band Bark spectrogram

FbandL = [0, 100, 200, 300, 400, 510, 630, 770, 920, 1080, 1270, 1480, 1720, 2000, 2320, 2700, 3150, 3700, 4400, 5300, 6400, 7700, 9500, 12000, 15500];
FbandL = transpose(FbandL);
FbandC = [50, 150, 250, 350, 450, 570, 700, 840, 1000, 1170, 1370, 1600, 1850, 2150, 2500, 2900, 3400, 4000, 4800, 5800, 7000, 8500, 10500, 13500, 19500];
FbandC=transpose(FbandC);
FbandR = [100, 200, 300, 400, 510, 630, 770, 920, 1080, 1270, 1480, 1720, 2000, 2320, 2700, 3150, 3700, 4400, 5300, 6400, 7700, 9500, 12000, 15500, 23500];
FbandR=transpose(FbandR);

window=zeros(24,1);
bark1=zeros(25,1);
bark2=zeros(25,1);
barkfreq1=zeros(1025,1);
barkfreq2=zeros(1025,1);

for j=1:24
    window(j,:)=FbandR(j,:)-FbandL(j,:);
end

figure(3)
[b,Fb,t]=SPECTROGRAM(x,HannWindow,(N/2),FbandC,Fs);
b=real(b);  %size: 25 by 13
SPECTROGRAM(x,HannWindow,(N/2),FbandC,Fs)
title('Bark Spectrogram Plot of Speech')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4.
% Power Spectrum

for i=1:13
    P(:,i)=real(s(:,i)).^2+imag(s(:,i)).^2;
end
Pdb=10*log10(P);

maxpow=0;

for i=1:1025
    for j=1:13
      maximum=Pdb(i,j);
      if maximum>maxpow
          maxpow=maximum;
          f1=i;
          t1=j;
      end
    end
end

onefourthmaxpow=10;

for i=1:1025
    for j=1:13
      if Pdb(i,j)>0 && Pdb(i,j)<onefourthmaxpow
          onefourthmaxpow=j;
          f2=i;
          t2=j;
      end
    end
end

freq = (1:N/2)/(N/2)*(Fs/2);  

figure(4)
plot(freq,Pdb(1:1024,t1))   %Power Spectrum at Maximum Signal Power
title('Maximum Signal Power of Signal Speech')
xlabel('Frequency (Hz)')
ylabel('SPL (dB)')

figure(5)
plot(freq,Pdb(1:1024,t2))    %Power Spectrum at 1/4th Maximum Signal Power
title('1/4th Maximum Signal Power of Signal Speech')
xlabel('Frequency (Hz)')
ylabel('SPL (dB)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Bark Spectrum for Maximum Signal Power

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j)
            bark1(j)=P(i,t1)+bark1(j);
        end
    end
end

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j);
        barkfreq1(i) = bark1(j);
        end
    end
end

barkfreq1dB=10*log10(barkfreq1);
figure(4)
hold on
plot(F,barkfreq1dB,'g') %Bark Spectrum at Maximum Signal Power

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Bark Spectrum for 1/4th Signal Power

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j)
            bark2(j)=P(i,t2)+bark2(j);
        end
    end
end

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j);
        barkfreq2(i) = bark2(j);
        end
    end
end

barkfreq2dB=10*log10(barkfreq2);
figure(5)
hold on
plot(F,barkfreq2dB,'g') %Bark Spectrum at 1/4th Signal Power

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Normalized Bark Spectrum for Maximum Signal Power

k=0;
for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j)
            bark1(j)=Pdb(i,t1)+bark1(j);
            k=k+1;
            count1(j)=k;
        end
    end
    k=0;
end

count1=transpose(count1);

for j=5:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j);
        barkoffset1(i) = bark1(j)/count1(j);
        end
    end
end

figure(4)
hold on
plot(F,barkoffset1,'y')  %Offset Bark Spectrum for Maximum Signal Power

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Normalized Bark Spectrum for 1/4th Signal Power

k=0;
for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j)
            bark2(j)=Pdb(i,t2)+bark2(j);
            k=k+1;
            count2(j)=k;
        end
    end
    k=0;
end

count2=transpose(count2);

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j);
        barkoffset2(i) = bark2(j)/count2(j);
        end
    end
end

figure(5)
hold on
plot(F,barkoffset2,'y') %Offset Bark Spectrum for 1/4th Signal Power

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Speading Function of Maximum Signal Power
spread_value=zeros(25,1);

for j=1:25
    for i=1:25
        spread_value(i,:)=sprdngf(j,i)+spread_value(i,:);
    end
end

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j);
        spread1(i,:) = spread_value(j,:).*bark1(j);
        end
    end
end

figure(4)
hold on
spread1=10*log10(spread1);
plot(F,spread1,'r')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Speading Function of 1/4th Signal Power

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j);
        spread2(i,:) = spread_value(j,:).*bark2(j);
        end
    end
end

figure(5)
hold on
spread2=10*log10(spread2);
plot(F,spread2,'r')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Offset Bark Spectrum of Maximum Signal Power

for i=1:1025
    geo1(i)=geomean(P(i,t1));
    art1(i)=mean(P(i,t1));
end

for i=1:1025
    SFM1(i)=10*log10(geo1(i)/art1(i));
    alpha1(i)=min((SFM1(i)/-60),1);
end

for i=1:1025
    O_i1(i)=alpha1(i)*(14.5+i)+(1-alpha1(i))*5.5;
end

O_i1=transpose(O_i1);

for i=1:1025
    Ti_1(i)=10^(log10(spread1(i))-(O_i1(i)/10));
end

figure(4)
hold on
plot(F,Ti_1,'c')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Offset Bark Spectrum of 1/4th Signal Power

for i=1:1025
    geo2(i)=geomean(P(i,t2));
    art2(i)=mean(P(i,t2));
end

for i=1:1025
    SFM2(i)=10*log10(geo2(i)/art2(i));
    alpha2(i)=min((SFM2(i)/-60),1);
end

for i=1:1025
    O_i2(i)=alpha2(i)*(14.5+i)+(1-alpha2(i))*5.5;
end

O_i2=transpose(O_i2);

for i=1:1025
    Ti_2(i)=10^(log10(spread2(i))-(O_i2(i)/10));
end

figure(5)
hold on
plot(F,Ti_2,'c')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Threshold of Human Hearing

figure(6)
    A=3.64*(F./1000).^(-.8);
    B=-6.5*exp((-.6)*((F./1000)-3.3).^2);
    C=10^(-3)*(F./1000).^4;
    D=db(1/(2^15.5));
    threshold=A+B+C+D;

plot(F,threshold,'r')
title('Absolute Threshold of Human Hearing')
xlabel('Frequency (Hz)')  %Put x-axis into Logarithmic scale
ylabel('SPL (dB)')

figure(8)
hold on
plot(F,threshold,'k')
plot(freq,Pdb(1:1024,t1))   %Power Spectrum at Maximum Signal Power
title('Maximum Signal Power of Signal Speech')
xlabel('Frequency (Hz)')
ylabel('SPL (dB)')

figure(9)
hold on
plot(F,threshold,'k')
plot(freq,Pdb(1:1024,t2))    %Power Spectrum at 1/4th Maximum Signal Power
title('1/4th Maximum Signal Power of Signal Speech')
xlabel('Frequency (Hz)')
ylabel('SPL (dB)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Bark Spectrum of Threshold of Human Hearing
bark_thres=zeros(25,1);
k=0;
for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j)
            bark_thres(j)=threshold(i)+bark_thres(j);
            k=k+1;
            count_thres(j)=k;
        end
    end
    k=0;
end

count_thres=transpose(count_thres);

for j=1:25
    for i=1:1025
        if F(i)>FbandL(j) && F(i)<FbandR(j);
        barkoffset_thres(i) = bark_thres(j)/count_thres(j);
        end
    end
end

figure(6)
hold on
plot(F,barkoffset_thres,'c') 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Final Threshold of Maximum Signal Power

for i=1:1025
   if Ti_1(i)>barkoffset_thres(i)
    abs_max1(i) = Ti_1(i);
   else
       abs_max1(i) = barkoffset_thres(i);
   end
end
abs_max1=abs_max1-20;
figure(8)
hold on
plot(F,abs_max1,'m')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Final Threshold of 1/4th Signal Power

for i=1:1025
   if Ti_2(i)>barkoffset_thres(i)
    abs_max2(i) = Ti_2(i);
   else
       abs_max2(i) = barkoffset_thres(i);
   end
end
abs_max2=abs_max2-20;
figure(9)
hold on
plot(F,abs_max2,'m')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Remaining Power Spectrum of Maximum Signal Power
g1=0;
for i=1:1025
    if Pdb(i,t1)>=abs_max1(i) 
        absolute1(i)=Pdb(i,t1);
        remaining1(i)=Pdb(i,t1);
        g1=g1+1;
    else
        absolute1(i)=nan;
        remaining1(i)=Pdb(i,t1)-300;
    end
end

percent_remain1=(g1/1025)*100

figure(8)
hold on
plot(F,absolute1,'c')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Remaining Power Spectrum of 1/4th Signal Power
g2=0;
for i=1:1025
    if Pdb(i,t2)>=abs_max2(i) 
        absolute2(i)=Pdb(i,t2);
         g2=g2+1;
    else
        absolute2(i)=nan;
    end
end

percent_remain2=(g2/1025)*100

figure(9)
hold on
plot(F,absolute2,'c')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Spectrogram of Remaining Power Spectrum of Maximum Signal Power
   
absolute1_ifft=ifft(remaining1,N);
WAVWRITE(absolute1_ifft,Fs,16,'SpeechCompressed')

[Q,Fs,bits] = wavread('SpeechCompressed');
t=[0:.5/2047:.5];
figure(10)
plot(t,Q)
title('Time Waveform of Compressed Speech')
xlabel('Time (s)')
ylabel('Amplitude')

figure(11)
SPECTROGRAM(Q)
title('Spectrogram Plot of Remaining Power Spectrum of Signal Speech')

%%%%%%%%%%%%%
%Speech: percent_remain1 = 8.9756, percent_remain2 = 25.8537
%Classical : percent_remain1 = 17.8537,  percent_remain2 = 12.3902
%Rock: percent_remain1 = 61.3659,  percent_remain2 = 58.3415