% Testing fourier Trasform for the oscilloscope
tmin=0;
time_step=1e-6; 
tmax=10e-3;
t=tmin:time_step:tmax;

N=length(t);
amp=2;
f1=20e3;
f2=50e3;
signal=amp.*sin(2*pi*f1.*t);
%signal=amp.*sin(2*pi*f1.*t)+2.*amp.*sin(2*pi*f2.*t);
% Facotor variable is use later to compute the delta between the 2 peacks.
% Sawtooth wave
% Signal has a DC component of -28 dB if amp=2 (using mean(signal))
%signal=amp.*sawtooth(2*pi*f1*t); factor=2;
%signal=amp.*sawtooth(2*pi*f1*t)+amp.*sin(2*pi*f2.*t);
% Symmetric Saw tooth wave (triangle)
%signal=amp.*sawtooth(2*pi.*f1.*t,0.5); factor=3;
% Square Wave
%signal=amp.*square(2.*pi.*f1.*t); factor=3;

% Amplitude modulation
%fbase=1e4;
%signal=amp.*square(2.*pi.*fbase.*t).*sin(2*pi*f1*t);
%signal=amp.*sawtooth(2.*pi.*fbase.*t).*sin(2*pi*f1*t);


%figure;
%plot(t,signal);

fourier=fftshift(fft(signal))./N;

xx=real(fourier);
yy=imag(fourier);
mod=abs(fourier);
phase=angle(fourier);
mod_db=20.*log10(mod);

freq=-1./(2.*time_step):1./tmax:1./(2.*time_step);
figure;
plot(freq,mod_db)
%xlim([0 max(freq)])

delta_ind=10;
% Compute the amplitude at f1
[aux,ind1]=min(abs(freq-f1));
amp_aux=mod_db(ind1-delta_ind:ind1+delta_ind);
[max_amp1,ind_max_aux]=max(amp_aux);
ind_max1=ind1-delta_ind+ind_max_aux;
% Compute the amplitude at factor*f1
[aux,ind2]=min(abs(freq-factor*f1));
amp_aux=mod_db(ind2-delta_ind:ind2+delta_ind);
[max_amp2,ind_max_aux]=max(amp_aux);
ind_max2=ind2-delta_ind+ind_max_aux;
hold on;
plot(freq(ind_max1),max_amp1,'go')
plot(freq(ind_max2),max_amp2,'go')
hold off;
delta_peacks=max_amp1-max_amp2
delta_harmonic_num=round((freq(ind_max2)-freq(ind_max1))/f1)