clear; 
%load('C:\Users\jaanp\Documents\Work\Teaching\Microwave EO Instrumentation\workshops\jaan_SDR_Radar.mat')
load(['CW_Doppler.mat'])
%load('D:\ancortek\Workshop\CW_Doppler.mat')
%load('D:\ancortek\mittaukset\cw1.mat')
whos

c = 299792458;                                              % speed of light (m/s) 
Fs = samplenumberpersweep/(SWEEPTIME/1000)                  %divided by 1000 because sweeptime is given in milliseconds
% Sampling frequency (Hz) 
L = size(DATA,2);                                           % Length of signal (# of samples)
t = (0:L-1)*(1/Fs);                                         % Time vector (s), just increasing vector by 1/samplefrequency


% take only the meaningful part of the data from the beginning
D=DATA(1:L);
% make the data 0-centered in terms of amplitude
M=mean(D);
D=D-M;

% plot I- and Q- parts of the data just to visualize 
figure(1);
plot(t,real(D),'r-');hold on;
plot(t,imag(D),'b-');
ylabel('relative amplitude');
xlabel('time (s)');


% make Fast Fourier Transform to entire measurement
Y=fft(D);
YS=fftshift(Y);               % and center zero by rearranging the data
P=abs(YS)/L;                  % normalize the amplitude by deviding by amount of measurements
fft_frequency = (-L/2:L/2-1)*(Fs/L); % Create frequency vector by knowing that fft bandwidth is sampling frequency (- Fs/2 to + Fs/s)



% visualize the frequency spectrum of the entire measurement
figure(2)
plot(fft_frequency,P,'-') 
title('Frequency Spectrum of entire measurement')
xlabel('f (Hz)')
ylabel('|P(f)|')
M=max(P(:));
axis([-100,100,0,M]); % let's window the x-axis close to +/- 100 Hz as 30 Hz is the dopler of a walking person

% let's divide the data vector to pieces. 


Ltimewindow=8000; % length of the time window for one fft calculation. 

Fres_Hz=Fs/Ltimewindow
timewindow_lenght_in_seconds=Ltimewindow/Fs

Rows=Ltimewindow;
Columns=floor((L/Ltimewindow)); %calculating, how many columns fits into the measurement with given timewindow size

D2=D(1:Rows*Columns);        % cutting lenght of measurement vector which fits into Rows and Columns
DM=reshape(D2,Rows,Columns);   %folding the measurement vector into rows and columns so that every column has a lengt of timewindow

figure(3)
imagesc(t(1:Columns)*Ltimewindow,t(1:Rows),abs(DM)); %  showin
title('Measured time signal folded to rows (time window) and columns ')
ylabel('time in window [s]')
xlabel('measurement time [s]')


% lets now make fft for every time window, which is done automatically if
% data is in rows anc columns

YM = fft(DM);         % FFT
YMS = fftshift(YM);   % putting 0 Hz to the center
fft_frequency2 = (-Ltimewindow/2:Ltimewindow/2-1)*(Fs/Ltimewindow)    % making the frequency vector, -half of the sampling frequency to             
                                                                      % plus half of the sampling frequecny
imagesc(t(1:Columns)*Ltimewindow,fft_frequency2,(abs(YMS)))             % showing the frequency matrix with frequency and time axis
xlabel('measurement time [s]')
ylabel('Frequency')
ylim([-120,120])
% remember, that our frequency resolution 



[a,b]=max(YMS);      %find max frequency in every timewindow (column).
                     %  in vector a is maximum value (not interesting here)
                     % in vector b we have matrix indec, where the max
                     % value is, this is frequency (and speed)
                     
figure(5)
Doppler=fft_frequency2(b);      % fft_frequency2 is the bin actual frequency values, so we make vector where b is pointer
plot(t(1:Columns)*Ltimewindow,Doppler,'*')
title('Doppler in every time winow')
xlabel('time (s)')
ylabel('Doppler frequency (Hz)')

figure(6)
Speed=(Doppler/CENTERFREQUENCY)*c;    % speed is Doppler freq divided by center frequency time speed of light
plot(t(1:Columns)*Ltimewindow,Speed,'r*')
title('Speed in every time winow')
xlabel('time (s)')
ylabel('Speed (m/s)')
                     