%% Class 19: 
% Alex Townsend, 22nd April 2015
%        ajt@mit.edu

%% Alien sound: 
fs = 1000; 
t = 0 : 1/fs : 1 - 1/fs; 
sound = @(t) 3*cos(2*pi*20*t + 0.2) + 1*cos(2*pi*30*t - 0.3) + 2*cos(2*pi*40*t+2.4); 
samples = sound(t); 

plot(t, samples, 'linewidth',2 )
set(gca,'fontsize',16)

%% Play the sound 
% chebtune( chebfun(sound,[0,1]), 5) % cannot play this without Chebfun.

%% FFT the samples: 
modes = fft( samples ); 
mag = abs( modes ); 
plot( 0:length(mag)-1, mag, 'linewidth',2 )
set(gca,'fontsize',16) 

%% Recover phase: 
ang = angle( modes ); 
ang( [21 31 41] ) 

%% Recover amplitude: 
2*abs(modes([21 31 41]))/N

%% Another example: Non-integer frequencies
fs = 1000; 
t = 0 : 1/fs : 1 - 1/fs; 
sound = @(t) 3*cos(2*pi*29.5*t + 0.2) + 1*cos(2*pi*30.5*t - 0.3); 
samples = sound(t); 

plot(t, samples, 'linewidth',2 )
set(gca,'fontsize',16)

%% Play the sound 
% chebtune( chebfun(sound,[0,1]), 5)

%% FFT the samples, not good this time.
modes = fft( samples ); 
mag = abs( modes ); 
plot( 0:length(mag)-1, mag, 'linewidth',2 )
set(gca,'fontsize',16) 

%% Zero padding and windowing improves things
samples = samples.*hanning(length(samples))';
samples = [samples  zeros(1,20000)];
modes = fft( samples ); 
mag = abs( modes ); 
plot( 0:length(mag)-1, mag, 'linewidth',2 )
set(gca,'fontsize',16)

%% Recover frequency, amplitude, and phase.
[mx, idx] = max(abs(modes));
freq = idx / 21000 * 1000
amp = mx/1000*4 
phase = angle( modes(idx) ) 

%% Why N samples, but only N/2 modes? Nyquist sampling rate.
t = linspace(0,1,1000); 
plot(t, cos(2*pi*t), 'linewidth',2 )
set(gca,'fontsize',16), hold on 

%% Aliasing: 
t1 = linspace(0,1,10); 
plot(t1, cos(2*10*pi*t1),'k.','markersize',30)

%% Aliasing:
plot(t, cos(2*10*pi*t), 'linewidth',2)

%% Aliasing:
t = linspace(0,1,1000); 
plot(t, cos(2*2*pi*t), 'linewidth',2 )
set(gca,'fontsize',16), hold on 

%% Aliasing:
t1 = linspace(0,1,10); 
plot(t1, cos(2*11*pi*t1),'k.','markersize',30)
%%
plot(t, cos(2*11*pi*t), 'linewidth',2)

%% Police car sound. Cannot play this with Chebfun. 
%x = chebfun('x',[0,1]); 
%f = 9 + 6*sin(92*x); 
%g = 7 + 10*sin(8*x);
%chebtune([f, g], 5);