Analysis_Sensitivity.py

'''
authored collaboratively by the authors of PUBLICATION_LINK

            #################################################
            ### Local Oscillator and Coupler Optimization ###
            #################################################

This script walks through using a calibrated applied electric field to find the electric field sensitivity in V/m/sqrt(Hz)

It is assumed you have access to the following variables:

Oscilloscope data [Volts] in np.arrays 'probe' and 'LIA' with an optical frequency array 'optical_freqs'
Spectrum Analyzer data [dBm] from the waveguide in np.array 'SA_dBms' with frequencies 'SA_freqs'
Applied Local Oscillator frequency [Hz] as float 'LO'
Applied Signal frequency [Hz] as float 'SIG'
'beatnote' [Hz] = np.abs(LO - SIG)

For the photon shot noise it is assumed that you have:
Total optical power [W] as float 'P_opt'
Responsivity of photodiode [A / W] at probe wavelength as float 'R_lambda'
Total transimpedance gain of photodiode [V / A] as 'gain'  
Probe optical frequency [Hz as float 'probe_frequency'

For a given measured Signal power, convert that to an electric field.  That can either be done via a free-space or atomic calibration factor.
By whatever method, that will result in a float Efield_cal_factor [Efield per applied sqrt(Watt)] = [V/m / sqrt(W)]
'''

#################
##### DUMMY ##### 
#################
# I'll delete this before publication, but this is to check for errors when using the variables below
probe = np.array([])
LIA = np.array([])
optical_freqs = np.array([])
SA_dBms = np.array([])
SA_freqs = np.array([])
LO = 0.; SIG = 0.; beatnote = 0.
freq_shifts = np.array([])
measured_dBms = np.array([])
LIA_peaks = np.array([])
LIA_peaks_freqs = np.array([])
Efield_cal_factor = 0.
t = freq_shifts = np.array([])
P_opt = 0.; R_lambda = 0.; gain = 0.; probe_frequency = 0.
sensitivities = np.array([])
#################
#################
#################

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif' # style choice
import numpy as np
import scipy.constants as cons
import scipy.signal as sig
    
    # varibales of convenience
size = 2**6 # useful to scale image text
figsize = (8,4.9) # usefult to scale image size

################################
### CHOOSE WHAT TO SHOW/SAVE ###
################################
show_plots = True
save_plots = False

###################################
### GET SPECTRUMA ANALYZER PEAK ###
###################################
    # find the Signal peak dBm
wanted_freq_idx = np.argmin([abs(x-SIG) for x in SA_freqs])
SA_dBm = SA_dBms[wanted_freq_idx]

if show_plots:    
    fig,ax=plt.subplots(figsize = figsize)
        # plot data
    ax.plot(SA_freqs,SA_dBms)
    ax.scatter(SA_freqs[wanted_freq_idx], SA_dBm, marker = '*', s = 2**8, color = 'red')   
        # labels
    ax.set_xlabel('Frequency (Hz)', fontsize = size)
    ax.set_ylabel('RF Power (dBm)', fontsize = size)
        # axes
    ax.tick_params(axis = 'both', labelsize = size)
    ax.minorticks_on()            
    ax.tick_params(which = 'major', length = 12)
        # figure layout
    fig.suptitle('Spectrum Analyzer Data', fontsize = size)
    fig.tight_layout()
    
    if save_plots:
        fig.savefig('SpectrumAnalyzer_data.png')

###################
### FIND EFIELD ###
###################
Efield = Efield_cal_factor * np.sqrt(10**((SA_dBm - 30) / 10)) # V/m

##########################################
### FIND PEAK OF POWER SPECTRUM [V**2] ###
##########################################
nperseg = 2**int(np.log2(len(probe) / 50)) # this gives nperseg ~ n_total / 32 in a power of 2 for efficient calculation
V2_freqAxis, V2_spectrum = sig.welch(probe, 1 / np.mean(np.diff(t)), nperseg = nperseg, return_onesided = True, window='hamming', scaling='spectrum') # power spectrum V**2 vs Hz
wanted_beatnote_idx = np.argmin([abs(x-beatnote) for x in V2_freqAxis])
V2_beatnote = np.max(V2_spectrum[max(0,wanted_beatnote_idx-50):min(wanted_beatnote_idx+50, len(V2_spectrum))]) # V**2 at the beatnote

