## Contents

%Tutorial on Fourier Transforms:

N=1024; %number of points in the time series

OutFile1='PureToneInteger.dat';
OutFile2='PureToneNonInteger.dat';
OutFile3='Aliasing.dat';
OutFile4='PureTone_Noise.dat';
OutFile5='Rectangular.dat';
OutFile6='IrregularSampling.dat';
OutFile7='Windowed.dat';

t=0:N-1;    %time bin index array
f=0:N-1;    %frequency bin index array


## Fourier transform of a pure cosine wave (integer frequency)

    figure(1)
f0=32;          %frequency of the cosine wave
y=cos(2*pi*t*f0/N);
z=fft(y);
Y=z.*conj(z)/(N*N);   %power spectrum

fid=fopen(OutFile1,'w');
fprintf(fid,'%f\n',y);
fclose(fid);

subplot(3,1,1)
plot(t,y)
ylim([-2 2])
xlim([0 N])
xlabel('time')
ylabel('y(t)')
title(['frequency, f_0= ' num2str(f0)])
set(gca,'FontSize',16,'LineWidth',2)
subplot(3,1,2)
plot(f,Y)
ylim([0 0.5])
xlim([0 N])
ylabel('Power Spectrum, P(f)')
xlabel('Frequency (Positive & Negative)')
subplot(3,1,3)
semilogy(f(1:N/2),Y(1:N/2),'.')
ylim([1e-12 1])
xlim([0 N/2])
legend('note log scale')
ylabel('log10(P(f))')
xlabel('Positive Frequency')


## Fourier transform of a pure cosine wave (non-integer frequency)

    figure(2)
f0=100.5;
y=cos(2*pi*t*f0/N);
z=fft(y);
Y=z.*conj(z)/(N*N/4);

fid=fopen(OutFile2,'w');
fprintf(fid,'%f\n',y);
fclose(fid);

subplot(3,1,1)
plot(t,y)
title(['f_0=' num2str(100.5)])
ylim([-2 2])
xlim([0 N])
xlabel('time')
ylabel('y(t)')
set(gca,'FontSize',16,'LineWidth',2)
subplot(3,1,2)
plot(f(1:N/2),Y(1:N/2))
xlim([0 N/2])
ylabel('P(f)')
xlabel('Positive Frequency')
subplot(3,1,3)
plot(f(80:120),Y(80:120),'+','MarkerSize',6)
xlabel('Positive Frequency')
ylabel('P(f)')
xlim([80 120])


## Aliasing (what happens if your input time series is high frequency)

    figure(3)
af0=100;
f0=N-af0;    %any frequency f_0>N/2 is aliased
y=cos(2*pi*t*f0/N);
z=fft(y);
Y=z.*conj(z)/(N*N/4);

fid=fopen(OutFile3,'w');
fprintf(fid,'%f\n',y);
fclose(fid);

subplot(2,1,1)
plot(t,y)
title(['f_0=' num2str(f0) ' which is' num2str(N) '-' num2str(af0)])
ylabel('y(t)=cos(2\pi f_0t/N)','FontSize',12)
xlabel('time')
xlim([0 N])

subplot(2,1,2)
plot(f(1:N/2),Y(1:N/2))
title(['Frequency is now aliased to freq bin ' num2str(af0)])
ylabel('Power')
xlabel('Positive Frequencies')
xlim([0 N/2])


## We now add Gaussian noise to an underlying cosine

    figure(4)
f0=64;  %frequency of cosine
a=1;    %rms of Gaussian noise  (rms)

y=cos(2*pi*t*f0/N);  %amplitude of cosine
n=a*randn(1,N);  %mean=0, rms=1, Gaussian noise
yn=y+n;

fid=fopen(OutFile4,'w');    %write out the signal
fprintf(fid,'%f\n',yn);
fclose(fid);

subplot(3,1,1)
plot(t,y)
ylim([-2 2])
title(['pure tone, f_0=' num2str(f0)])
xlabel('time')
ylabel('y(t)')
xlim([0 N])

subplot(3,1,2)
plot(t,yn)
ylim([-2 2])
xlim([0 N/2])
ylabel('y(t)+n(t)')
title('added unit variance Gaussian noise: ratty signal')
subplot(3,1,3)
z=fft(yn);
Yn=z.*conj(z)/(N*N/4);
semilogy(f(1:N/2),Yn(1:N/2))
title('switched to log. Notice the noise floor')
xlim([0 N/2])
ylim([1e-3 1])
ylabel('P(f)')


## Periodic signals which are not pure sine (or cosine)

Square wave
    figure(5)
f0=16.001;
y=GenerateSignal(N,f0,0);  %my code to generate rect signal

fid=fopen(OutFile5,'w');    %write out the signal
fprintf(fid,'%f\n',y);
fclose(fid);

subplot(3,1,1)
plot(t,y);
ylim([-2 2])
ylabel('amplitude')
xlabel('time')
subplot(3,1,2)
z=fft(y);
Y=z.*conj(z)/(N*N/4)+eps;
plot(f(1:N/2),Y(1:N/2))
legend('Note the presence of odd harmonics')
ylabel('P(f)')
xlabel('f')
subplot(3,1,3)
semilogy(f(1:N/2),Y(1:N/2))
ylim([1e-3 1])
legend('Note log scale')
ylabel('P(f)')
xlabel('f')


## Irregular sampling

    figure(6)
f0=16;
cutoff=0.8;
x=cos(2*pi*t*f0/N);
a=(sign(rand(1,N)-cutoff)+1)/2;  %select only 1/10 data points
y=a.*x;

fid=fopen(OutFile6,'w');    %write out the signal
fprintf(fid,'%f\n',y);
fclose(fid);

z=fft(y);
Y=z.*conj(z)/(N*N/4)+eps;
subplot(3,1,1)
plot(t,x,'.')
xlim([1 N])
xlabel('time')
ylabel('amplitude')
title('pure tone')
subplot(3,1,2)
plot(t,y,'.')
nselect=sum(y~=0);
legend(['# points=' num2str(nselect)]);
xlabel('t')
ylabel('amplitude')
xlim([0 N])
title('sparse sampling (1 in 5)')
subplot(3,1,3)
semilogy(f(1:N/2),Y(1:N/2))
xlabel('f')
ylabel('P(f)')
xlim([0 N/2])
ylim([1e-4 1])
legend('note log scale')


## Windowing a sinusoidal signal

    figure(7)
f0=70;
F0=4;
y=WindowedSineWave(N,f0,F0);

fid=fopen(OutFile7,'w');    %write out the signal
fprintf(fid,'%f\n',y);
fclose(fid);

z=fft(y);
Y=z.*conj(z)/(N*N/4)+eps;
subplot(3,1,1)
plot(t,y)
xlabel('t')
ylabel('amplitude')
title(['f_0=' num2str(f0) ' F_0=' num2str(F0)]);
subplot(3,1,2)
plot(f(1:N/2),Y(1:N/2))
xlabel('f')
ylabel('P(f)')
subplot(3,1,3)
plot(f(1:round(2*f0)),Y(1:round(2*f0)),'+-')
xlabel('f')
ylabel('P(f)')
title('Side lobes from window function')


## GenerateSignal

%Here you will find all the routines used in this program

% <include>GenerateSignal.m</include>

% <include>WindowedSineWave</include>