% 2007-11-27
% FFT in matlab
[d,sr]=wavread('mpgr1_sx419.wav');
tic;sd = fft(d(1:32768));toc
% (times are from a slow CPU, for illustration)
%% elapsed_time =
%%     0.0732
tic;sd = fft(d(1:32767));toc
%% elapsed_time =
%%     0.3352
% i.e. if length is 2^N, computation is about 5x faster than 2^N-1
% 32767 is still a composite number:
factor(32767)
%ans =
%    7    31   151
% Closest *prime* to 32768 is 32749
tic;sd = fft(d(1:32749));toc
%% elapsed_time =
%%     1.2743
% i.e. 4x slower again

% Compare to direct matrix implementation
% Get the actual matrix for a DFT calculation
% (can't get much larger than 1024x1024 on my machine)
DM = dftmtx(1024);
% Calculate DFT by matrix multiply
tic; D1 = DM*d(1:1024); toc
%% Elapsed time is 0.006674 seconds.
% Compare to fft routine
tic; D2 = fft(d(1:1024)); toc
%% Elapsed time is 0.000201 seconds.
% i.e. fft is 34x faster
max(abs(D1-D2))
%% ans =
%%    1.8240e-17
% .. for the same answer (down to roundoff error)


% Plotting FFT outputs
sd = fft(d(1:32768))
subplot(211)
plot(sd)
% Not nice - plotting a complex value plots re against im
plot(abs(sd))
% OK, can now see |X(e^jw)| - but for all the way around the 
% circle, so mirror image.  Plot first half for just 
% the freqs 0 to fnyq, and plot it in dB:
plot(20*log10(abs(sd(1:16384))))
% Plot it against the actual frequencies in Hz
% The original spectrum uniformly divides 0 to 2pi (i.e. 0 to SR in Hz)
frqs = [0:16383]/length(sd)*sr;
plot(frqs, 20*log10(abs(sd(1:16384))))
% A dB range of more than 60-80 dB is pretty meaningless
axis([0 sr/2 -20 40])
grid

% STFT
% These routines (stft.m and istft.m) are not great, use them 
% at your own risk!!
% Original clean speech
soundsc(d,sr)
% Show the spectrogram using our routine
% STFT with 256 point FFT, 256 point Hamming window, 128 point hop
D = stft(d',256,256,128);
% See it
subplot(311)
imagesc(20*log10(abs(D)));
axis xy
% (specgram(d,256,sr) does more or less the same thing)
% Adjust color axis to sensible range of dB
caxis([-60 20])
colorbar % provides decode of colors

% Add white noise corruption
ln = length(d);
e = d'+randn(1,ln)*0.01;
soundsc(e,sr)
% Take STFT, pass in a rectangular window
E = stft(e,256,ones(1,256),128);
% Look at it:
subplot(312)
imagesc(20*log10(abs(E)));
axis xy
caxis([-60 20])
colorbar
% Only the very peaks of the speech energy poke through
% Make a mask to select those bits
M=(20*log10(abs(E))>-10);
subplot(313)
imagesc(M)
axis xy
% Resynthesize just those t-f cells
% Use Hamming window for reconstruction to smooth edges
% (*didn't* apply it during analysis)
fm = istft(E.*M, 256, 256, 128);
% Compare 'noise reduced' version to original:
soundsc(e,sr)
soundsc(fm,sr)
% Classic "musical noise" in resynthesis