%Fast Fourier transform (FFT) 
%MATLAB tutorial series (Part 2.1)

%% Example 2.1: pure tone

% Sampling frequency
Fs=44100;
% Sampling period
Ts=1/Fs;
% Signal length (in second)
N_sec=5;
% Signal length (in sample)
N=N_sec*Fs;
% Maximum time
Tmax=(N-1)*Ts;
% Time vector
t=0:Ts:Tmax;
% pure tone freq (e.g. 300, 600, 2000, ...)jh
F0=600;
% Signal
x=sin(2*pi*F0.*t);
% Play the sound
sound(x,Fs)
% Plot the sound (show from zero to 60 msec)
figure(1)
plot(t,x,'LineWidth',2.5)
xlim([0 0.06])
% take FFT (without shift)
X1=fft(x);
% plot result
figure(2)
plot(abs(X1),'LineWidth',2.5);
% take FFT (with shift)
X2=fftshift(fft(x));
% Frequency range
F=-Fs/2:Fs/(N-1):Fs/2;
figure(3)
plot(F,abs(X2)/N,'LineWidth',2.5);
xlabel('Frequency (Hz)')
title('Double sided magnitude response')


%% Example 2.2: single side spectrum   

% % Sampling frequency
% Fs=44100;
% % Sampling period
% Ts=1/Fs;
% % Signal length (in second)
% N_sec=5;
% % Signal length (in sample)
% N=N_sec*Fs;
% % Maximum time
% Tmax=(N-1)*Ts;
% % Time vector
% t=0:Ts:Tmax;
% % pure F0, F1 and F3
% F0=600;
% F1=1300;
% F2=2000;
% % Signal
% x=sin(2*pi*F0.*t)+ 0.5*sin(2*pi*F1.*t)+0.2*sin(2*pi*F2.*t);
% % Play the sound
% sound(x,Fs)
% % Plot the sound (show from zero to 60 msec)
% figure(1)
% plot(t,x,'LineWidth',2.5)
% xlim([0 0.06])
% % Compute fft
% X=fft(x);
% % Take abs and scale it
% X2=abs(X/N);
% % Pick the first half
% X1=X2(1:N/2+1);
% % Multiply by 2 (except the DC part), to compenseate 
% % the removed side from the spectrum.
% X1(2:end-1) = 2*X1(2:end-1);
% % Frequency range
% F = Fs*(0:(N/2))/N;
% % Plot single-sided spectrum
% figure(2)
% plot(F,X1,'LineWidth',2.5) 
% title('Single-Sided Amplitude Spectrum')
% xlabel('f (Hz)');


%% Example 2.3: siren   

% % Sampling frequency
% Fs=44100;
% % Sampling period
% Ts=1/Fs;
% % Signal length (in second)
% N_sec=15;
% % Signal length (in sample)
% N=N_sec*Fs;
% % Maximum time
% Tmax=(N-1)*Ts;
% % Time vector
% t=0:Ts:Tmax;
% % pure F0, F1 and F3
% F0=1300;
% F1=200;
% F2=1400;
% B0=100;
% B1=100;
% B2=500;
% % Signal
% x=sin(2*pi*F0.*t+B0*pi*t.^2)+sin(2*pi*F1.*t+B1*pi*t.^2)...
% +sin(2*pi*F2.*t+B2*pi*t.^2);
% % Play the sound
% sound(x,Fs)
% % Plot the sound (show from zero to 600 msec)
% figure(1)
% plot(t,x,'LineWidth',2.5)
% xlim([0 0.06])
% % Compute fft
% X=fft(x);
% % Take abs and scale it
% X2=abs(X/N);
% % Pick the first half
% X1=X2(1:N/2+1);
% % Multiply by 2 (except the DC part), to compenseate 
% % the removed side from the spectrum.
% X1(2:end-1) = 2*X1(2:end-1);
% % Frequency range
% F = Fs*(0:(N/2))/N;
% % Plot single-sided spectrum
% figure(2)
% plot(F,X1,'LineWidth',2.5) 
% title('Single-Sided Amplitude Spectrum')
% xlabel('f (Hz)');

%% Example 2.4: voice 

% % read audio file .wav
% [x,Fs]=audioread('adult_female_speech.wav');
% % play the sound
% sound(x,Fs)
% % Sampling period
% Ts=1/Fs;
% % Length of signal
% N=length(x);
% % Maximum time
% Tmax=(N-1)*Ts;
% % Time vector
% t=0:Ts:Tmax;
% % Plot the sound )
% figure(1)
% plot(t,x,'LineWidth',2.5)
% % Compute fft
% X=fft(x);
% % Take abs and scale it
% X2=abs(X/N);
% % Pick the first half
% X1=X2(1:N/2+1);
% % Multiply by 2 (except the DC part), to compenseate 
% % the removed side from the spectrum.
% X1(2:end-1) = 2*X1(2:end-1);
% % Frequency range
% F = Fs*(0:(N/2))/N;
% % Plot single-sided spectrum
% figure(2)
% plot(F,X1,'LineWidth',2.5) 
% title('Single-Sided Amplitude Spectrum')
% xlabel('f (Hz)');

%% Example 2.5: voice 

% % read audio file .wav
% [x,Fs]=audioread('Ghost.wav');
% % play the sound
% sound(x,Fs)
% % Sampling period
% Ts=1/Fs;
% % Length of signal
% N=length(x);
% % Maximum time
% Tmax=(N-1)*Ts;
% % Time vector
% t=0:Ts:Tmax;
% % Plot the sound 
% figure(1)
% plot(t,x,'LineWidth',2.5)
% % Compute fft
% X=fft(x);
% % Take abs and scale it
% X2=abs(X/N);
% % Pick the first half
% X1=X2(1:N/2+1);
% % Multiply by 2 (except the DC part), to compenseate 
% % the removed side from the spectrum.
% X1(2:end-1) = 2*X1(2:end-1);
% % Frequency range
% F = Fs*(0:(N/2))/N;
% % Plot single-sided spectrum
% figure(2)
% plot(F,X1,'LineWidth',2.5) 
% title('Single-Sided Amplitude Spectrum')
% xlabel('f (Hz)');