diff --git a/Time-Series-Analyzer/computeAndcompareRIN.m b/Time-Series-Analyzer/computeAndcompareRIN.m index 1083007..c0002bc 100644 --- a/Time-Series-Analyzer/computeAndcompareRIN.m +++ b/Time-Series-Analyzer/computeAndcompareRIN.m @@ -22,13 +22,13 @@ dirACData = 'C:\\Users\\Karthik\\Documents\\GitRepositories\\Calculations\\Time- %------------------------------------------------------------------------- % bgsignal = load( [ dirACData '20240915_Background_AC'] ); % -dcsignal = load( [ dirDCData '20240915_1V_-3Mod_DC_PID_Inactive'] ); % -acsignal_1 = load( [ dirACData '20240915_1V_-3Mod_AC_PID_Inactive'] ); % -acsignal_2 = load( [ dirACData '20240915_1V_-3Mod_AC_PID_Active'] ); % +dcsignal = load( [ dirDCData '20240915_5V_1.1Mod_DC_PID_Inactive'] ); % +acsignal_1 = load( [ dirACData '20240915_5V_1.1Mod_AC_PID_Inactive'] ); % +acsignal_2 = load( [ dirACData '20240915_5V_1.1Mod_AC_PID_Active'] ); % label_0 = 'Background'; % -label_1 = 'Power = 1 V, Modulation = -3.0 V, PID Inactive'; % -label_2 = 'Power = 1 V, Modulation = -3.0 V, PID Active'; % +label_1 = 'Power = 5 V, Modulation = 1.1 V, PID Inactive'; % +label_2 = 'Power = 5 V, Modulation = 1.1 V, PID Active'; % %------------------------------------------------------------------------- % %% Read out the important parameters @@ -38,10 +38,10 @@ acdata_1 = acsignal_1.A; acdata_2 = acsignal_2.A; bgdata = bgsignal.A; -N = length(dcdata); % #samples -f_s = 1/dcsignal.Tinterval; % Sample Frequency -delta_f = f_s/N; % step size in frequency domain -delta_t = 1/f_s; % time step +N = length(dcdata); % #samples +f_s = 1/dcsignal.Tinterval; % Sample Frequency +delta_f = f_s/N; % step size in frequency domain +delta_t = 1/f_s; % time step %% Custom Control Parameters @@ -53,12 +53,12 @@ delta_t = 1/f_s; % time step %% Computes the RIN % compute the average power (voltage^2) -average_P = mean(dcdata.*dcdata); +average_P = mean(dcdata.*dcdata); % compute the power spectrum density FFT(A) x FFT*(A)/N^2 of the source & the bg psd_src_1 = fft(acdata_1) .* conj(fft(acdata_1))/N^2; psd_src_2 = fft(acdata_2) .* conj(fft(acdata_2))/N^2; -psd_bg = fft(bgdata) .* conj(fft(bgdata))/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 @@ -67,23 +67,26 @@ for i = 1 : N/2+1 if i>1 spsd_src_1(i) = 2*psd_src_1(i); spsd_src_2(i) = 2*psd_src_2(i); - spsd_bg(i) = 2*psd_bg(i); + spsd_bg(i) = 2*psd_bg(i); else spsd_src_1(i) = psd_src_1(i); spsd_src_2(i) = psd_src_2(i); - spsd_bg(i) = psd_bg(i); + spsd_bg(i) = psd_bg(i); end end % smooths the spsd by doing a moving average spsd_src_smooth_1 = smooth(spsd_src_1,span,'moving'); spsd_src_smooth_2 = smooth(spsd_src_2,span,'moving'); -spsd_bg_smooth = smooth(spsd_bg, 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_1 = 10*log10(spsd_src_smooth_1/(average_P*delta_f)); -RIN_src_smooth_2 = 10*log10(spsd_src_smooth_2/(average_P*delta_f)); -RIN_bg_smooth = 10*log10(spsd_bg_smooth /(average_P*delta_f)); +spsd_src_smooth_rel_1 = spsd_src_smooth_1/(average_P*delta_f); +spsd_src_smooth_rel_2 = spsd_src_smooth_2/(average_P*delta_f); +spsd_bg_smooth_rel = spsd_bg_smooth /(average_P*delta_f); +RIN_src_smooth_1 = 10*log10(spsd_src_smooth_rel_1); +RIN_src_smooth_2 = 10*log10(spsd_src_smooth_rel_2); +RIN_bg_smooth = 10*log10(spsd_bg_smooth_rel); % creates an array for the frequencies up to half the sampling frequency f = f_s/2 * linspace(0,1,N/2+1); @@ -125,6 +128,68 @@ grid on % saveas(f_,'FileName','png'); % %------------------------------------------% +%% + +compute_eFoldingTime(f_smooth, spsd_src_smooth_rel_2, 100, 1000, label_2) + +%% + +compute_eFoldingTime(f_smooth, spsd_src_smooth_rel_2, 1000, 5000, label_2) + +%% + +function compute_eFoldingTime(faxis, Skk, fstart, fend, data_str) + % computes the energy e-folding time: time to increase the energy by a factor e in sec. + + [~,idx_ini] = min(abs(faxis-(2 * fstart))); % Find closest index for 2*fstart + [~,idx_fin] = min(abs(faxis-(2 * fend))); % Find closest index for 2*fend + Skk = Skk(idx_ini:idx_fin); % Slice Skk array + freqs = linspace(fstart, fend, length(Skk)); + + % Calculate e-folding time + e_folding_time = arrayfun(@(idx) 1 / (pi^2 * freqs(idx)^2 * Skk(idx)), 1:length(freqs)); + + % Plotting + figure('Position', [100, 100, 1200, 800]); + semilogx(freqs, e_folding_time, 'r', 'LineWidth', 2, 'DisplayName', data_str); + hold on; + + grid on; + set(gca, 'GridLineStyle', '-'); % Thin grid lines + + % Create gridlines for multiples of 10 + x_multiples_of_10 = 10.^floor(log10(min(freqs(freqs > 0))):log10(max(freqs(freqs > 0)))); + for val = x_multiples_of_10 + xline(val, 'k-', 'LineWidth', 2); % Thick lines for multiples of 10 + end + + legend('show', 'Location', 'southwest', 'FontSize', 12); + xlabel('Frequency [Hz]', 'FontSize', 14); + ylabel('$T_e$ [s]', 'Interpreter', 'latex', 'FontSize', 14); + title(sprintf('Lower Bound= %.5f s', min(e_folding_time)), 'FontSize', 14); + + set(gca, 'XScale', 'log'); + set(gca, 'FontSize', 12); + + hold off; + +end + + + + + + + + + + + + + + + + diff --git a/Time-Series-Analyzer/computeRIN.m b/Time-Series-Analyzer/computeRIN.m index d916fd5..db9fffc 100644 --- a/Time-Series-Analyzer/computeRIN.m +++ b/Time-Series-Analyzer/computeRIN.m @@ -20,16 +20,16 @@ dirACData = 'C:\\Users\\Karthik\\Documents\\GitRepositories\\Calculations\\Time- % - 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 = load( [ dirDCData 'IPG1064_100W'] ); % -acsignal = load( [ dirACData 'IPG1064_100W'] ); % -bgsignal = load( [ dirACData 'Background'] ); % -picosignal = load( [ dirACData 'Picoscope_Background'] ); % -label_0 = 'Picoscope noise floor'; % -label_1 = 'Background (Picoscope + InGaAs PIN PD + Transimp amp + Buffer amp + power supply)'; % -label_2 = 'Laser output power=100 W (Incident on PD=1mW)'; % -label_3 = 'Shot-Noise limit @ 1 mW incident power'; % -%-------------------------------------------------------------------------% +%--------------------------------------------------------------------------------------------------% +dcsignal = load( [ dirDCData 'IPG1064_100W'] ); % +acsignal = load( [ dirACData 'IPG1064_100W'] ); % +bgsignal = load( [ dirACData 'Background'] ); % +picosignal = load( [ dirACData 'Picoscope_Background'] ); % +label_0 = 'Picoscope noise floor'; % +label_1 = 'Background (Picoscope + InGaAs PIN PD + Transimp amp + Buffer amp + power supply)'; % +label_2 = 'Laser output power=100 W (Incident on PD=1mW)'; % +label_3 = 'Shot-Noise limit @ 1 mW incident power'; % +%------------------------------------------------------------------------------------------------- % %% Read out the important parameters