% Script to compute the Relative Intensity Noise of a laser by recording the y-t signal
% by Mathias Neidig 2012_09_11

% The RIN is defined as 
%
%   RIN = 10* log10 [Single-sided power spectrum density / (average power)]
%
% and is given in [RIN] = dB/Hz

clear all
close all

%% Set the directory where the data is

dirDCData = ['C:\\Users\\Karthik\\Documents\\GitRepositories\\Calculations\\Time-Series-Analyzer\\Time-Series-Data\\20240807\\DC Coupling\\'];
dirACData = ['C:\\Users\\Karthik\\Documents\\GitRepositories\\Calculations\\Time-Series-Analyzer\\Time-Series-Data\\20240807\\AC Coupling\\'];

%% Load the files which contain: - the DC coupled y-t signal to obtain the averaged power
%                                 - the AC coupled y-t signal to obtain the fluctuations
%                                 - the AC coupled y-t signal with the beam blocked to obtain the background fluctuations

%-------------------------------------------------------------------------%
dcsignal = readmatrix( [ dirDCData 'P7.0_M3.0_OOL.csv'] );                %
acsignal = readmatrix( [ dirACData 'P7.0_M3.0_OOL.csv'] );                %
bgsignal = readmatrix( [ dirACData 'Bkg_OOL.csv'] );                      %
%-------------------------------------------------------------------------%

%% Read out the important parameters

time_increment = 2E-6;
dctime         = dcsignal(1:end, 1) .* time_increment;
actime         = acsignal(1:end, 1) .* time_increment;
bgtime         = bgsignal(1:end, 1) .* time_increment;

dcdata         = dcsignal(1:end, 2);
acdata         = acsignal(1:end, 2);
bgdata         = bgsignal(1:end, 2);

N       = length(actime);               % #samples
f_s     = 1/time_increment;             % Sample Frequency
delta_f = f_s/N;                        % step size in frequency domain
delta_t = 1/f_s;                        % time step

%% Custom Control Parameters

% Choose smoothing parameter; has to be odd
%----------------%  
    span = 21;   %  
%----------------%

%% Computes the RIN

% compute the average power (voltage^2)
average_P = mean(dcdata.*dcdata);  

% compute the power spectrum density FFT(A) x FFT*(A)/N^2 of the source & the bg
psd_src   = fft(acdata) .* conj(fft(acdata))/N^2;
psd_bg    = fft(bgdata) .* conj(fft(bgdata))/N^2;

% converts the psd to the single-sided psd --> psd is symmetric around zero --> omit
% negative frequencies and put the power into the positive ones --> spsd

for i = 1 : N/2+1
    if i>1
         spsd_src(i) = 2*psd_src(i);
         spsd_bg(i)  = 2*psd_bg(i);
    else spsd_src(i) = psd_src(i);
         spsd_bg(i)  = psd_bg(i);
    end
end

% smooths the spsd by doing a moving average
spsd_src_smooth = smooth(spsd_src,span,'moving');
spsd_bg_smooth  = smooth(spsd_bg, span,'moving');

% calculates the RIN given in dB/Hz; the factor delta_f is needed to convert from dB/bin into dB/Hz
RIN_src_smooth = 10*log10(spsd_src_smooth/(average_P*delta_f));
RIN_bg_smooth  = 10*log10(spsd_bg_smooth /(average_P*delta_f));

% creates an array for the frequencies up to half the sampling frequency
f        = f_s/2 * linspace(0,1,N/2+1);
f_smooth = smooth(f,span,'moving');

% Plots the RIN vs frequency
f_ = clf;
figure(f_);
semilogx(f_smooth,RIN_bg_smooth,'k-')
hold on
semilogx(f_smooth,RIN_src_smooth,'r-')
xlabel('Frequency [Hz]')
ylabel('RIN [dB/Hz]')
xlim([min(f) max(f)]);
title('\bf Relative Intensity Noise of ODT Arm 1')
legend('Background PD Box', 'Power:7 V, Mod:-3.0 V','Location','NorthEast');
% text(1e5,-95,['\bf MovingAverage = ' num2str(span) ]);
grid on
% optional: save the picture without editing wherever you want
%------------------------------------------%     
%     saveas(f_,'FileName','png');         %
%------------------------------------------%