% 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