clear; clc;

%%
param = parametersDeclaration();

%% Reading csv file

[Idata_all, Qdata_all] = dataFile(param);

%% Balancing IQ signals using S.W.Ellingson algorithm

[Idataall_bal, Qdataall_bal]=SWEllingson_(Idata_all,Qdata_all);

%reshaping the data for signal processing
counter = 0;
for i=1:param.datalength  %datalength

    for j=1:param.numSamplesperChirp % samples per chirp
        
        counter = counter+1;
        Idata_OP(i,j) = Idataall_bal(counter);
        Qdata_OP(i,j) = Qdataall_bal(counter);

    end

end

%% PLOTTING Balanced I AND Q SIGNALS

figure
t = (0:length(Idata_OP(1,:))-1)*param.numSamplesperChirp/length(Idata_OP(1,:));
plot(t,Idata_OP(1,:))
title("Balanced IQ signals")
hold on
t = (0:length(Qdata_OP(1,:))-1)*param.numSamplesperChirp/length(Qdata_OP(1,:));
plot(t,Qdata_OP(1,:))
xlabel("samples per chirp")
ylabel("amplitude")
legend('Idata','Qdata')
hold off

%% Balanced Constellation Diagram

figure
plot(Idataall_bal, Qdataall_bal, 'r*');
title('Balanced IQ constellation diagram');
xlabel('real part');
ylabel('imaginary part');

%% Range FFT

[rangeProfile, peakAmplitudes, maxFrequencies] = rangeFFT(Idata_OP,Qdata_OP,param);

%% object range frequecy for each chirp

[peakFreq] = maxFreq(peakAmplitudes, maxFrequencies);

%% Range for every chirp using the respective beat signal
[target_range] = Range_(peakFreq, rangeProfile, param);

%% Rangebin to be found

cmplx = complex(Idata_OP,Qdata_OP);
[rangebin, databin] = rangeBin_(target_range, rangeProfile, cmplx, param);

%% phase and phase-based range measurement
% 6000 chirps were trasmitted which means 6000/20 number of frames which is 
% the same number of seconds the radar was active. After doing the math, the radar 
% recorded for 300 seconds. Taking 1200 samples covers 60 seconds which is used 
% for measuring the phase and phase-based range measurement.

[phi, phase_range] = extendedDACMalgorithm_(rangebin, param);

%% Phase difference

phi_new = phaseDifference(phi,param);

phi_cut = phi_new(param.startsample:param.endsample);
phi_cut = ImpulseDenoise(phi_cut,param);


%% Filtering both HeartBeat and BreathingRate out

[filterBR, filterHR] = filterDesign_(param);

breathing_signal = filter(filterBR,phi_new);
heart_signal = filter(filterHR,phi_new);

figure
subplot(2,1,1)
plot(linspace(param.startsec,param.endsec,param.samplerange),breathing_signal(param.startsample:param.endsample))
title("Breathing Signal vs Time")
xlabel("time(sec)")
ylabel("Amplitude")
subplot(2,1,2)
plot(linspace(param.startsec,param.endsec,param.samplerange),heart_signal(param.startsample:param.endsample))
title("Heart Signal vs Time")
xlabel("time(sec)")
ylabel("Amplitude")


%% Autocorrelation estimation in time domain

[estimatedBRauto_filter, estimatedHRauto_filter] = autoEstimation(breathing_signal(param.startsample:param.endsample), heart_signal(param.startsample:param.endsample),param);


%% Compressing and Reconstructing the Breathing Signal Using CS-OMP algorithm

[reconstructed_breathingSignal, brPeakAmps, brPeakFreqs] = signalReconstruction_(breathing_signal(param.startsample:param.endsample),param);
[reconstructed_heartSignal, hrPeakAmps, hrPeakFreqs] = signalReconstruction_(heart_signal(param.startsample:param.endsample),param, 20);

[estimatedBRauto_CS, estimatedHRauto_CS] = autoEstimation(reconstructed_breathingSignal, reconstructed_heartSignal,param);
[estimatedBRfft_CS, estimatedHRfft_CS] = CalculateCS_VitalSigns(brPeakAmps,brPeakFreqs,hrPeakAmps,hrPeakFreqs);


%% Use waveletAnalyzer toolbox to decompose and reconstruct the signals just using sinusoidal wave for easy approximation 

Wavelet_reconstructionTool(breathing_signal,1,param);
Wavelet_reconstructionTool(heart_signal,10,param);


%% Wavelet Decomposition

[peakAmp, locFreq]=waveletTrial(phi_cut,param);
 
[estimatedBR_WT, estimatedHR_WT] = Calculate_VSWT(peakAmp,locFreq);