%--------------------------------------------------------------------------
%                    INTRODUCTION TO FFT ANALYSIS with MATLAB PHYS 460/660
%--------------------------------------------------------------------------


%                          ON-LINE HELP 

help fft

% Figure C.2 of the textbook:

dt=0.1;
t=0:dt:12.7;
signal_1=sin(2*pi*0.2362*t);
figure(1),plot(t,signal_1,'-o');

N=128;
ft_signal_1=fft(signal_1);
F=[0:N/2-1]/(N*dt);
figure(2),plot(F(1:64),imag(ft_signal_1(1:64)),'o')


% Figure C.3 of the textbook:

dt=0.1;
t=0:dt:12.7;
signal_2=sin(2*pi*0.2362*t+pi/4);

N=128;
ft_signal_2=fft(signal_2,N);

F=[0:N/2-1]/(N*dt);
figure(3),
subplot(3,1,1)
plot(t,signal_2,'-o'),title('Original Signal')
subplot(3,1,2)
plot(F(1:32),real(ft_signal_2(1:32)),'o'),title('Real part of FFT')
subplot(3,1,3) 
plot(F(1:32),imag(ft_signal_2(1:32)),'o'),title('Imaginary part of FFT')


% Figure C.4 of the textbook:

dt=0.1;
t=0:dt:12.7;
signal_3=sin(2*pi*0.5*t)+3*sin(2*pi*1.5*t+pi/8);

N=128;
ft_signal_3=fft(signal_3,N);

F=[0:N/2-1]/(N*dt);
figure(4),
subplot(3,1,1)
plot(t,signal_3,'-o'),title('Original Signal')
subplot(3,1,2)
plot(F(1:32),real(ft_signal_3(1:32)),'o'),title('Real part of FFT')
subplot(3,1,3) 
plot(F(1:32),imag(ft_signal_3(1:32)),'o'),title('Imaginary part of FFT')



%Figure C.5 of the textbook:

dt=0.1;
t=0:dt:6.3;
signal_4=sin(2*pi*4*t);

N=64;
ft_signal_4=fft(signal_4);


F=[0:N/2-1]/(N*dt);
figure(6),
subplot(3,1,1)
plot(t,signal_4,'-o')
subplot(3,1,2)
plot(F(1:32),real(ft_signal_4(1:32)),'o')
subplot(3,1,3) 
plot(F(1:32),imag(ft_signal_4(1:32)),'o')


%Figure C.6 of the textbook:

dt=0.2;
t=0:dt:6.3;
signal_5=sin(2*pi*4*t);

N=32;
ft_signal_5=fft(signal_5);


t1=0:0.05:6.3;
signal_4_full=sin(2*pi*4*t1);

F=[0:N/2-1]/(N*dt);
figure(7),
subplot(3,1,1)
plot(t,signal_5,'-o',t1,signal_4_full)
subplot(3,1,2)
plot(F(1:16),real(ft_signal_5(1:16)),'o')
subplot(3,1,3) 
plot(F(1:16),imag(ft_signal_5(1:16)),'o')

% COMMENT: Note that "true" frequency of the signal is f_true=4 Hz, 
% while f_Nyquist=1/(2*dt)=2.5 Hz; So, reflected frequency is f_true-f_Nyquist=1.5 Hz 
% below f_Nyquist!


% Figure C.7 of the textbook:

power_signal_3=abs(ft_signal_3).^2/128;

N=128;
F=[0:N/2-1]/(N*dt);
figure(8),
plot(F(1:64),power_signal_3(1:64),'-o')




%                          POWER SPECTRUM ANALYSIS WITH MATLAB FFT

% A common use of Fourier transforms is to find the frequency components of a signal buried in a 
% noisy time domain signal. Consider data sampled at 1000 Hz. Form a signal containing 50 Hz and 120 Hz: 
t = 0:0.001:0.6;
x = sin(2*pi*50*t)+sin(2*pi*120*t);

% Corrupt the signal with it some zero-mean random noise:
y = x + 2*randn(size(t));

help randn

% Plot the final signal in time domain:
figure(3),plot(y(1:600))

% It is difficult to identify the frequency components from looking at the original signal. 
% Converting to the frequency domain, the discrete Fourier transform of the noisy signal y 
% is found by taking the 512-point fast Fourier transform (FFT):
Y = fft(y,512);

% The power spectral density, a measurement of the energy at various frequencies, is:
Pyy = abs(Y).^2 / 512;

% The first 256 points (the other 256 points are symmetric) can be graphed on a 
% frequency axis where we input length of the physical time interval 0.001s:
f = 1000*(0:255)/512;
figure(4),plot(f,Pyy(1:256)), xlabel('Frequency (Hz)'), ylabel('Power Spectrum')



%                POWER SPECTRUM ANALYSIS WITH MATLAB FFT on EXTERNAL FILE

% Load external file (containing, e.g., time vs. position of particle 27 in
% our nonlinear chain of 100 atoms) into a matrix nonlinear_xy.dat

load signal_xy.dat;

% Copy first column of this matrix into the time vector
t=signal_xy(:,1);

% Copy second column of this matrix into the fpu_signal vector  
fpu_signal=signal_xy(:,2);

figure(5),plot(t,fpu_signal)


Y = fft(fpu_signal,512);

% The power spectral density, a measurement of the energy at various frequencies, is:
Pyy = abs(Y).^2 / 512;

% The first 256 points (the other 256 points are symmetric) can be graphed on a 
% frequency axis where we input length of the physical time interval 0.001s:
f = 1000*(0:255)/512;
figure(6),plot(f,Pyy(1:256)), xlabel('Frequency (Hz)'), ylabel('Power Spectrum')




%                FFT and ARTIFACTS of WINDOWS CONVOLUTED WITH A FUNCTION


%Vary the number of repetitions of the fundamental period
n=0:30;
signal_3=cos(2*pi*n/10); %3 periods
signal_6=[signal_3 signal_3]; %6 periods
signal_9=[signal_3 signal_3 signal_3]; %9 periods

N=512;

ampl_s1=abs(fft(signal_3,N));
ampl_s2=abs(fft(signal_6,N));
ampl_s3=abs(fft(signal_9,N));

%A real number signal should have a transform amplitude that is symmetrical for
% positive and negative frequencies

help fftshift

ampl_s1=fftshift(ampl_s1);
ampl_s2=fftshift(ampl_s2);
ampl_s3=fftshift(ampl_s3);

% Spectrum now goes from -f_s/2 to f_s/2 where f_s is the sampling
% frequency and f_s/2 is the Nyqist frequency

F=[-N/2:N/2-1]/N;

figure(2)
subplot(3,1,1)
plot(F,ampl_s1),title('3 periods'),
subplot(3,1,2)
plot(F,ampl_s2),title('6 periods'),
subplot(3,1,3)
plot(F,ampl_s3),title('9 periods'), xlabel('frequency/f_s')





%                                MORE ABOUT ALIASING IN DISCRETE FOURIER TRANSFORM 


% MATLAB example of a continuous-time signal being sampled at various
% frequencies, illustrating the problem of aliasing caused by sampling
% at too low of a frequency.

% Sampling periods and sampling frequencies to be used.

% ws1 = 125.6637, ws2 = 62.8319, ws3 = 31.4159

%

fig_size = [232 84 774 624];

Ts1 = 0.05; Ts2 = 0.1; Ts3 = 0.2;

ws1 = 2*pi/Ts1; ws2 = 2*pi/Ts2; ws3 = 2*pi/Ts3;

%

% Frequencies for the continous-time signal and the time vector.

%

w1 = 7; w2 = 23;

t = [0:0.005:2];

%

% Original continuous-time signal is the sum of two cosines, with

% frequencies of 7 r/s and 23 r/s, respectively.

%

x = cos(w1*t) + cos(w2*t);

figure(5),clf,plot(t,x),grid,xlabel('Time (s)'),ylabel('Amplitude'),...

title('Continuous-Time Signal; x(t) = cos(7t) + cos(23t)'),...

set(gcf,'Position',fig_size)

%

% Sampling the continuous-time signal with a sampling period Ts = 0.05 s.
% The sampled signal is exactly equal to the continuous-time signal at the
% sample times, and the samples accurately model the original signal in
% the following respect: if you look at the samples by themselves and
% wanted to guess what the continuous-time signal looks like, you would be
% able to get pretty close. Note that ws1 is approximately 5.5*w2.

%

t1 = [0:Ts1:2];

xs1 = cos(w1*t1) + cos(w2*t1);

figure(6),clf,stem(t1,xs1);grid,hold on,plot(t,x,'r:'),hold off,...

xlabel('Time (s)'),ylabel('Amplitude'),...

title('Sampled Version of x(t) with T_s = 0.05 s'),...

set(gcf,'Position',fig_size)

%

% Sampling the continuous-time signal with a sampling period Ts = 0.1 s.
% The sampled signal is exactly equal to the continuous-time signal at the
% sample times. The samples are a less accurate representation of the
% original signal than with the smaller Ts (higher sampling frequency ws).
% Note that ws2 is approximately 2.7*w2.

%

t2 = [0:Ts2:2];

xs2 = cos(w1*t2) + cos(w2*t2);

figure(7),clf,stem(t2,xs2);grid,hold on,plot(t,x,'r:'),hold off,...

xlabel('Time (s)'),ylabel('Amplitude'),...

title('Sampled Version of x(t) with T_s = 0.1 s'),...

set(gcf,'Position',fig_size)

%

% Sampling the continuous-time signal with a sampling period Ts = 0.2 s.
% The sampled signal is exactly equal to the continuous-time signal at the
% sample times. The samples now are not a good representation of the
% original signal at all. Note that ws3 is approximately 1.37*w2.

%

t3 = [0:Ts3:2];

xs3 = cos(w1*t3) + cos(w2*t3);

figure(8),clf,stem(t3,xs3);grid,hold on,plot(t,x,'r:'),hold off,...

xlabel('Time (s)'),ylabel('Amplitude'),...

title('Sampled Version of x(t) with T_s = 0.2 s'),...

set(gcf,'Position',fig_size)

%

% Since ws3 < 2*w2, the Nyquist Sampling Theorem is violated, and x(t)
% could not be recovered from the samples obtained with Ts3 using an ideal
% low-pass filter. Aliasing has occurred. The samples of the original
% x(t) using a sampling period Ts3 have exactly the same values that the
% signal x1(t) = cos(w1*t) + cos((w2-ws3)*t) would have when sampled with
% a sampling period Ts3. w2 - w3 = -8.4159 r/s.

%

w2s3 = w2 - ws3;

x1 = cos(w1*t) + cos(w2s3*t);

%

figure(9),clf,stem(t3,xs3);grid,hold on,plot(t,x,'k:',t,x1,'r:'),...

hold off,xlabel('Time (s)'),ylabel('Amplitude'),...

title('Sampling x(t) and x_1(t) with T_s = 0.2 s'),...

set(gcf,'Position',fig_size),...

text(1.13,1.2,'x(t)'),text(0.1,1.6,'x_1(t)')