Corrected and updated analysis scripts.
This commit is contained in:
parent
235a436436
commit
557d84649e
@ -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;
|
||||
|
||||
|
||||
|
@ -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
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user