%Lab 6

%Fourier transform lab

%% 1 Plotting sinusoids
figure %Open figure
fs = 1000; %Sampling frequency: 1000 Hz
dur = 1; %Over 1 second
freq = 2; %Signal frequency 
t = 0:1/fs:dur; %Timebase: Sampling over all time in steps of sampling rate
sinW = sin(2*pi*freq*t); %Code for a sine wave of a given freq over time
cosW = cos(2*pi*freq*t); %Cosine sinusoid
%Changing phase: Adding or subtracting a number from 0 to 2 pi
%cosW = cos(2*pi*freq*t); %Cosine sinusoid-pi/2
plot(t,sinW,'color','k','linewidth',4)
hold on
plot(t,cosW,'color','b','linewidth',4)

%% 2 Building up a complex signal from sine waves and showing power
%Then representing it in frequency space
figure
dur = 4; %New duration
t = 0:1/fs:dur; %New time base
freq = 10; %New frequency
signal = sin(2*pi*freq*t); %Another sine wave
subplot(3,1,1)
plot(t,signal,'color','k','linewidth',2)

freq2 = 7; %2nd frequency
signal2 = 0.7.*sin(2*pi*freq2*t); %Another sine wave, at 70% amplitude, ~50% power 
subplot(3,1,2)
plot(t,signal2,'color','k','linewidth',2)

jointSignal = signal+signal2; %Linear addition
subplot(3,1,3)
plot(t,jointSignal,'color','k','linewidth',2)
shg

%Make the signal more realistic: Add noise
n = length(t); %How much noise?
noise = randn(1,n); %Just random noise from a normal distribution
signalAndNoise = jointSignal + noise; %Linear addition
figure
subplot(3,1,1)
plot(t,jointSignal,'color','k','linewidth',2)
subplot(3,1,2)
plot(t,signalAndNoise,'color','k','linewidth',2)

%Recovering the original signal frequencies even in relatively
%simple signals - 2 sine waves plus a little bit of noise is
%hopeless from visual inspection. 
%But no problem for the Fourier transform:

nyquist = fs/2; %Half the sampling frequency
fSpaceSignal = fft(signalAndNoise)/(length(t)/2); %Scaled fft
fBase = linspace(0,nyquist,floor(length(signalAndNoise)/2)+1);%Matlab indexing off by 1, so we need to correct for this
powerPlot = subplot(3,1,3); %Open the 3rd subplot
%As you recall, the fft gives you the fft of the entire signal, but
%half of that (all beyond Nyquist) will be redundant
%So we plot only half that
halfTheSignal = fSpaceSignal(1:length(fBase));
complexConjugate = conj(halfTheSignal); %We want to plot power!
pow = halfTheSignal.*complexConjugate; %This is power
h = plot(fBase,pow,'color','k','linewidth',2);
xlim([0 20]); ylim([0 1])
powerPlot.XTick = 0:20; 
xlabel('Frequency in Hz')
ylabel('Power')

%% 3 Recovering the Phase information of the signal
%In addition to amplitude, the fft returns information about the phase
%of the signals. 

%Defining a new simple signal of variable phase but constant frequency
figure 
for theta = -2*pi:0.05:2*pi %Doing a couple of loops
dur = 1; %New duration
t = 0:1/fs:dur; %New time base
freq = 20; %New frequency
signal = sin(2*pi*freq*t+theta); %Another sine wave
fftSignal = fft(signal);
fBase = linspace(0,nyquist,floor(length(signal)/2)+1);%Matlab indexing off by 1, so we need to correct for this

subplot(3,1,1)
plot(t,signal)
title('Signal in the time domain')
xlabel('time in seconds')
ylabel('Amplitude')
subplot(3,1,2)
pow = fftSignal(1:length(fBase)).*conj(fftSignal(1:length(fBase)));
plot(fBase,pow)
xlim([0 50])
title('Signal power in the frequency domain')
xlabel('Frequency')
ylabel('Power')
subplot(3,1,3)
[val,ind] = max(pow); %Recover the right index (the one with max power)
phaseAngle = angle(fftSignal(ind)); %The range space of this function is -pi to pi - it wraps around
plot(phaseAngle,1,'o','markerfacecolor','k')
xlim([-pi pi])
ylim([0 1.1])
title(['Phase angle of the signal where power is maximal: ', num2str(phaseAngle,'%1.2f')])
xlabel('Phase angle')
pause(0.1)
end

%To recover the phase information of a given frequency, we can use its
%index, but we need to account for MATLAB's indexing by 1. 

%% 4 Aliasing and sampling
%The important thing to keep in mind when thinking about
%sampling and aliasing is that there are 3 relevant frequencies
%The sampling rate of reality, that of the signal and that of your sampling

realSampling = 1000; %Let's say reality updates 1000 times a second
t = 0:1/realSampling:1; %Our new time base - one second
signalFrequency = 10; %10 Hz
samplingRate = [3 4 8 10 11 12 20 41 200]; 
signal = sin(2*pi*signalFrequency*t);
figure
for ii = 1:length(samplingRate)
subplot(3,3,ii)
plot(t,signal)
hold on
sampledSignal = round(linspace(1,length(t),samplingRate(ii)));
q = plot(t(sampledSignal), signal(sampledSignal),'color','r');
q.Marker = 'o'; %Little circles
q.MarkerFaceColor = 'k'; %Black fill
q.MarkerEdgeColor = 'y'; %Yellow edges
hT = title(['fs = ', num2str(samplingRate(ii))])
hT.FontSize = 26;
%Samples = Measures (of the signal)
end

%% 5 Complex sine waves
figure
trace = exp(i*2*pi*freq2*t) %This will be our trace. This is the term that went into the Fourier transform for one frequency
comet3(t,real(trace),imag(trace)); %Let's do a comet
xlabel('t'); ylabel('real'); zlabel('imaginary')

%Complex exponential = slinky
%Looking at it by turning it on its ear:
%Circle (looking in) or sine wave (from the side)
pause
set(gca,'view',[0 0])
pause
set(gca,'view',[90 0])
axis square