if show_plots:
    fig,ax=plt.subplots(figsize = figsize)
        # plot data
    ax.loglog(V2_freqAxis, V2_spectrum) 
    ax.scatter(V2_freqAxis[wanted_beatnote_idx], V2_beatnote, marker = '*', color = 'r', s = 2**8)
        # labels    
    ax.set_xlabel('Frequency (Hz)', fontsize = size)    
    ax.set_ylabel('$V_{RMS}$', fontsize = size)    
        # axes
    ax.tick_params(axis = 'both', labelsize = size)
    ax.minorticks_on()            
    ax.tick_params(which = 'major', length = 12)
        # figure layout
    fig.suptitle('Probe Power Spectrum', fontsize = size)
    fig.tight_layout()
    
    if save_plots: # save plot
        fig.savefig('Probe_Power_Spectrum.png')

#################################################
### MEASURED PHOTODIODE RMS VOLTAGE TO EFIELD ###
#################################################
# Now we know that an RMS Voltage of sqrt(V2_beatnote) corresponds to an electric field of Efield
V_to_Efield = Efield / np.sqrt(V2_beatnote) # V/m / V

###############################
### FIND V2 SNR AT BEATNOTE ###
###############################
# one could use this as a flag to omit data where the signal is too weak to be detected with sufficient SNR (eg SNR < 2)
Hz_window = 10000 # +/- (arbitrary)
pk_rgn = (np.abs(V2_freqAxis-beatnote) < Hz_window)

below = (V2_freqAxis > (beatnote - 2 * Hz_window)) & (V2_freqAxis < (beatnote - Hz_window))
above = (V2_freqAxis > (beatnote + Hz_window)) & (V2_freqAxis < (beatnote + 2 * Hz_window))

V2_bkgd = np.mean(np.concatenate((V2_spectrum[below],V2_spectrum[above])))

V2_SNR = V2_beatnote / V2_bkgd

################################################################
### CALCULATE PROBE POWER SPECTRAL DENSITY [V**2 / sqrt(Hz)] ###
################################################################
PSD_freqAxis, PSD_spectrum = sig.welch(probe, 1 / np.mean(np.diff(t)), nperseg = nperseg, return_onesided = True, window='hamming') # PSD V**2/sqrt(Hz) vs Hz

if show_plots:
    fig,ax=plt.subplots(figsize = figsize)
        # plot data
    ax.loglog(PSD_freqAxis, PSD_spectrum) 
    ax.scatter(V2_freqAxis[wanted_beatnote_idx], V2_beatnote, marker = '*', color = 'r', s = 2**8)
        # peak region and regions to be integrated
    ax.axvspan(PSD_freqAxis[pk_rgn][0], PSD_freqAxis[pk_rgn][-1], color='red', alpha=0.3)
    ax.axvspan(PSD_freqAxis[below][0], PSD_freqAxis[below][-1], color='green', alpha=0.3)
    ax.axvspan(PSD_freqAxis[above][0], PSD_freqAxis[above][-1], color='green', alpha=0.3)
        # labels    
    ax.set_xlabel('Frequency (Hz)', fontsize = size)    
    ax.set_ylabel('$V_{RMS}/\sqrt{Hz}$', fontsize = size)    
        # axes
    ax.tick_params(axis = 'both', labelsize = size)
    ax.minorticks_on()            
    ax.tick_params(which = 'major', length = 12)
        # figure layout
    fig.suptitle('Probe Power Spectral Density', fontsize = size)
    fig.tight_layout()
    
    if save_plots: # save plot
        fig.savefig('Probe_Power_spectral_Density.png')

###################
### SENSITIVITY ###
###################
# scale PSD to field units 
spectrum_sensitivity = V_to_Efield * np.sqrt(PSD_spectrum) # V/m / sqrt(Hz)
 # integrate noise floor around beatnote to obtain sensitivity
sensitivity = np.mean(np.concatenate([spectrum_sensitivity[below],spectrum_sensitivity[above]]))

