% stfttest.m
%
% An example Matlab script to demonstrate block FFT processing of an input
% sound and subsequent reconstruction via the IFFT and overlap adding.
%
% Gary Scavone, McGill University, 2016.

clear;
[y, fs] = audioread( 'test.wav' );
N = length( y );

nFFT = 1024;
hop = floor( nFFT / 4 );
nBlocks = floor( N / hop ) - 1; % approximately
F = zeros( nFFT, nBlocks );
w = hanning( nFFT );

for n = 1:nBlocks
  iStart = (n-1) * hop;
  if iStart+nFFT > N, break; end
  F(:, n) = fft( w .* y(iStart+1 : iStart+nFFT) );
end

% One frequency-domain processing option ... lowpass filter
fCut = 4000; % Hz
tBin = ceil( fCut * nFFT / fs );
for n = 1:nBlocks
  F(tBin:end, n) = 0; % zero values above tBin
  F(nFFT/2+2:end, n) = conj(flipud(F(1:nFFT/2-1, n))); % make conjugate symmetric
end
plot((1:nFFT/2) * fs / nFFT, 20*log10(abs(F(1:nFFT/2, 8))));
xlabel( 'Frequency (Hz)')
ylabel( 'Magnitude (dB)')
xlim( [0 10000] );

ynew = zeros( N, 1 );
for n = 1:nBlocks
  iStart = (n-1) * hop;
  if iStart+nFFT > N, break; end
  ynew(iStart+1 : iStart+nFFT) = ynew(iStart+1 : iStart+nFFT) + w .* (real(ifft( F(:, n) )));
end

soundsc( [ y; ynew], fs );