% More reading: https://www.dsprelated.com/showarticle/192.php

clear all
close all

%% Decleration of input Signal, Starts %%%
 fs = 125*10^6; % Sampling frequency (samples per second) 
 dt = 1/fs; % seconds per sample 
 samples = 8000; %8192
 StopTime = (samples-1)*dt; % seconds 
 t_range = (0:dt:StopTime); % seconds 
 f_range = linspace(-fs/2,fs/2,samples); 
  
 F_in = 12.5*10^6; % Sine wave frequency (hertz) 
 A_in = 0.001; % Amplitude (V)
 P_in = pi/5; %pi/5
 signal_in = A_in*cos(2*pi*F_in*t_range + P_in);

 figure
 subplot(4,2,1)
 plot(t_range,signal_in)
 xlim([0 2*10^(-6)])
 ylim([-0.002 0.002]) 
 set(gca,'TickDir','out');
 title('input signal that we want to reconstruct')
 xlabel('time (s)')
 ylabel('Amplitude (V)')

Y = fft(signal_in); %Y is complex
Yshift = fftshift(Y)./samples;
subplot(4,2,3) 
plot(f_range, abs(Yshift)) 
xlim([min(f_range) max(f_range)])
ylim([-0.5*10^-4 12*10^-4])
set(gca,'XTick',[-60*10^6, -40*10^6, -20*10^6, -12.5*10^6,  0, 12.5*10^6, 20*10^6, 40*10^6, 60*10^6])
set(gca,'xticklabel',({'-60', '-40', '-20', '-12.5', '0', '12.5', '20', '40', '60'}))
set(gca,'TickDir','out');
title('Spectrum of input signal')
xlabel('frequency (MHz)')
ylabel('|P1(f)|') 
%%% Decleration of input Signal, Finishes %%%
 
%% Adding Noise and Other Signals on top of our important signal, Starts %%%
DC = 0.0007;
noise = 0.001*((randi(1000,1,samples)-500)./500);
data_s1 = 0.002*cos(2*pi*(3*10^6)*t_range + 0); 
data_s2 = 0.0005*sin(2*pi*(20*10^6)*t_range + pi/4); 
signal_in_noisy = signal_in + data_s1 + data_s2 + noise + DC;
%signal_in_noisy = signal_in;

subplot(4,2,2)
hold on
plot(t_range,signal_in, 'b')
plot(t_range,signal_in_noisy, 'r')
xlim([0 2*10^(-6)])
ylim([-0.01 0.01]) 
set(gca,'TickDir','out');
set(gca,'Box','on');
title('our important signal (blue) with added noise, DC and two other signals (red)')
xlabel('time (s)')
ylabel('Amplitude (V)')
hold off 
 
Y = fft(signal_in_noisy); %Y is complex
Yshift = fftshift(Y)./samples;
subplot(4,2,4)
plot(f_range, abs(Yshift)) 
xlim([min(f_range) max(f_range)])
ylim([-0.5*10^-4 12*10^-4]) 
set(gca,'XTick',[-60*10^6, -40*10^6, -20*10^6, -12.5*10^6,  0, 12.5*10^6, 20*10^6, 40*10^6, 60*10^6])
set(gca,'xticklabel',({'-60', '-40', '-20', '-12.5', '0', '12.5', '20', '40', '60'}))
set(gca,'TickDir','out');
title('Spectrum of input noisy signal')
xlabel('frequency (MHz)')
ylabel('|P1(f)|')  
%%% Adding Noise and Other Signals on top of our important signal, Finishes %%%
 
%% IQ Demodulation Starts %%%% 
 F_LO = F_in; % Sine wave frequency (hertz) 
 A_LO = 1.000; % Amplitude (V)
 P_LO = 0;
 I_LO_signal = 2*A_LO*cos(2*pi*F_LO*t_range + P_LO);
 Q_LO_signal = 2*A_LO*sin(2*pi*F_LO*t_range + P_LO);

 I_before_filter = I_LO_signal.*signal_in_noisy;
 Q_before_filter = Q_LO_signal.*signal_in_noisy;
 Signal_IQ = I_before_filter + 1i.*Q_before_filter;
 Y_IQ = fft (Signal_IQ);
 Y_IQshift = fftshift(Y_IQ)./samples;
 
 subplot(4,2,5)
 plot(t_range, I_LO_signal)
 xlim([0 1*10^(-6)]) 
 ylim([-2.5 2.5])
 set(gca,'TickDir','out');
 title('LO signal')
 xlabel('time (s)')
 ylabel('Amplitude (V)')
  
 subplot(4,2,6)
 plot(f_range, abs(Y_IQshift))  
 xlim([min(f_range) max(f_range)])
 ylim([-0.1*10^-3 2.5*10^-3]) 
 set(gca,'XTick',[-60*10^6, -40*10^6, -20*10^6, -12.5*10^6,  0, 12.5*10^6, 20*10^6, 40*10^6, 60*10^6])
 set(gca,'xticklabel',({'-60', '-40', '-20', '-12.5', '0', '12.5', '20', '40', '60'}))
 set(gca,'TickDir','out');
 title('Spectrum after Digital Down Conversion')
 xlabel('frequency (MHz)')
 ylabel('|P1(f)|')
  %%%% IQ Demodulation Finishes %%%% 
  
