
% the script is made for Matlab, but it works without modification as well in Octave

load('FMCW_sawtooth.mat');

DATE
WAVEFORM
BANDWIDTH                 % bandwidth (Hz)
CENTERFREQUENCY           % frequency of the radar (Hz)
samplenumberpersweep      % samples per sweep 
SWEEPTIME=SWEEPTIME/1000  % sweep time (s) (the value in the file is given in milliseconds, therefore /1000)
speedoflight=299792458    % (m/s)
DATA=DATA(256*3000:256*9000-1); % I shorten here the data vector, otherwise Ocatave imagesc does not work

sweeps=floor(length(DATA)./samplenumberpersweep)
Fs = samplenumberpersweep/(SWEEPTIME) %sampling frequency

B=DATA-mean((DATA(:)));    %Remove the digitalization offset

% Reshape the data to stack the sweeps in columns

D=reshape(B,samplenumberpersweep,floor(sweeps)); %arrange every sweep as one column

% create time vector by knowing how many samples we had in sweep lenght
sweeptimevector=(1:samplenumberpersweep).*(SWEEPTIME./samplenumberpersweep);

% Show the waveform of a sweep, time on x-axis, real part on y-axis. 
figure;
plot(sweeptimevector,real(D(:,500))); % plotting the 
xlabel('time (s)');
ylabel('real (I) part');
title('Real (I) of the sweep 500');

%Show the time domain data waterfall (sweeps as columns) 
figure;
imagesc(real(D))
xlabel(['sweep number']);
ylabel('sample');
title('Measured DATA real part');


% make the Fast Fourier transform for all the data
% note that fft acts column-wise, in the resulting matrix you have instead
% of amplitudes, the sweep frequencies as columns


F=abs(fft(D));
FF=fftshift(F,1); % change left and right part of FFT output, to get 0 -frequency in the midle


% Show the FFT transformed frequency map
% create time vector (start time for every sweep)
time=(0:sweeps-1)*SWEEPTIME;

% create frequencies vector, - half of sampling frequency to + half of
% sampling frequency
samplingfrequencystep=Fs/samplenumberpersweep;
frequencies=(-(samplenumberpersweep/2):((samplenumberpersweep/2)-1))*samplingfrequencystep;

% show the frequency waterfall map, time on x and frequencies on y axis

figure;
imagesc(time,frequencies,FF);
xlabel('time (s)');
ylabel('frequency (Hz)');
title('Fourier transform of the measured DATA');


% Show the frequency spectrunm of one sweep

figure;
plot(frequencies,FF(:,500));
xlabel('frequency (Hz)');
ylabel('frequency amplitude');
title('Frequency spectrum of sweep 500');

% find max frequency in every column
[a,b]=max(FF);  % note that a gives max amplitud and b is the index where max was found, this tells the frequency

% calculate range calibration constand for easy calculation of distance
rangecalib=(speedoflight*SWEEPTIME)./(2*BANDWIDTH);
% calculate distance
distance=frequencies(b)*rangecalib;

% Plot the distance as function of starting time of every sweep
figure;
plot(time,distance,'*')
xlabel('time (s)');
ylabel('range (m)');
title('Distance of the target from the radar');

% Note, that the distance measurement seems to have poor resolution.
% Unfortunately this cannot be avoided. By studying the distance
% calculation equation and taking into account the frequency resolution
% % equation one can see, that distance resolution for this type of radar is
% deltaD=speedoflight/(2*bandwidth), independently of other radar
% parameters. In our case, for 400 MHz, the deltaD=37 cm.
% therefore, smoother curves are only possible with greater distances