if show_plots: # recreates version of Figure 6 from PUBLICATION_LINK
    fig,ax=plt.subplots(figsize = (16,9))
        # plot data
    ax.loglog(PSD_freqAxis, spectrum_sensitivity)        
    wanted_beatnote_idx = np.argmin([abs(x-beatnote) for x in PSD_freqAxis])
    ax.scatter(PSD_freqAxis[wanted_beatnote_idx], PSD_freqAxis[wanted_beatnote_idx], marker='*', s=2**10, color='red') # points out beatnote
        # Photon Shot Noise
    PSN_optical = np.sqrt(2 * cons.h * probe_frequency * P_opt) # PSN in optical Watts [W / sqrt(Hz)]
    PSN_Voltage = PSN_optical * (R_lambda * gain) # PSN in electrical voltage [V / sqrt(Hz)]            
    PSN_Efield = V_to_Efield * PSN_Voltage # PSN in electric field [V/m / sqrt(Hz)]
    ax.axhline(y = PSN_Efield, linewidth = 2**2, color = 'red', linestyle = '--')    
        # labels
    ax.set_xlabel('Frequency (Hz)', fontsize = size)
    ax.set_ylabel(r'Sensitivity (V/m/$\sqrt{Hz}$)', fontsize = size)
        # axes
    ax.tick_params(axis = 'both', labelsize = size)
    ax.minorticks_on()            
    ax.tick_params(which = 'major', length = 12)    
        # figure layout
    fig.suptitle('Electric_Field_Sensitivity', fontsize = size)
    fig.tight_layout()
    
    if save_plots:
        fig.savefig('Electric_Field_Sensitivity.png')

'''
##################        
### SAVE DATA ###
##################
Save the following data to a file:
    SA_dBm
    sensitivity    

One can now perform statistics on many sensitivity measurements

Assume you have np.arrays of:
    Measured applied RF SIG powers 'measured_dBms' populated with each SA_dBm from above
    Measured electric field sensitivites 'sensitivites' populated with each sensitivity from above
'''

###########################
### AVERAGE SENSITIVITY ###
###########################
good_lower_bound = (V2_SNR > 2) # dBm where beatnote is below noise floor
good_upper_bound = (SA_dBms < -15) # example dBm where heterodyne conditions have been broken
good_data = good_lower_bound & good_upper_bound

S_ave = np.average(sensitivities[good_data])
S_std = np.std(sensitivities[good_data])   

if show_plots: # recreates version of Figure 6 from PUBLICATION_LINK
    import matplotlib.patches as patches

    fig,ax=plt.subplots(figsize = (16,9))
        # plot data
    ax.scatter(SA_dBms,sensitivities, s=2**8)
    xmin = np.min(SA_dBms[good_data]); xmax = np.max(SA_dBms[good_data]) # convenience
    ax.hlines(y = S_ave, xmin = xmin, xmax = xmax, color='red', linewidth=2**3, linestyle = '--')    
    rect0 = patches.Rectangle((xmin, (S_ave - S_std)), np.abs(xmin-xmax), 2 * S_std, 
                         linewidth=2, 
                         edgecolor='red', 
                         facecolor='red', 
                         alpha=0.3) # standard deviation of good data
    ax.add_patch(rect0)
    rect1 = patches.Rectangle((xmin, np.min(sensitivities[good_data])*0.75), np.abs(xmin-xmax),
                         np.max(sensitivities[good_data]*1.25) - np.min(sensitivities[good_data])*0.75,
                         linewidth=2, 
                         edgecolor='green', 
                         facecolor='green', 
                         alpha=0.3) # good data
    ax.add_patch(rect1)
        # labels
    ax.set_xlabel('SIG Input Power (dBm)', fontsize = size)
    ax.set_ylabel('Sensitivity ($\mu V/m/\sqrt{Hz}$)', fontsize = size)
        # axes
    ax.tick_params(axis = 'both', labelsize = size)
    ax.minorticks_on()            
    ax.tick_params(which = 'major', length = 12)
        # figure layout
    fig.suptitle('Spectrum Analyzer Data', fontsize = size)
    fig.tight_layout()

    if save_plots:
        fig.savefig('Electric_Field_Sensitivities.png')