%% Boxcar Filter Starts %%%%%%%
 windowSize = 128; 
 a = 1; %normalization coefficient
 b = (1/windowSize)*ones(1,windowSize); %filter coefficients
 I_boxcar = filter(b,a,I_before_filter);
 Q_boxcar = -1.*(filter(b,a,Q_before_filter));
 Signal_IQ_boxcar = I_boxcar + 1i.*Q_boxcar;
 Y_IQ_boxcar = fft (Signal_IQ_boxcar);
 Y_IQ_boxcarshift = fftshift(Y_IQ_boxcar)./samples;
  
 subplot(4,2,8)
 plot(f_range, abs(Y_IQ_boxcarshift))
 xlim([min(f_range) max(f_range)])
 ylim([-0.1*10^-3 2*10^-3]) 
 set(gca,'XTick',[-60*10^6, -40*10^6, -20*10^6, -12.5*10^6,  0, 12.5*10^6, 20*10^6, 40*10^6, 60*10^6])
 set(gca,'xticklabel',({'-60', '-40', '-20', '-12.5', '0', '12.5', '20', '40', '60'}))
 set(gca,'TickDir','out');
 title('Spectrum after LowPass Filter')
 xlabel('frequency (MHz)')
 ylabel('|P1(f)|')
  
     %Below is to see the fft of Boxcar, it is not used on input signal
     filter_in_time= [ones(1,windowSize) zeros(1,(samples-windowSize))];
     Y_filter_in_time = fft (filter_in_time);
     Y_filter_fftshift = fftshift(Y_filter_in_time)./samples;
     subplot(4,2,7)
     plot(f_range, abs(Y_filter_fftshift))
     xlim([min(f_range) max(f_range)])
     ylim([-0.001 0.02]) 
     set(gca,'XTick',[-60*10^6, -40*10^6, -20*10^6, -12.5*10^6,  0, 12.5*10^6, 20*10^6, 40*10^6, 60*10^6])
     set(gca,'xticklabel',({'-60', '-40', '-20', '-12.5', '0', '12.5', '20', '40', '60'}))
     set(gca,'TickDir','out');
     title('Spectrum of LowPass Filter')
     xlabel('frequency (MHz)')
     ylabel('|P1(f)|')
%%%%%%% Boxcar Filter Finishes %%%%%%%
  
%% Reconstruct Signal Starts %%%%%%%
 A_reconstruct = sqrt(I_boxcar.^2 + Q_boxcar.^2);
 P_reconstruct = atan2(Q_boxcar,I_boxcar);
 
 recovered_signal = A_reconstruct.*cos(2*pi*F_LO*t_range + P_reconstruct);
 
 figure
 subplot(5,1,1)
 hold on
 plot(t_range, I_boxcar, 'b')
 plot(t_range, Q_boxcar, 'r')
 xlim([min(t_range) max(t_range)])
 set(gca,'TickDir','out');
 set(gca,'Box','on');
 title('I (blue) and Q (red) signals')
 xlabel('time (s)')
 ylabel('Amplitude (V)')
 hold off
 
 subplot(5,1,2)
 hold on
 plot(t_range, A_reconstruct,'r')
 xlim([min(t_range) max(t_range)])
 ylim([-0.002 0.002])  
 set(gca,'TickDir','out');
 set(gca,'Box','on');
 title('the measured amplitude')
 xlabel('time (s)')
 ylabel('Amplitude (V)')
 hold off
 
 subplot(5,1,3)
 hold on
 plot(t_range, P_reconstruct.*180/pi,'r')
 xlim([min(t_range) max(t_range)])
 ylim([-180 180]) 
 set(gca,'TickDir','out');
 set(gca,'Box','on');
 title('the measured phase')
 xlabel('time (s)')
 ylabel('Phase (degree)')
 hold off
 
 subplot(5,1,4)
 hold on
 plot(t_range, recovered_signal,'r')
 plot(t_range, signal_in,'b')
 xlim([min(t_range) max(t_range)])
 ylim([-0.002 0.002]) 
 set(gca,'TickDir','out');
 set(gca,'Box','on');
 title('input signal and reconstruct signal')
 xlabel('time (s)')
 ylabel('Amplitude (V)')
 hold off
 
 subplot(5,1,5)
 hold on
 plot(t_range, recovered_signal,'r')
 plot(t_range, signal_in,'b')
 xlim([10*10^(-6) 12*10^(-6)])
 ylim([-0.002 0.002]) 
 set(gca,'TickDir','out');
 set(gca,'Box','on');
 title('input signal and reconstruct signal (Zoom)')
 xlabel('time (s)')
 ylabel('Amplitude (V)')
 hold off
%%%%%%% Reconstruct Signal Finishes %%%%%%%


