%	dft_chop
%	show effect of chopping DFT coefficients

%file = 'handel';
file = 'train.mat';	% located in matlab-6.5/toolbox/matlab/audio
dat = load(file);	% load .wav file
x = dat.y;		% samples of audio signal
fs = dat.Fs;		% sampling frequency

clf, pl = 420;
N = length(x);
n = [1:N];
t = n / fs;

subplot(pl+1), plot(n, x), axis([0 N -1 1])
ylabel x[n], %xlabel n
title(sprintf('Original audio signal, Fs=%g', fs))

subplot(pl+3), plot(n, x), axis([8000 8200 -1 1])
ylabel x[n]

X = (1/N) * fft(x);	% compute DFT of x using FFT

subplot(pl+2), plot(0:[N-1], abs(X)), axis([0 N 0 0.08])
xlabel k, ylabel |X[k]|, title 'Magnitude spectrum'

subplot(pl+4), plot(0,0), axis off
text(0,0,sprintf('%d nonzero DFT coefficients', sum(X~=0)), 'horiz', 'center')

% YOU WILL EXPERIMENT WITH THE VALUE BELOW
Y = X;				% copy DFT coefficients
threshold = 0.08 * max(abs(X));	% try 8% of maximum DFT coefficient
Y(abs(X) < threshold) = 0;	% zero out the "small" DFT coefficients

subplot(pl+6), plot(0:[N-1], abs(Y)), axis([0 N 0 0.08])
xlabel k, ylabel |Y[k]|, title 'Sparsified magnitude spectrum'

% synthesize signal from modified spectrum using inverse DFT via IFFT
y = real(N * ifft(Y));

subplot(pl+5), plot(n, y), axis([0 N -1 1])
ylabel y[n], title 'Decompressed audio signal'

subplot(pl+8), plot(0,0), axis off
text(0,0,sprintf('%d nonzero DFT coefficients', sum(Y~=0)), 'horiz', 'center')
text(0,-0.3,sprintf('with |X[k]| > %g', threshold), 'horiz', 'center')

subplot(pl+7), plot(n, y), axis([8000 8200 -1 1])
xlabel n, ylabel y[n]

% compute MSE between original "x" and reconstructed "y"
%%% MSE = ?		% PUT YOUR CODE HERE

% combine the "before" and "after" signals for listening
z = [zeros(500,1); x; zeros(1000,1); y];
sound(z)
%wavwrite(z, fs, 8, 'dft_chop.wav');
%Print('dft_chop')
