Skip to content
Snippets Groups Projects
Commit 19d89ac8 authored by osherso2's avatar osherso2
Browse files

Computes the noise pre and post aliasing of a generic TDM readout system

parent 012e5a23
No related branches found
No related tags found
1 merge request!1S4 design
Pipeline #184838 canceled
import numpy
import matplotlib.pyplot as plt
font = {'size' : 12}
from matplotlib import rc
rc('font', **font)
hbar = 6.62607015e-34 # J*s
kb = 1.3806503e-23 # J/K
dpi = 2.*numpy.pi
def butterworthed(NEI, cut_off=60., n_poles=4):
return lambda f: NEI(f)/(1. + (f/cut_off)**(2.*n_poles))
def alias(NEI, new_fs, old_fs):
assert(old_fs > new_fs)
nyq_old = old_fs/2.
nyq_new = new_fs/2.
def aliased(f):
N_set = numpy.arange(numpy.ceil((-nyq_old-f)/new_fs), numpy.floor((nyq_old-f)/new_fs)+1)
f_set = numpy.abs(f + N_set*new_fs)
return (f <= nyq_new)*numpy.sum(NEI(f_set)**2.)**.5
return numpy.vectorize(aliased)
def band_rms(NEI, band=[1e-3, 30.], N=10000):
band_min = band[0]; band_max = band[1]
assert(band_min < band_max)
f_arr = numpy.linspace(band_min, band_max, N)
return numpy.trapz(numpy.nan_to_num(NEI(f_arr))**2., f_arr)**.5
class noise_model:
# Need a few extra parameters to model readout
def __init__(self, det, readout_plateau, pink_knee, roll_off, m, dwell_by_tau):
self.det = det
self.pink_knee = pink_knee
self.roll_off = roll_off
self.plateau = readout_plateau
self.dwell_by_tau = dwell_by_tau
self.m = m
def get_Fp(self):
det = self.det
T = det.T
Tb = det.Tb
betaG = det.betaG
b1 = betaG + 1.; b2 = 2.*betaG + 3.
Fp = (b1/b2) * (1. - (Tb/T)**b2)/(1. - (Tb/T)**b1)
return Fp
### Responsivity stuff ###
def responsivity_matrix(self):
det = self.det
R = det.R
I = det.I
L = det.L
Rsh = det.Rsh
G = det.get_G()
C = det.get_C()
beta = det.get_beta()
loop_gain = det.get_loop_gain()
t0 = C/G; tel = L/R
M11 = lambda w: G*(1. - loop_gain + 1.j*w*t0)
M12 = lambda w: -I*R*(2. + beta)*(w/w)
M21 = lambda w: G*loop_gain*(w/w)
M22 = lambda w: I*(Rsh + R*(1. + beta))*(1. + 1.j*w*tel)
M = (M11, M12, M21, M22)
prefactor = lambda w: 1./(M11(w)*M22(w) - M12(w)*M21(w))
S11 = lambda w: prefactor(w)*M22(w)
S12 = lambda w: -prefactor(w)*M12(w)
S21 = lambda w: -prefactor(w)*M21(w)
S22 = lambda w: prefactor(w)*M11(w)
S = (S11, S12, S21, S22)
return M, S
def M_response(self, dI, dT):
M, _ = self.responsivity_matrix()
M11, M12, M21, M22 = M
dPj = lambda w: M21(w)*dT(w) + M22(w)*dI(w)
dPq = lambda w: M11(w)*dT(w) + M12(w)*dI(w)
return dPj, dPq
def S_response(self, dPj, dPq):
_, S = self.responsivity_matrix()
S11, S12, S21, S22 = S
dI = lambda w: S21(w)*dPq(w) + S22(w)*dPj(w)
dT = lambda w: S11(w)*dPq(w) + S12(w)*dPj(w)
return dI, dT
def NEI_to_NEP(self, dI):
_, S = self.responsivity_matrix()
_, _, S21, _ = S
return lambda w: dI(w)/S21(w)
def NEI_to_absNEP(self, dI):
_, S = self.responsivity_matrix()
_, _, S21, _ = S
return lambda w: numpy.absolute(dI(w)/S21(w))
def dP_to_NEP(self, dPj, dPq):
_, S = self.responsivity_matrix()
S11, S12, S21, S22 = S
dP = lambda w: dPq(w) + S22(w)*dPj(w)/S21(w)
return dP
def dP_to_absNEP(self, dPj, dPq):
_, S = self.responsivity_matrix()
S11, S12, S21, S22 = S
dP = lambda w: numpy.absolute(dPq(w) + S22(w)*dPj(w)/S21(w))
return dP
### Noise sources and models ###
def readout_noise(self):
# Readout noise is NOT an input to which the detector responds. It is `downstream'.
I = self.det.I
plateau = self.plateau
pink_knee = self.pink_knee*dpi
roll_off = self.roll_off*dpi
pink_coeff = plateau * pink_knee
pink_noise = lambda w: pink_coeff / (1e-3 + w)
low_pass = lambda w: 1. / numpy.sqrt((1. + (w / roll_off)**2.))
dI = lambda w: (pink_noise(w) + plateau) * low_pass(w)
return dI
def photon_noise(self):
det = self.det
Popt = det.Popt
center_freq = det.center_freq
bandwidth = det.bandwidth
shot_noise = (4. * numpy.pi * hbar * center_freq * Popt)
p = lambda w: shot_noise + (Popt**2. / (2. * bandwidth))*(w/w)
dPj = lambda w: 0.*w
dPq = lambda w: numpy.sqrt(p(w))
return dPj, dPq
def phonon_noise(self):
det = self.det
T = det.T
G = det.get_G()
Fp = self.get_Fp()
dPj = lambda w: 0.*w
dPq = lambda w: numpy.sqrt(4. * kb * (T**2.) * G * Fp * (w/w))
return dPj, dPq
def johnson_sh_noise(self):
det = self.det
num_rows = det.num_rows
Tb = det.Tb
Rb = det.Rb
Rsh = det.Rsh
R = det.R
I = det.I
num_rows = det.num_rows
johnson_sh = lambda w: 4. * kb * Tb * Rsh * (w/w)
dPj = lambda w: I*johnson_sh(w)**.5
dPq = lambda w: 0.*w
return dPj, dPq
def johnson_det_noise(self):
det = self.det
R = det.R
I = det.I
T = det.T
beta = det.get_beta()
johnson = lambda w: 4. * kb * T * R * (1. + 2.*beta) * (w/w)
dPj = lambda w: numpy.sqrt(johnson(w))*I
dPq = lambda w: -numpy.sqrt(johnson(w))*I
return dPj, dPq
def excess_noise(self):
m = self.m
dPj_n, dPq_n = self.johnson_det_noise()
dPj = lambda w: (m**2.)*dPj_n(w)
dPq = lambda w: (m**2.)*dPq_n(w)
return dPj, dPq
def total_NEP(self, plot=False):
def quad_add(*args):
def func(w):
res = sum(numpy.absolute(a(w))**2. for a in args[:])
return res**.5
return func
dP_a = self.NEI_to_NEP(self.readout_noise())
dP_b = self.dP_to_NEP(*self.photon_noise())
dP_c = self.dP_to_NEP(*self.phonon_noise())
dP_d = self.dP_to_NEP(*self.johnson_sh_noise())
dP_e = self.dP_to_NEP(*self.johnson_det_noise())
dP_f = self.dP_to_NEP(*self.excess_noise())
dP = quad_add(dP_a, dP_b, dP_c, dP_d, dP_e, dP_f)
if plot:
max_f = self.roll_off * 2.
f_arr = numpy.logspace(0, numpy.log10(dpi*self.roll_off), 10000)
plt.plot(f_arr, 1e18*numpy.absolute(dP_a(f_arr*dpi)), label="readout")
plt.plot(f_arr, 1e18*numpy.absolute(dP_b(f_arr*dpi)), label="photon")
plt.plot(f_arr, 1e18*numpy.absolute(dP_c(f_arr*dpi)), label="phonon")
plt.plot(f_arr, 1e18*numpy.absolute(dP_d(f_arr*dpi)), label="johnson_sh")
plt.plot(f_arr, 1e18*numpy.absolute(dP_e(f_arr*dpi)), label="johnson_det")
plt.plot(f_arr, 1e18*numpy.absolute(dP_f(f_arr*dpi)), label="excess")
plt.plot(f_arr, 1e18*dP(f_arr*dpi), label="total")
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r"$\nu$ [Hz]")
plt.ylabel(r"$|\delta P|$ [aW/sqrt(Hz)]")
plt.xlim((1, max_f))
plt.legend(ncol=2, prop={'size': 8})
return dP
def total_NEI(self, plot=False):
def quad_add(*args):
def func(w):
res = sum(numpy.absolute(a(w))**2. for a in args[:])
return res**.5
return func
dI_a = self.readout_noise()
dI_b = self.S_response(*self.photon_noise())[0]
dI_c = self.S_response(*self.phonon_noise())[0]
dI_d = self.S_response(*self.johnson_sh_noise())[0]
dI_e = self.S_response(*self.johnson_det_noise())[0]
dI_f = self.S_response(*self.excess_noise())[0]
dI = quad_add(dI_a, dI_b, dI_c, dI_d, dI_e, dI_f)
if plot:
max_f = self.roll_off * 2.
f_arr = numpy.logspace(0, numpy.log10(dpi*self.roll_off), 10000)
plt.plot(f_arr, 1e12*numpy.absolute(dI_a(f_arr*dpi)), label="readout")
plt.plot(f_arr, 1e12*numpy.absolute(dI_b(f_arr*dpi)), label="photon")
plt.plot(f_arr, 1e12*numpy.absolute(dI_c(f_arr*dpi)), label="phonon")
plt.plot(f_arr, 1e12*numpy.absolute(dI_d(f_arr*dpi)), label="johnson_sh")
plt.plot(f_arr, 1e12*numpy.absolute(dI_e(f_arr*dpi)), label="johnson_det")
plt.plot(f_arr, 1e12*numpy.absolute(dI_f(f_arr*dpi)), label="excess")
plt.plot(f_arr, 1e12*dI(f_arr*dpi), label="total")
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r"$\nu$ [Hz]")
plt.ylabel(r"$|\delta I|$ [pA/sqrt(Hz)]")
plt.xlim((1, max_f))
plt.legend(ncol=2, prop={'size': 8})
return dI
### Sampling and aliasing ###
def multiplexer_bandwidth(self):
sq_bw = self.roll_off # SQUID bandwidth in NOT omega
factr = self.dwell_by_tau * self.det.num_rows # number of taus between resampling
return sq_bw/factr
def total_multiplexed_noise(self, science_freq=1., plot_NEI=False, plot_NEP=False, title=""):
multi_f = self.multiplexer_bandwidth()
multi_w = multi_f * dpi
max_f = self.roll_off * 5. # extra factor for extra padding in aliasing
max_w = max_f * dpi
science_omega = science_freq * dpi
if plot_NEI and plot_NEP:
plt.subplot(211)
plt.xlabel("")
tNEI = self.total_NEI(plot=plot_NEI)
rNEI = self.readout_noise()
alias_tNEI = alias(tNEI, multi_w, max_w)
alias_rNEI = alias(rNEI, multi_w, max_w)
if plot_NEI and plot_NEP: plt.subplot(212)
tNEP = self.total_NEP(plot=plot_NEP)
rNEP = self.NEI_to_absNEP(rNEI)
alias_tNEP = self.NEI_to_absNEP(alias_tNEI)
alias_rNEP = self.NEI_to_absNEP(alias_rNEI)
lowf_NEP_tot = alias_tNEP(science_omega)
lowf_NEP_out = alias_rNEP(science_omega)
lowf_NEP_det = numpy.sqrt(lowf_NEP_tot**2. - lowf_NEP_out**2.)
NEP_frac = lowf_NEP_tot/lowf_NEP_det
entry1 = (1e18*lowf_NEP_tot, 100.*NEP_frac, science_freq)
title+= "%.1f aW/sqrt(Hz) total NEP (%.1f%% det NEP) at %.1f Hz\n" % entry1
if plot_NEI and plot_NEP: plt.suptitle(title)
elif plot_NEI or plot_NEP: plt.title(title)
if plot_NEI:
if plot_NEI and plot_NEP: plt.subplot(211)
f_arr = numpy.logspace(-3, numpy.floor(10.*numpy.log10(multi_f))/10., 10000)
plt.plot(f_arr, 1e12*alias_tNEI(dpi*f_arr), lw=3, ls=":", alpha=.5, label="aliased total")
plt.plot(f_arr, 1e12*alias_rNEI(dpi*f_arr), lw=3, ls=":", alpha=.5, label="aliased readout")
plt.axvline(multi_f/2., color=(.2, .2, .2), label="nyquist frequency")
plt.legend(ncol=2, prop={'size': 8})
if plot_NEP:
if plot_NEI and plot_NEP:
plt.subplot(212)
plt.gca().get_legend().remove()
f_arr = numpy.logspace(-3, numpy.floor(10.*numpy.log10(multi_f))/10., 10000)
plt.plot(f_arr, 1e18*alias_tNEP(dpi*f_arr), lw=3, ls=":", alpha=.5, label="aliased total")
plt.plot(f_arr, 1e18*alias_rNEP(dpi*f_arr), lw=3, ls=":", alpha=.5, label="aliased readout")
if not plot_NEI: plt.legend(ncol=2, prop={'size': 8})
plt.axvline(multi_f/2., color=(.2, .2, .2))
return lowf_NEP_tot, lowf_NEP_det, lowf_NEP_out
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment