%----------------- Fourier Analysis -----------------
%
% Technical University of Sofia - MUSSPI lab course
% Created by Viktor Stoev
% October 2016
%----------------------------------------------------

clear all;
close all;

plength = -pi:0.01:pi;
length = plength(1,1:end-1);
s = cat(1,-1*ones(size(length,2)/2,1),ones(size(length,2)/2,1));

figure(1);
plot(length,s,'LineWidth',2);
xlim([-3.14 3.14]);
ylim([-1.1 1.1]);
title('Square shaped signal with period ot 2pi');
ylabel('Amplitude');
xlabel('Period of 2pi');

%% Fourier Series - approximate signal "s" with a series of sinusoids

while 1
    number = input('Choose the order of approximation (max = 20): ');
    if (number > 20)
        display('Enter a smaller number (max = 20)');
        number = input('Enter a correct number: ');
    else
        break
    end
end

% approximation formulae: f(x) = 4/pi*(sin(x)/1 + sin(3x)/3 + sin(5x)/5 + ...)
vec = 1:2:2*number;
sin_mat = zeros(number,size(length,2));  % allocate memory for the variable, speeds up the process

for ii = 1:number
    sin_mat(ii,:) = sin(vec(ii)*length)./vec(ii);
    sin_sum = (4/pi)*sum(sin_mat,1);
    figure(2);
    plot(length,s,'LineWidth',2);
    hold on;
    plot(length,sin_sum,'r','LineWidth',2);
    hold off;
    xlim([-3.14 3.14]);
    ylim([-1.5 1.5]);
    title('Square shaped signal and fourier series approximation');
    ylabel('Amplitude');
    xlabel('Period of 2pi');
    pause
end

%% Fourier Transform

clear all;
close all;

fs = 100; % Hz
T = 1/fs; % s 
t = 0:0.0001:5*T; % time vector

f = fs*((0:(size(t,2))/2)/(T*(size(t,2)))); % frequency vector

s1 = sin(2*pi*fs*t);
s2 = sin(2*pi*1.2*fs*t);
s3 = sin(2*pi*4*fs*t);
s4 = sin(2*pi*16*fs*t);

 %noise = input('Add gaussian noise ? [Y/N]: ','s');
 %if (strcmp(noise,'Y'))
     %amp = input('Choose noise magnitude ((play with different values and observe the results)) [0.0:2]: ');
     %s = s1+s2+s3+s4+amp*randn(1,size(t,2));
 %else
    %s = s1+s2+s3+s4;
 %end
%FT = abs(fft(s));
 %FT = real(fft(s));
% FT2 = imag(fft(s));

figure(3);
subplot(2,1,1);
plot(t,s);
title('Sum of four sinusoids in time domain');
ylabel('Amplitude');
xlabel('Time,[s]');
subplot(2,1,2);
plot(f,FT(1:end/2+1));
title('Fourier transform of the sum of four sinusoids in time domain');
ylabel('Magnitude');
xlabel('Frequency, [Hz]');

%% Zero-padding

z = 4; % zero-padding length (times signal length), play with this parameter as well
f_zero = fs*((0:(size(t,2)*z)/2)/(T*size(t,2)*z/2));

FT_zero = abs(fft(s,z*size(s,2)));

figure(4);
subplot(2,1,1);
plot(f,FT(1:end/2+1));
title('Fourier transform of the sum of four sinusoids in time domain');
ylabel('Magnitude');
xlabel('Frequency, [Hz]');
subplot(2,1,2);
plot(f_zero(1,1:size(FT_zero(1:end/2+1),2)),FT_zero(1:end/2+1));
title('Zero Padding result');
ylabel('Magnitude');
xlabel('Frequency, [Hz]');


%% Leakage

%wintool % create windows

t1 = 110;
t2 = 390;
zero_vec = zeros(1,size(s,2));
zero_vec(1,t1:t2) = 1;
s_t = s.*zero_vec;
FT_t = abs(fft(s_t));

figure(5);
subplot(3,1,1)
plot(t,s);
title('Sum of four sinusoids in time domain');
ylabel('Amplitude');
xlabel('Time,[s]');
hold on;
plot(t,s_t,'r');
title('Truncated signal in time domain');
ylabel('Amplitude');
xlabel('Time,[s]');
subplot(3,1,2)
plot(f,FT(1:end/2+1));
title('Fourier transform of the sum of four sinusoids in time domain');
ylabel('Magnitude');
xlabel('Frequency, [Hz]');
subplot(3,1,3)
plot(f,FT_t(1:end/2+1));
title('Truncated FT');
ylabel('Magnitude');
xlabel('Frequency, [Hz]');

wintool % create windows with a tool, try out different windows and
% observe the results
w = hann(size(s,2)).';
FT_leak = abs(fft(w.*s_t));

figure(6);
subplot(2,1,1)
plot(f,FT_t(1:end/2+1));
title('Fourier transform without windows');
ylabel('Magnitude');
xlabel('Frequency, [Hz]');
subplot(2,1,2)
plot(f,FT_leak(1:end/2+1));
title('Fourier transform with applied window');
ylabel('Magnitude');
xlabel('Frequency, [Hz]');
