close all;

% ***************get a plot of Hz vs. Bark***************
%f=[1:1:20000];
%b=hz2bark(f);
%plot(f,b);
%xlabel('Frequency (Hz)');
%ylabel('Frequency (Bark)');
%title('Relationship between Hertz and Bark Frequencies');

% ***************open a file***************
[y, fs, nbits]=wavread('wavfiles/sub1.wav');

% ***************necessary constants***************
FFTlength=512;                     % FFT length, naturally
f=[1:FFTlength/2]*(fs/FFTlength);  % array of Hz corresponding to bins
b=hz2bark(f);                      % array of Bark corresponding to bins
% get absolute threshold of hearing with given freq.
ATH=3.64.*(f./1000).^(-.8) - 6.5.*exp(-0.6.*(f./1000-3.3).^2) + 0.001.*(f./1000).^4;
%figure;
%plot(f,ATH);
%xlabel('Frequency (Hz)');
%ylabel('SPL (dB)');
%title('Absolute Threshold of Hearing (Hz scale)');
%figure;
%plot(b,ATH);
%xlabel('Frequency (Bark)');
%ylabel('SPL (dB)');
%title('Absolute Threshold of Hearing (Bark scale)');

% ***************practice run***************
y=y(86001:86256);

figure;
plot(y);
xlabel('Time (sample)');
ylabel('Magnitude');
title('sub1.wav');

% step 1: normalize the signal
%y=normalize(y, FFTlength, nbits);
%figure;
%plot(y);
%xlabel('Time (sample)');
%ylabel('Magnitude');
%title('sub1.wav NORMALIZED');

% step 2: get power density spectrum
powernormalizationconstant=90.302;
power=psd(y, FFTlength, powernormalizationconstant);
% only need first half since signal is real
power=power(1:256);

figure;
subplot(2,1,1);
plot(f, power);
hold on;
plot(f, ATH, 'r--');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Power Density Spectrum of sub1.wav');

subplot(2,1,2);
plot(b, power);
hold on;
plot(b, ATH, 'r--');
xlabel('Frequency (Bark)');
ylabel('Magnitude');
title('Power Density Spectrum of sub1.wav');

% step 3: find tones
tone_psd=find_tones(power);

figure;
toneloc=find(tone_psd);
subplot(2,1,1);
stem(f(toneloc), tone_psd(toneloc), 'x');
hold on;
plot(f, power, 'b');
plot(f, ATH, 'r--');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Tone Maskers');

subplot(2,1,2);
stem(b(toneloc), tone_psd(toneloc), 'x');
hold on;
plot(b, power, 'b');
plot(b, ATH, 'r--');
xlabel('Frequency (Bark)');
ylabel('Magnitude');
title('Tone Maskers');

% step 4: find noise maskers within critical band
% start with 0-1
noise_psd=zeros(1,length(tone_psd));
lowbin=1;
highbin=max(find(b<1));
% remainder of critical bands can be done with loop
for band = 1:24,
   [noise_psd_at_loc, loc]=noise_masker(power, tone_psd, lowbin, highbin);
   if (loc ~= -1)
      noise_psd(floor(loc))=noise_psd_at_loc;
   end
   lowbin=highbin;
   highbin=max(find(b<(band+1)));
end

figure;
noiseloc=find(noise_psd);
subplot(2,1,1);
stem(f(noiseloc), noise_psd(noiseloc), 'o');
hold on;
plot(f, power, 'b');
plot(f, ATH, 'r--');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Noise Maskers');

subplot(2,1,2);
stem(b(noiseloc), noise_psd(noiseloc), 'o');
hold on;
plot(b, power, 'b');
plot(b, ATH, 'r--');
xlabel('Frequency (Bark)');
ylabel('Magnitude');
title('Noise Maskers');

% step 5: remove maskers that are near each other
[tone_psd, noise_psd]=check_maskers(tone_psd, noise_psd, ATH, b);

figure;
toneloc=find(tone_psd);
subplot(2,1,1);
stem(f(toneloc), tone_psd(toneloc), 'x');
hold on;
plot(f, power, 'b');
plot(f, ATH, 'r--');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Checked Tone Maskers');

subplot(2,1,2);
stem(b(toneloc), tone_psd(toneloc), 'x');
hold on;
plot(b, power, 'b');
plot(b, ATH, 'r--');
xlabel('Frequency (Bark)');
ylabel('Magnitude');
title('Checked Tone Maskers');

figure;
noiseloc=find(noise_psd);
subplot(2,1,1);
stem(f(noiseloc), noise_psd(noiseloc), 'o');
hold on;
plot(f, power, 'b');
plot(f, ATH, 'r--');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Checked Noise Maskers');

subplot(2,1,2);
stem(b(noiseloc), noise_psd(noiseloc), 'o');
hold on;
plot(b, power, 'b');
plot(b, ATH, 'r--');
xlabel('Frequency (Bark)');
ylabel('Magnitude');
title('Checked Noise Maskers');

% step 6: now that we have important maskers, calculate global threshold
thres=global_threshold(ATH, tone_psd, noise_psd, b);

figure;
plot(b,thres);
hold on;
plot(b(find(tone_psd)), tone_psd(find(tone_psd)), 'x');
plot(b(find(noise_psd)), noise_psd(find(noise_psd)), 'o');
plot(b, power, 'g');
plot(b, ATH, 'r--');
xlabel('Frequency (Bark)');
ylabel('Magnitude');
title('Global Masking Threshold');
legend('Overall Threshold', 'Tones', 'Noise', 'Original PSD', 'ATH');