%Lecture_041612.m
%
%1-D Fourier Analysis
%Based on Lesson 11

%% Define a sinusoid 

dur = 1;
dt = .01; 
nCycles = 10; 
phase = 0;
amp = 1;

t = 0:dt:(dur-dt);

y = amp*cos(nCycles*2*pi*t-phase*pi/360);



figure(1)
clf
subplot(1,2,1)
plot(t,y);
xlabel('time (s)');

%% take its FFT

fftY = fft(y);

%% convert the output into sensible numbers using complex2real

Y = complex2real(fft(y),t);

%plot the amplitudes
subplot(1,2,2)
stem(Y.freq,Y.amp);
set(gca,'XLim',[1,Y.freq(end)]);

%% fft of the sum of two sinusoids

y = cos(2*pi*4*t) + sin(2*pi*7*t) + 10*t;

figure(2)
clf
subplot(1,2,1)
plot(t,y);
xlabel('time (s)');

Y = complex2real(fft(y),t);
subplot(1,2,2)
stem(Y.freq,Y.amp);
set(gca,'XLim',[1,Y.freq(end)]);

%% reconstructing our time series by adding up sinuoids of its fft 

yhat = Y.dc*ones(size(t));
for i=1:length(Y.freq)
    yhat = yhat + Y.amp(i)*cos(2*pi*t*Y.freq(i)-Y.ph(i)*pi/180);
end

figure(3)
plot(t,yhat,'b-');
hold on
plot(t,y,'r.');

%% the fft of some friends

y = t<.5;  %step function

%noise
%y = randn(size(t));

figure(2)
clf
subplot(1,2,1)
plot(t,y);
xlabel('time (s)');

Y = complex2real(fft(y),t);
subplot(1,2,2)
stem(Y.freq,Y.amp);
set(gca,'XLim',[1,Y.freq(end)]);

%%

%mess with the phases

Y.ph = 360*rand(size(Y.ph));


yhat = ifft(real2complex(Y));

figure(3)
clf
plot(t,yhat)

%% the Gabor!

nCycles = 15;
sig = .05;
amp = 1;
phase  =0;
center = .5;

% A Gabor is the product of a Gaussian and a sinusoid
% The Gaussian is the 'envelope' and the sinusoid is the 'carrier'

y = amp*sin(nCycles*2*pi*t-phase*pi/180) .* ...
    exp(-(t-center).^2/(2*sig^2));

figure(4)
clf
subplot(1,2,1)
plot(t,y);
xlabel('time (s)');


%the FFT of a Gabor is a Gaussian centered at the carrier frequency.
%the width of the Guassian is inversely proportional to the width of the
%envelope. 
Y = complex2real(fft(y),t);
subplot(1,2,2)
stem(Y.freq,Y.amp);
set(gca,'XLim',[1,Y.freq(end)]);

%% Filtering sound
clear all

%all files:
%chirp.mat     handel.mat    splat.mat     
%gong.mat      laughter.mat  train.mat     


load train

%sound(y,Fs)
maxT = length(y)/Fs;

t = linspace(0,maxT,length(y))';
figure(1)
clf
plot(t,y);
xlabel('Time (s)');

%Take and plot the FFT
Y = complex2real(fft(y),t);

figure(2)
clf
plot(Y.freq,Y.amp,'-');

xlabel('Frequency (Hz)');

%hard edged filter keeping amplitudes between these frequencies:
keepFreqs = [1000,1400];

Y.amp(Y.freq<keepFreqs(1) | Y.freq>keepFreqs(2)) = 0;

hold on
plot(Y.freq,Y.amp,'r-');

yhat = ifft(real2complex(Y));
sound(y,Fs)
sound(yhat,Fs)

%% create a spectrograph

dt = .01;
sig = .025;

tPlot =0:dt:maxT;

S = zeros(length(Y.freq),length(tPlot));

for i=1:length(tPlot)
    %mulitply the time series by a Gaussian centered at t(i)
    
    Gauss = exp(-(t-tPlot(i)).^2/(2*sig^2));
    Y = complex2real(fft(y.*Gauss),t);
    S(:,i) = Y.amp;
end
%
figure(5)
clf
imagesc(t,Y.freq,S)
xlabel('Time (s)');
ylabel('Frequency (Hz)');
set(gca,'YDir','normal');
set(gca,'YLim',[0,1500]);

cmap=hot(64);
cmap = cmap.^.5;
colormap(cmap);
colorbar

    












