karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Corrected analysis, added Matthias Niedig script used in Selim Jochim's group.

Browse Source
This commit is contained in:
Karthik 2024-09-04 17:06:13 +02:00
parent c0e96006d7
commit c5fa89ed5d
3 changed files with 411 additions and 572 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

328
Time-Series-Analyzer/Analyzer.py
View File

@ -1,6 +1,6 @@
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec from scipy import integrate
from scipy.signal import find_peaks, resample, detrend from scipy.signal import find_peaks, resample, detrend
import csv import csv
@ -49,13 +49,18 @@ def compute_autocorrelation(data):
acor = np.correlate(yunbiased, yunbiased, "same")/ynorm acor = np.correlate(yunbiased, yunbiased, "same")/ynorm
return acor return acor
def pre_process_data(data, new_sampling_rate): def pre_process_data(data, re_sample = False, new_sampling_rate = 10000):
# Resample the time series # Resample the time series
# Define the new sampling rate and calculate the new time values # Define the new sampling rate and calculate the new time values
n_points = int((len(data[:, 0]) * (data[:, 0][1] - data[:, 0][0])) * new_sampling_rate) if re_sample:
t = np.linspace(data[:, 0][0], data[:, 0][-1], n_points) n_points = int((len(data[:, 0]) * (data[:, 0][1] - data[:, 0][0])) * new_sampling_rate)
x_array = resample(data[:, 1], n_points) t = np.linspace(data[:, 0][0], data[:, 0][-1], n_points)
x = detrend(x_array - x_array.mean()) x_array = resample(data[:, 1], n_points)
x = detrend(x_array - x_array.mean())
else:
t = data[:, 0]
x_array = data[:, 1]
x = detrend(x_array - x_array.mean())
return np.column_stack((t, x)) return np.column_stack((t, x))
@ -63,41 +68,35 @@ def smoothen_data(data, window_size):
# Smooth data by doing a moving average # Smooth data by doing a moving average
return np.convolve(data, np.ones(window_size, dtype=int)/window_size, mode='valid') return np.convolve(data, np.ones(window_size, dtype=int)/window_size, mode='valid')
def compute_psd(data, new_sampling_rate, window_size = 21): def compute_psd(data, re_sample = False, new_sampling_rate = 10000, window_size = 21):
""" """
A power spectral density (PSD) takes the amplitude of the FFT, multiplies it by its complex conjugate and normalizes it to the frequency bin width. A power spectral density (PSD) takes the amplitude of the FFT, multiplies it by its complex conjugate and normalizes it to the frequency bin width.
""" """
processed_data = pre_process_data(data, new_sampling_rate) processed_data = pre_process_data(data, re_sample, new_sampling_rate)
t, x = processed_data[:, 0], processed_data[:, 1] t, x = processed_data[:, 0], processed_data[:, 1]
dt = t[1] - t[0] # Define the sampling interval N = t.shape[0] # Total number of data points
N = x.shape[0] # Define the total number of data points
T = N * dt # Define the total duration of the data
# Calculate fft # Calculate fft
print("Calculating power spectrum...") print("Calculating power spectrum...")
fft_ts = np.fft.fft(x) # Compute Fourier transform of x fft_ts = np.fft.fft(x) # Compute Fourier transform of x
Sxx = 2 * (dt ** 2 / T) * (fft_ts * fft_ts.conj()) # Compute spectrum Sxx = (fft_ts * fft_ts.conj()) / N**2 # Compute spectrum
Sxx = Sxx[:int(len(x) / 2)] # Ignore negative frequencies, we have accounted for this by the scaling factor of 2 in the previous step Sxx = 2 * Sxx[:int(len(x) / 2)] # Ignore negative frequencies, we have accounted for this by the scaling factor of 2
return processed_data, smoothen_data(Sxx.real, window_size) return processed_data, smoothen_data(Sxx.real, window_size)
def compute_RIN(time, voltages, Sxx_smooth): def compute_RIN(dc_voltages, delta_f, Sxx_smooth):
dt = time[1] - time[0] # Define the sampling interval
N = voltages.shape[0] # Define the total number of data points
T = N * dt # Define the total duration of the data
df = 1 / T.max()
# Compute the average power # Compute the average power
average_P = np.mean(np.squared(voltages)) average_P = np.mean(np.square(dc_voltages))
# Calculate the RIN # Calculate the RIN
RIN_Sxx_smooth = 10 * np.log10(Sxx_smooth / (average_P * df)) Skk = Sxx_smooth / (average_P * delta_f)
RIN_Sxx_smooth = 10 * np.log10(Skk)
return RIN_Sxx_smooth return Skk, RIN_Sxx_smooth
def find_noise_peaks(psd, faxis, freq_range, threshold): def find_noise_peaks(psd, faxis, freq_range, threshold):
""" """
@ -215,43 +214,50 @@ def extract_data(filepath):
return data_array return data_array
def plot_analysis(data, data_bkg, Sxx, Sxx_bkg, data_str, bkg_str, peak_find_threshold, window_size = 21, plot_only_psd = True): def closest(list, num):
idx = int(min(range(len(list)), key=lambda i: abs(list[i]-num)))
return idx
def plot_psd(data, Sxx, data_str, peak_find = False, peak_find_threshold = -100, window_size = 21):
time, voltages = data[:, 0], data[:, 1] if isinstance(data, list) and isinstance(Sxx, list) and isinstance(data_str, list):
time_bkg, voltages_bkg = data_bkg[:, 0], data_bkg[:, 1]
dt = time[1] - time[0] # Define the sampling interval
N = voltages.shape[0] # Define the total number of data points
T = N * dt # Define the total duration of the data
df = 1 / T.max()
fNQ = 1 / dt / 2 # Determine Nyquist frequency
faxis = smoothen_data(np.linspace(0,fNQ,N//2), window_size) # Construct frequency axis
""" Noise levels in units of Vrms^2/Hz"""
# resolution_bandwidth = (1/dt)/N
# broadband_noise_level = compute_noise_level(Sxx, resolution_bandwidth) # Integrates across PSD from DC to Nyquist frequency, gives result in in units of Vrms^2/Hz
# noise_floor = np.mean(Sxx_bkg)
freq_range = (50, max(faxis))
threshold = 10**(peak_find_threshold/10)
peak_powers, peak_frequencies = find_noise_peaks(Sxx, faxis, freq_range, threshold)
if plot_only_psd:
plt.figure(figsize=(12, 8)) plt.figure(figsize=(12, 8))
# Plot Power Spectrum in dB for idx in range(len(data)):
plt.semilogx(faxis, 10 * np.log10(Sxx_bkg), color='orange', linewidth=0.5, label = bkg_str)
plt.semilogx(faxis, 10 * np.log10(Sxx), color='green', linewidth=2, label = data_str)
# plt.axhline(y=10 * np.log10(broadband_noise_level), color='red', linewidth=2, linestyle='--', label=f'Broadband cumulative noise level: {10 * np.log10(broadband_noise_level):.1f} dB') dat = data[idx]
# plt.axhline(y=10 * np.log10(noise_floor), color='blue', linewidth=2, linestyle='--', label=f'Broadband noise floor: {10 * np.log10(noise_floor):.1f} dB') psd = Sxx[idx]
plt.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles time = dat[:, 0]
for freq, power in zip(peak_frequencies, peak_powers): dt = time[1] - time[0] # Sampling interval / time step
plt.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency N = time.shape[0] # Total number of data points
f_s = 1 / dt # Sampling frequency
fNQ = f_s / 2 # Determine Nyquist frequency
faxis = smoothen_data(np.linspace(0,fNQ,N//2), window_size) # Construct frequency axis
""" Noise peak finder"""
if peak_find:
freq_range = (50, max(faxis))
threshold = 10**(peak_find_threshold/10)
peak_powers, peak_frequencies = find_noise_peaks(psd, faxis, freq_range, threshold)
""" Noise levels in units of V/sqrt(Hz)"""
idx_ini = closest(faxis, 100)
idx_fin = closest(faxis, 1000)
noiseat1hunHz = np.sqrt(psd[idx_ini])
noiseat1kHz = np.sqrt(psd[idx_fin])
cum_noise = np.sqrt(integrate.cumtrapz(psd[idx_ini:idx_fin], initial=0))
print('Noise level @ 100 Hz = {}'.format(noiseat1hunHz))
print('Noise level @ 1 kHz = {}'.format(noiseat1kHz))
print('Integrated Noise from 100 Hz to 1 kHz = {}'.format(cum_noise[-1]))
# Plot Power Spectrum in dB
plt.semilogx(faxis, 10 * np.log10(psd), linewidth=2, label = data_str[idx])
if peak_find:
plt.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles
for freq, power in zip(peak_frequencies, peak_powers):
plt.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency
plt.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid plt.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
plt.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid plt.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
@ -262,15 +268,63 @@ def plot_analysis(data, data_bkg, Sxx, Sxx_bkg, data_str, bkg_str, peak_find_thr
for val in x_multiples_of_10: for val in x_multiples_of_10:
plt.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10 plt.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10
f_sig_idx = np.argmax(Sxx)
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / np.sum(np.delete(Sxx, f_sig_idx)))
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / noise_floor)
plt.xlim([min(faxis), max(faxis)]) plt.xlim([min(faxis), max(faxis)])
# plt.ylim([-100, 10]) # plt.ylim([-150, -80])
plt.legend(loc = 3, fontsize=12) plt.legend(loc = 3, fontsize=12)
plt.xlabel('Frequency [Hz]', fontsize=14) plt.xlabel('Frequency [Hz]', fontsize=14)
plt.ylabel('Power Spectral Density [dBV/Hz]', fontsize=14) plt.ylabel('PSD [dB/Hz]', fontsize=14)
# plt.title('SNR= %.2f dB' % (SNR_f), fontsize=14)
# Adjust layout
plt.tight_layout()
else:
time = data[:, 0]
dt = time[1] - time[0] # Sampling interval / time step
N = time.shape[0] # Total number of data points
f_s = 1 / dt # Sampling frequency
fNQ = f_s / 2 # Determine Nyquist frequency
faxis = smoothen_data(np.linspace(0,fNQ,N//2), window_size) # Construct frequency axis
""" Noise peak finder"""
if peak_find:
freq_range = (50, max(faxis))
threshold = 10**(peak_find_threshold/10)
peak_powers, peak_frequencies = find_noise_peaks(Sxx, faxis, freq_range, threshold)
""" Noise levels in units of V/sqrt(Hz)"""
idx_ini = closest(faxis, 100)
idx_fin = closest(faxis, 1000)
noiseat1hunHz = np.sqrt(Sxx[idx_ini])
noiseat1kHz = np.sqrt(Sxx[idx_fin])
cum_noise = np.sqrt(integrate.cumtrapz(Sxx[idx_ini:idx_fin], initial=0))
print('Noise level @ 100 Hz = {}'.format(noiseat1hunHz))
print('Noise level @ 1 kHz = {}'.format(noiseat1kHz))
print('Integrated Noise from 100 Hz to 1 kHz = {}'.format(cum_noise[-1]))
plt.figure(figsize=(12, 8))
# Plot Power Spectrum in dB
plt.semilogx(faxis, 10 * np.log10(Sxx), color='green', linewidth=2, label = data_str)
if peak_find:
plt.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles
for freq, power in zip(peak_frequencies, peak_powers):
plt.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency
plt.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
plt.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
# Calculate the x-axis values for multiples of 10
x_multiples_of_10 = [10**i for i in range(int(np.log10(min(faxis[faxis > 0]))), int(np.log10(max(faxis[faxis > 0]))) + 1)]
# Add thick lines for multiples of 10
for val in x_multiples_of_10:
plt.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10
plt.xlim([min(faxis), max(faxis)])
# plt.ylim([-50, 0])
plt.legend(loc = 3, fontsize=12)
plt.xlabel('Frequency [Hz]', fontsize=14)
plt.ylabel('PSD [dB/Hz]', fontsize=14)
# plt.title('SNR= %.2f dB' % (SNR_f), fontsize=14) # plt.title('SNR= %.2f dB' % (SNR_f), fontsize=14)
# Adjust layout # Adjust layout
@ -279,55 +333,143 @@ def plot_analysis(data, data_bkg, Sxx, Sxx_bkg, data_str, bkg_str, peak_find_thr
# Show plot # Show plot
plt.show() plt.show()
else: def plot_rin(data, rin, data_str, peak_find = False, peak_find_threshold = -100, window_size = 21):
# Create subplots
if isinstance(data, list) and isinstance(rin, list) and isinstance(data_str, list):
plt.figure(figsize=(12, 8)) plt.figure(figsize=(12, 8))
gs = gridspec.GridSpec(2, 3, width_ratios=[1, 1, 1], height_ratios=[1, 1])
# Plot 1: Time vs Voltage for idx in range(len(data)):
axs1 = plt.subplot(gs[0, 0:])
axs1.plot(time_bkg, voltages_bkg, marker='o', color='orange', linewidth=0.5, ms=1, label = bkg_str)
axs1.plot(time, voltages, marker='o', color='green', linewidth=0.5, ms=1, label = data_str)
axs1.set_ylim([-0.5, 0.5])
axs1.set_xlabel('Time (s)', fontsize=14)
axs1.set_ylabel('Voltage (V)', fontsize=14)
axs1.legend(loc = 1, fontsize=12)
axs1.autoscale(tight=True)
axs1.grid(True)
# Plot 2: Power Spectrum in dB dat = data[idx]
axs2 = plt.subplot(gs[1, 0:]) rn = rin[idx]
axs2.semilogx(faxis, 10 * np.log10(Sxx_bkg), color='orange', linewidth=0.5, label = bkg_str)
axs2.semilogx(faxis, 10 * np.log10(Sxx), color='green', linewidth=2, label = data_str)
# axs2.axhline(y=10 * np.log10(broadband_noise_level), color='red', linewidth=2, linestyle='--', label=f'Broadband cumulative noise level: {10 * np.log10(broadband_noise_level):.1f} dB') time = dat[:, 0]
# axs2.axhline(y=10 * np.log10(noise_floor), color='blue', linewidth=2, linestyle='--', label=f'Broadband noise floor: {10 * np.log10(noise_floor):.1f} dB') dt = time[1] - time[0] # Sampling interval / time step
N = time.shape[0] # Total number of data points
f_s = 1 / dt # Sampling frequency
fNQ = f_s / 2 # Determine Nyquist frequency
faxis = smoothen_data(np.linspace(0,fNQ,N//2), window_size) # Construct frequency axis
axs2.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles """ Noise peak finder"""
for freq, power in zip(peak_frequencies, peak_powers): if peak_find:
axs2.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency freq_range = (50, max(faxis))
threshold = 10**(peak_find_threshold/10)
peak_powers, peak_frequencies = find_noise_peaks(rn, faxis, freq_range, threshold)
# Plot RIN
plt.semilogx(faxis, rn, linewidth=2, label = data_str[idx])
if peak_find:
plt.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles
for freq, power in zip(peak_frequencies, peak_powers):
plt.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency
plt.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
plt.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
axs2.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
axs2.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
# Calculate the x-axis values for multiples of 10 # Calculate the x-axis values for multiples of 10
x_multiples_of_10 = [10**i for i in range(int(np.log10(min(faxis[faxis > 0]))), int(np.log10(max(faxis[faxis > 0]))) + 1)] x_multiples_of_10 = [10**i for i in range(int(np.log10(min(faxis[faxis > 0]))), int(np.log10(max(faxis[faxis > 0]))) + 1)]
# Add thick lines for multiples of 10 # Add thick lines for multiples of 10
for val in x_multiples_of_10: for val in x_multiples_of_10:
axs2.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10 plt.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10
f_sig_idx = np.argmax(Sxx) plt.xlim([min(faxis), max(faxis)])
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / np.sum(np.delete(Sxx, f_sig_idx))) # plt.ylim([-150, -80])
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / noise_floor) plt.legend(loc = 3, fontsize=12)
plt.xlabel('Frequency [Hz]', fontsize=14)
plt.ylabel('RIN [dB/Hz]', fontsize=14)
# plt.title('SNR= %.2f dB' % (SNR_f), fontsize=14)
axs2.set_xlim([min(faxis), max(faxis)]) # Adjust layout
# axs2.set_ylim([-100, 10]) plt.tight_layout()
axs2.legend(loc = 3, fontsize=12) else:
axs2.set_xlabel('Frequency [Hz]', fontsize=14) time = data[:, 0]
axs2.set_ylabel('Power Spectral Density [dBV/Hz]', fontsize=14)
# axs2.set_title('SNR= %.2f dB' % (SNR_f), fontsize=14) dt = time[1] - time[0] # Sampling interval / time step
N = time.shape[0] # Total number of data points
f_s = 1 / dt # Sampling frequency
fNQ = f_s / 2 # Determine Nyquist frequency
faxis = smoothen_data(np.linspace(0,fNQ,N//2), window_size) # Construct frequency axis
""" Noise peak finder"""
if peak_find:
freq_range = (50, max(faxis))
threshold = 10**(peak_find_threshold/10)
peak_powers, peak_frequencies = find_noise_peaks(rin, faxis, freq_range, threshold)
plt.figure(figsize=(12, 8))
# Plot Power Spectrum in dB
plt.semilogx(faxis, rin, color='green', linewidth=2, label = data_str)
if peak_find:
plt.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles
for freq, power in zip(peak_frequencies, peak_powers):
plt.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency
plt.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
plt.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
# Calculate the x-axis values for multiples of 10
x_multiples_of_10 = [10**i for i in range(int(np.log10(min(faxis[faxis > 0]))), int(np.log10(max(faxis[faxis > 0]))) + 1)]
# Add thick lines for multiples of 10
for val in x_multiples_of_10:
plt.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10
plt.xlim([min(faxis), max(faxis)])
# plt.ylim([-50, 0])
plt.legend(loc = 3, fontsize=12)
plt.xlabel('Frequency [Hz]', fontsize=14)
plt.ylabel('RIN [dB/Hz]', fontsize=14)
# plt.title('SNR= %.2f dB' % (SNR_f), fontsize=14)
# Adjust layout # Adjust layout
plt.tight_layout() plt.tight_layout()
# Show plot # Show plot
plt.show() plt.show()
def plot_e_folding_time(data, Skk, fstart, fend, data_str, window_size = 21):
time, voltages = data[:, 0], data[:, 1]
dt = time[1] - time[0] # Define the sampling interval
N = voltages.shape[0] # Define the total number of data points
T = N * dt # Define the total duration of the data
df = 1 / T.max()
fNQ = 1 / dt / 2 # Determine Nyquist frequency
faxis = smoothen_data(np.linspace(0,fNQ,N//2), window_size) # Construct frequency axis
idx_ini = closest(faxis, 2*fstart)
idx_fin = closest(faxis, 2*fend)
Skk = Skk[idx_ini:idx_fin]
freqs = np.linspace(fstart, fend, len(Skk))
e_folding_time = [1/(np.pi**2 * freqs[idx]**2 * Skk[idx]) for idx in range(len(freqs))]
plt.figure(figsize=(12, 8))
# Plot Power Spectrum in dB
plt.semilogx(freqs, e_folding_time, color='red', linewidth=2, alpha=0.75, label = data_str)
plt.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
plt.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
# Calculate the x-axis values for multiples of 10
x_multiples_of_10 = [10**i for i in range(int(np.log10(min(freqs[freqs > 0]))), int(np.log10(max(freqs[freqs > 0]))) + 1)]
# Add thick lines for multiples of 10
for val in x_multiples_of_10:
plt.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10
# plt.xlim([min(freqs), max(freqs)])
# plt.ylim([-50, 0])
plt.legend(loc = 3, fontsize=12)
plt.xlabel('Frequency [Hz]', fontsize=14)
plt.ylabel('$T_e$ [s]', fontsize=14)
plt.title('Lower Bound= %.5f s' % (min(e_folding_time)), fontsize=14)
# Adjust layout
plt.tight_layout()
# Show plot
plt.show()

506
Time-Series-Analyzer/Intensity Noise Analysis.ipynb
View File

File diff suppressed because one or more lines are too long

117
Time-Series-Analyzer/computeRIN.m Normal file
View File

@ -0,0 +1,117 @@
% 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'); %
%------------------------------------------%