% aud3.m
% See FFTs of music
% 2011/05/05  Written. (B.Baas)
% 2011/05/09  Cleaned up, "mode" variable at top added. (B.Baas)
% 2011/05/20  Auto-tune portion is getting closer. (B.Baas)

clear;
%===== set these variables
% 0: speedup auto-tune fail
% 1: see spectrogram
% 2: real time fft. Set "axis" command parameters to freeze movie view.
%    Resize figure 2 window to fill screen to see maximum detail.
mode = 0;
length_clip = 9;   % in seconds, instead of full clip

%===== read sound file
[all,fs] = wavread('gaga.wav');
%[a,fs] = m4aread('SomebodysCrying.m4a');
%a = mp3read('MustBeLove.mp3');
a = all(1:fs*length_clip);
len_a = length(a);
fprintf('Sample rate fs = %f samples/sec\n', fs);
fprintf('Number of samples = %i\n', length(a));
fprintf('Seconds of music = %f\n', len_a/fs);
        
%===== Speedup with ifft/fft
% Not yet working. I think this can be fixed by keeping the full
% frequency-domain representation and making sure the spectrum is
% symmetric about fs/2 before doing the ifft().
if(mode == 0)
  x1 = fft(a);
  x2 = x1(1:(length(x1)/2));  len_x2 = length(x2);
  shiftup = 100;
  %x3 = [zeros(shiftup,1) ; x2(1:(len_x2-shiftup))];  % col vector
  x3 = [zeros(1,shiftup) x2(1:(len_x2-shiftup)).'];
  x4 = [x3 x3(end:-1:1)];
  a_higher = ifft(x4);
  %     abs(ifft([zeros(shiftup,1) ; x1(1:(len_a-shiftup))]));
  keyboard
  sound(real(a_higher), fs);
  junk = input('  <Press return>');
  clear playsnd;
end

%===== Specgram
if(mode == 1)
  figure(3);
  specgram(a);
end

%===== Real-time FFT displays
len_fft = 20000/1;  % Set this variable to determine 1) how much resolution
                    % each displayed FFT has, and 2) how often the display
                    % is updated. Remember there are fs samples/sec which
                    % is typically 44,100.
                    % Choose len_fft so that (len_a/len_fft) is an integer.

if(mode == 2)
  time_fft = len_fft/fs;
  figure(2); %clf;
  for k=1:(floor(len_a/len_fft))
    a2 = a((k-1)*len_fft+1 : k*len_fft);
    ff1 = fft(a2);               % take fft of a section
    ff2 = ff1(1:(len_fft/2));    % extract left half only
    plot(abs(ff2));
    %grid on;
    axis([0 1000 0 400]);        % <== Comment and/or adjust as needed!
    sound(a2,fs);
    pause(time_fft*0.6);
  end
end