From 557d84649e8b209258d4665eeff93682b7b6428e Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sun, 15 Sep 2024 17:02:06 +0200 Subject: [PATCH] Corrected and updated analysis scripts. --- .../extractIRF.m | 89 +++++++++++-------- Time-Series-Analyzer/computeRIN.m | 79 +++++++++------- 2 files changed, 99 insertions(+), 69 deletions(-) diff --git a/Imaging-Response-Function-Extractor/extractIRF.m b/Imaging-Response-Function-Extractor/extractIRF.m index b98c5b0..74b9918 100644 --- a/Imaging-Response-Function-Extractor/extractIRF.m +++ b/Imaging-Response-Function-Extractor/extractIRF.m @@ -1,31 +1,35 @@ %% Parameters -groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", "/images/Vertical_Axis_Camera/in_situ_absorption"]; +groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", "/images/Vertical_Axis_Camera/in_situ_absorption"]; -folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/IRF/0044/"; +folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Imaging-Response-Function-Extractor/0127/"; -cam = 5; +cam = 5; -angle = 90 + 51.5; -center = [1700, 2300]; -span = [255, 255]; -fraction = [0.1, 0.1]; +% angle = 90 + 51.5; +% center = [1700, 2300]; +angle = 90; +center = [2035 1250]; +span = [255, 255]; +fraction = [0.1, 0.1]; NA = 0.6; pixel_size = 4.6e-6; lambda = 421e-9; -d = lambda/2/pi/NA; -k_cutoff = NA/lambda/1e6; +d = 0.61*lambda/NA; +k_cutoff = 2*NA/lambda/1e6; % in units of 1/µm) + +removeFringe = true; %% Compute OD image, rotate and extract ROI for analysis % Get a list of all files in the folder with the desired file name pattern. + filePattern = fullfile(folderPath, '*.h5'); files = dir(filePattern); refimages = zeros(span(1) + 1, span(2) + 1, length(files)); absimages = zeros(span(1) + 1, span(2) + 1, length(files)); - for k = 1 : length(files) baseFileName = files(k).name; fullFileName = fullfile(files(k).folder, baseFileName); @@ -40,15 +44,24 @@ for k = 1 : length(files) absimages(:,:,k) = subtract_offset(crop_image(calculate_OD(atm_img, bkg_img, dark_img), center, span), fraction); end + %% Fringe removal -optrefimages = removefringesInImage(absimages, refimages); -absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :); - -nimgs = size(absimages_fringe_removed,3); -od_imgs = cell(1, nimgs); -for i = 1:nimgs - od_imgs{i} = absimages_fringe_removed(:, :, i); +if removeFringe + optrefimages = removefringesInImage(absimages, refimages); + absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :); + + nimgs = size(absimages_fringe_removed,3); + od_imgs = cell(1, nimgs); + for i = 1:nimgs + od_imgs{i} = absimages_fringe_removed(:, :, i); + end +else + nimgs = size(absimages(:, :, :),3); + od_imgs = cell(1, nimgs); + for i = 1:nimgs + od_imgs{i} = absimages(:, :, i); + end end %% Compute the Density Noise Spectrum @@ -63,8 +76,8 @@ density_noise_spectrum = cell(1, length(od_imgs)); dx = pixel_size; dy = pixel_size; -xvals = (1:Nx)*dx*1e6; -yvals = (1:Ny)*dy*1e6; +xvals = (1:Nx)*dx*1e6; +yvals = (1:Ny)*dy*1e6; Nyq_k = 1/dx; % Nyquist dk = 1/(Nx*dx); % Wavenumber increment @@ -77,15 +90,15 @@ ky = -Nyq_k/2:dk:Nyq_k/2-dk; % wavenumber ky = ky * dy; % wavenumber (in units of 1/dy) % Create Circular Mask -n = 2^8; % size of mask -mask = zeros(n); -I = 1:n; -x = I-n/2; % mask x-coordinates -y = n/2-I; % mask y-coordinates -[X,Y] = meshgrid(x,y); % create 2-D mask grid -R = 32; % aperture radius -A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R -mask(A) = 1; % set mask elements inside aperture to 1 +n = 2^8; % size of mask +mask = zeros(n); +I = 1:n; +x = I-n/2; % mask x-coordinates +y = n/2-I; % mask y-coordinates +[X,Y] = meshgrid(x,y); % create 2-D mask grid +R = 32; % aperture radius +A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R +mask(A) = 1; % set mask elements inside aperture to 1 % Calculate Power Spectrum and plot @@ -106,7 +119,7 @@ for k = 1 : length(od_imgs) ylabel('µm', 'FontSize', 16) axis equal tight; colorbar - colormap (flip(jet)); + colormap jet; % (flip(jet)) % set(gca,'CLim',[0 100]); set(gca,'YDir','normal') title('Single-shot image', 'FontSize', 16); @@ -119,7 +132,7 @@ for k = 1 : length(od_imgs) ylabel('µm', 'FontSize', 16) axis equal tight; colorbar - colormap (flip(jet)); + colormap jet; % (flip(jet)) % set(gca,'CLim',[0 100]); set(gca,'YDir','normal') title('Averaged density image', 'FontSize', 16); @@ -132,7 +145,7 @@ for k = 1 : length(od_imgs) ylabel('µm', 'FontSize', 16) axis equal tight; colorbar - colormap (flip(jet)); + colormap jet; % (flip(jet)) % set(gca,'CLim',[0 100]); set(gca,'YDir','normal') title('Image noise = Single-shot - Average', 'FontSize', 16); @@ -145,7 +158,7 @@ for k = 1 : length(od_imgs) ylabel('µm', 'FontSize', 16) axis equal tight; colorbar - colormap (flip(jet)); + colormap jet; % (flip(jet)) % set(gca,'CLim',[0 100]); set(gca,'YDir','normal') title('Masked Noise', 'FontSize', 16); @@ -158,7 +171,7 @@ for k = 1 : length(od_imgs) ylabel('1/dy', 'FontSize', 16) axis equal tight; colorbar - colormap (flip(jet)); + colormap jet; % (flip(jet)) % set(gca,'CLim',[0 100]); set(gca,'YDir','normal') title('DFT', 'FontSize', 16); @@ -171,7 +184,7 @@ for k = 1 : length(od_imgs) ylabel('1/dy', 'FontSize', 16) axis equal tight; colorbar - colormap (flip(jet)); + colormap jet; % (flip(jet)) % set(gca,'CLim',[0 100]); set(gca,'YDir','normal') title('Density Noise Spectrum = |DFT|^2', 'FontSize', 16); @@ -256,7 +269,7 @@ end subplot('Position', [0.55, 0.1, 0.4, 0.8]) % Adjusted position % [r, Zr] = radial_profile(averagePowerSpectrum, 1); -% Zr = (Zr - min(Zr))./(max(Zr) - min(Zr)); +% Zr = (Zr - min(Zr))./(max(Zr) - min(Zr)); % plot(r, Zr, 'o-', 'MarkerSize', 4, 'MarkerFaceColor', 'none'); % set(gca, 'XScale', 'log'); % Setting x-axis to log scale @@ -267,7 +280,7 @@ ks = sqrt(kx.^2 + ky.^2); profile = profile(length(profile)/2:end); ks = ks(length(ks)/2:end); -n = 0.15; +n = 0.05; [val,slice_idx]=min(abs(ks-n)); ks = ks(1:slice_idx); profile = profile(1:slice_idx); @@ -275,7 +288,9 @@ plot(ks, profile, 'b*-'); % plot(profile, 'b*-'); grid on; % xlim([min(ks) max(ks)]) -title('Radial average of Density Noise Spectrum', 'FontSize', 16); +xlabel('k (1/µm)', 'FontSize', 16) +ylabel('Normalised amplitude', 'FontSize', 16) +title('Radial profile', 'FontSize', 16); grid on; diff --git a/Time-Series-Analyzer/computeRIN.m b/Time-Series-Analyzer/computeRIN.m index 427a696..d916fd5 100644 --- a/Time-Series-Analyzer/computeRIN.m +++ b/Time-Series-Analyzer/computeRIN.m @@ -1,5 +1,6 @@ % Script to compute the Relative Intensity Noise of a laser by recording the y-t signal -% by Mathias Neidig 2012_09_11 +% by Mathias Neidig in 2012 +% modified for DyLab use by Karthik in 2024 % The RIN is defined as % @@ -12,32 +13,33 @@ 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\\']; +dirDCData = 'C:\\Users\\Karthik\\Documents\\GitRepositories\\Calculations\\Time-Series-Analyzer\\Time-Series-Data\\20240915\\DC_Coupling\\'; +dirACData = 'C:\\Users\\Karthik\\Documents\\GitRepositories\\Calculations\\Time-Series-Analyzer\\Time-Series-Data\\20240915\\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'] ); % +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 -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.A; +acdata = acsignal.A; +bgdata = bgsignal.A; +picodata = picosignal.A; -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 +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 @@ -56,31 +58,37 @@ 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; +psd_pico = fft(picodata) .* conj(fft(picodata))/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); + spsd_src(i) = 2*psd_src(i); + spsd_bg(i) = 2*psd_bg(i); + spsd_pico(i) = 2*psd_pico(i); + else + spsd_src(i) = psd_src(i); + spsd_bg(i) = psd_bg(i); + spsd_pico(i) = psd_pico(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'); +spsd_src_smooth = smooth(spsd_src,span,'moving'); +spsd_bg_smooth = smooth(spsd_bg, span,'moving'); +spsd_pico_smooth = smooth(spsd_pico, 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)); +RIN_src_smooth = 10*log10(spsd_src_smooth/(average_P*delta_f)); +RIN_bg_smooth = 10*log10(spsd_bg_smooth /(average_P*delta_f)); +RIN_pico_smooth = 10*log10(spsd_pico_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'); - +% % Calculates the shot noise limit of the used PD given the wavelength of the light source and % incident average power PlanckConstant = 6.62607015E-34; @@ -88,21 +96,28 @@ SpeedOfLight = 299792458; WavelengthOfLaserLight = 1064E-9; FrequencyOfLaserLight = SpeedOfLight / WavelengthOfLaserLight; QuantumEfficiencyOfPD = 1; -ShotNoiseLimit = 10*log10((2 * PlanckConstant * FrequencyOfLaserLight / QuantumEfficiencyOfPD) * average_P); +AverageIncidentPower = 0.001; % (in W) +ShotNoiseLimit = 10*log10((2 * PlanckConstant * FrequencyOfLaserLight) / (QuantumEfficiencyOfPD * AverageIncidentPower)); + %% Plots the RIN % Plots the RIN vs frequency f_ = clf; figure(f_); -semilogx(f_smooth,RIN_bg_smooth,'k-') +semilogx(f_smooth, RIN_pico_smooth, LineStyle = "-", Color = [.7 .7 .7]) hold on +semilogx(f_smooth, RIN_bg_smooth, LineStyle = "-", Color = [.0 .0 .0]) semilogx(f_smooth,RIN_src_smooth,'r-') yline(ShotNoiseLimit,'--b'); -xlabel('Frequency [Hz]') -ylabel('RIN [dB/Hz]') -xlim([10 max(f)]); -title('\bf Relative Intensity Noise of ODT Arm 1') -legend('Detector Noise', 'Power:7 V, Mod: 100%, with PID ON', 'Shot-Noise limit','Location','NorthWest'); +ax = gca(f_); +ax.XAxis.FontSize = 14; +ax.YAxis.FontSize = 14; +xlabel('Frequency [Hz]', FontSize=16) +ylabel('RIN [dBc/Hz]', FontSize=16) +xlim([10 5E6]); +ylim([-175 -55]); +title('\bf Relative Intensity Noise of IPG 1064', FontSize=16) +legend(label_0, label_1, label_2, label_3, 'Location','NorthEast', FontSize=16); % text(1e5,-95,['\bf MovingAverage = ' num2str(span) ]); grid on