% acomp.m 
% threshhold an audio signal using global DFT

% parameters 
% desired compression rate, approximated by % of Fourier coefficients to be kept
p = 70

% load sound, y is the sound vector, Fs is the sample rate

load('splat')
%load('train')

ly = length(y);
secs = ly/Fs;
y = y/max(abs(y)); % normalize

% get DFT of y
z=fft(y);
maz =  max(abs(z));

% plot y and |dft|
figure(1)  
clf
subplot(1,2,1)
plot((0:(ly-1))/Fs,y,'.r'); 
axis([0,secs,-1,1]);
title('original signal')
subplot(1,2,2)
plot(0:(length(z)-1),abs(z),'.r');
axis([0,ly,-0.1*maz,maz]);     
title('original |dft|');
drawnow
disp('Press a key for thresholded signals'), pause

% clip out small stuff
saz = sort(abs(z));
cutlev = saz(round((1-p/100)*length(saz))); 
J = find(abs(z) < cutlev);
cz = z;
cz(J) = 0;
cy = real(ifft(cz));

% plot thresholded  |dft| and reconstructed signal
figure(2)
clf
subplot(1,2,1)
plot(0:(length(z)-1),abs(cz),'.r');
hold on;
plot(0:(length(z)-1),cutlev*ones(size(cz)),'-y'); 
axis([0,ly,-0.1*maz,maz]);    
title('thresholded |dft|');
subplot(1,2,2)
plot((0:(ly-1))/Fs,cy,'.r');
axis([0,secs,-1,1]);
title('reconstructed signal');

% listen to the sounds
disp('Press a key for original sound'); 
pause
sound(y)
disp('Press a key for reconstructed sound'); 
pause
sound(cy)

% Compute compression and distortion
PercentRetainedCoefficients = 100*(ly-length(J))/ly
PercentDistortion  = 100*energy(y-cy)/energy(y)

% compute logarithm of cutoff level
  logcutlev = -log10(cutlev/maz)