% M14-fft.diary
% FFT in matlab
% 2007-11-27

[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
% (ram usage gets severe quite fast for > 4096!))
DM = dftmtx(4096);
% Calculate DFT by matrix multiply
tic; D1 = DM*d(1:4096); toc
%% Elapsed time is 0.073838 seconds.
% Compare to fft routine
tic; D2 = fft(d(1:4096)); toc
%% Elapsed time is 0.000217 seconds.
% i.e. fft is 340x faster
max(abs(D1-D2))
%% ans =
%%    4.0131e-15
% .. 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

% Modulated sine as sum of sines of near frequencies
% time base (2^N points long for N=14)
sr = 8192;
t = [0:16383]/sr;
% three cosines that exactly fit, with near frequencies
c1 = cos(2*pi*1000*t);
c2 = cos(2*pi*1002*t);
c3 = cos(2*pi*998*t);
% each is just a steady tone
soundsc(c1,sr)
soundsc(c2,sr)
soundsc(c3,sr)
% .. but they sum to a modulated carrier
y = c1 - 0.5*(c2+c3);
soundsc(y,sr)
plot(y)
% The DFT of y is just the three points
stem(abs(fft(y)))
% .. but that's not the most natural way to describe what we hear.

% Interactive changing of window length
% funcslider is a very slick wrapper to re-evaluate and display a function 
% with arguments changed by sliders
[d,sr]=wavread('mpgr1_sx419.wav');
funcslider(@(x) specgram(d,1024,sr,round(1024*x),round(896*x)), 'imagesc(flipud(max(-60,20*log10(abs(z1)))))')

% 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, use 1/4 window hop so cos * cos window still summs OK
E = stft(e, 256, 256, 64);
% 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))>-15);
subplot(313)
imagesc(M)
axis xy
% Resynthesize just those t-f cells
% (We apply Hamming window on both analysis and synthesis
% but it's OK because we have 75% overlap)
fm = istft(E.*M, 256, 256, 64);
% Compare 'noise reduced' version to original:
soundsc(e,sr)
soundsc(fm,sr)
% Classic "musical noise" in resynthesis