import numpy as np
import scipy as scp
from pydiana import ROOT
import pydiana.tools.math_tools as mts
from pydiana.tools import graph_tools as gt
import plotly.express as px
import plotly.graph_objects as go
from scipy import constants as const
from array import array
[docs]
def make_NEP(ap,an,ADC2Amp:float=None,fs:float=None,CalibConst:float=None):
"""
Function for calculating the noise equivalent power.
Parameters
----------
ap: tuple with the average pulse in ADC counts
an: noise abs(FFT)^2
ADC2Amp: Conversion from the ADC values to signal amplitude value (if not given defaults to 1)
fs: sampling frequency in Hz (if not given defaults to 1)
CalibConst: calibration constant for passing from Amp-unit to eV. Disabled if ADC2Amp not provided
Returns
-------
Frequencies: in Hz if sampling freq is provided
NEP: in Watts/Hz^0.5 if ADC2Amp,CalibConst and sampling freq are provided
"""
print("CHECK AP NORMALIZATION")
if fs is None:
fs=1
if CalibConst is None or ADC2Amp is None:
CalibConst = 1
else:
CalibConst *=const.physical_constants['electron volt-joule relationship'][0]*1e18 # move to aJ/mV
if ADC2Amp is None:
ADC2Amp=1
#make frequency
freqs=np.arange(len(ap)//2+1)*fs/len(ap)#np.fft.fftfreq(n=len(ap),d=1/fs)
#make AP fft
ap = ap-ap[0]
ap /=np.max(np.abs(ap))
apFFT = np.abs(np.fft.fft(ap,norm='backward'))[:len(freqs)]
apFFT*=2/fs
#apFFT *= 2/(fs*len(ap))#non sono sicuro
#make AN PSD
_,anPSD = make_PSD_from_FFT(an=an,fs=fs,conv = ADC2Amp*CalibConst)#noise is in J^2/Hz
anPSD = anPSD[:len(freqs)]
#Calulate NEP
NEP = np.sqrt(anPSD)/apFFT #nep is in Watt/Hz^0.5
return freqs,NEP
[docs]
def make_NEP_plot(ap,an,ADC2Amp:float=None,fs:float=None,CalibConst:float=None,unitAmp:str=None,showplot:bool=True,log_x:bool=True,log_y:bool=True):
"""
Function for plotting NEP calculated from average pulse and average noise
Parameters
----------
ap: tuple with the average pulse in ADC counts
an: noise abs(FFT)^2
ADC2Amp: Conversion from the ADC values to signal amplitude value (if not given defaults to 1)
fs: sampling frequency in Hz (if not given defaults to 1)
CalibConst: calibration constant for passing from Amp-unit to eV. Disabled if ADC2Amp not provided
unitAmp: string with the amplitude unit of the signal (defaults to mV). It is used only if calib const is not provided.
showplot: boolean for enabling plot showing
log_x: boolean for using log scale on x axis
log_y: boolean for using log scale on y axis
Returns
-------
Figure instance with NEP plot
"""
freq,NEP=make_NEP(ap=ap,an=an,ADC2Amp=ADC2Amp,fs=fs,CalibConst=CalibConst)
unit_freq='idx'
def GetBaselineResolution(NEP,freq):
"""
"""
df = np.diff(freq)[1]#get frequency resolution
variance = 1/(np.sum(1/NEP**2)*df)
sigma=np.sqrt(variance)#in Joule
sigma*=const.physical_constants['joule-electron volt relationship'][0]*1e-18#conversion to aJ
return sigma
print(f"Baseline Resolution: {GetBaselineResolution(NEP,freq):.4g} eV")
if ADC2Amp is None:
unitAmp = 'ADCunits'
if CalibConst is not None and ADC2Amp is not None:
unitAmp = 'aW'
if unitAmp is None:
unitAmp = 'mV'
if fs is not None:
unit_freq='Hz'
ylabel = f"{unitAmp} {unit_freq}/√{unit_freq}"
if unit_freq=='Hz' and unitAmp=='aW':
ylabel = f"{unitAmp}/√{unit_freq}"
elif unitAmp=='aW':
ylabel = f"aJ {unit_freq}/√{unit_freq}"
ylabel="NEP ["+ylabel+']'
xlabel = "Frequency"
if unit_freq=='Hz':
xlabel +=' [Hz]'
else:
xlabel +=' Index [idx]'
if log_x:
freq = freq[1:]
NEP = NEP[1:]
fig = go.Figure()
fig.add_trace(go.Scatter(x=freq,y=NEP,mode='lines'))
fig.update_layout(title='Noise Equivalent Power')
fig.update_xaxes(title=xlabel,type='log' if log_x else 'linear')
fig.update_yaxes(title=ylabel,type='log' if log_y else 'linear')
if showplot:
fig.show()
return fig
[docs]
def make_PSD_from_FFT(an,fs:float=1,conv:float=1):
"""
Function for calculating and plotting Power Spectral Density from the absolute values of the FFT
Parameters
----------
an: array with the abs(fft)
fs: sampling frequency (defaults to 1).
conv: normalization conversion for the amplitude of the FFT (defaults to 1)
Returns
-------
x:array with frequency bins
y:array with frequency bin height
"""
nbin = len(an)//2+1
res = fs/len(an)
x = np.arange(nbin)*res
y= 2*np.array(an)*((conv/len(an))**2)/res
return x,y
[docs]
def make_PSD_plot(an,fs:float=None,conv:float=None,unitf:str=None,unitA:str=None,usesqrt:bool=False,showplot:bool=False):
"""
Function for calculating and plotting Power Spectral Density from the absolute values of the FFT
Parameters
----------
an: array with the abs(fft)
fs: sampling frequency (defaults to 1).
conv: normalization conversion for the amplitude of the FFT (defaults to 1)
unitf: string with the unit for the frequency axis (defaults to Hz)
unitA: string with the unit of the the abs(FFT) values (defaults to mV)
showplot: bool for enabling the plotting of the figure
Returns
-------
Plotly figure with PSD
"""
x,y=make_PSD_from_FFT(an=an,
fs=fs if fs is not None else 1,
conv=conv if conv is not None else 1,
)
if fs is None:
unitf = 'Frequency Bins'
elif unitf is None:
unitf = 'Hz'
if conv is None:
unitA = 'ADC'
elif unitA is None:
unitA = 'mV'
unitf2=unitf
if not usesqrt:
unitA=f"{unitA}<sup>2</sup>"
else:
unitf2=f"√{unitf}"
y=np.sqrt(y)
fig = gt.make_scatter(x=x[1:],y=y[1:],
xlab=f'Frequency [{unitf}]',
ylab=f'Power [{unitA}/{unitf2}]',
mode='l',logx=True,logy=True,showplot=showplot,
)
fig.update_layout(title='Noise Power Spectral Density')
fig.update_xaxes(autorange=True)
fig.update_yaxes(autorange=True)
return fig
[docs]
def make_periodogram(traces,window='cosine',fs=1,ADC2unit=1):
"""
Function for making the periodogram of traces with windowing. It will also convert the trace from ADC to a specified unit.
Parameters
----------
traces: 2D array with the traces
window: see https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window
fs: sampling frequency
ADC2unit: float with the conversion value from ADC units to a wanted unit like mV
Returns
-------
2D array with the periodogram
1D array with the corresponding frequencies
"""
traces = np.array(traces)
f,pxx = scp.signal.welch(noise_traces*ADC2unit,
fs=fs,
window=window,
nperseg=traces.shape[-1]-1,
scaling='spectrum',
)
return pxx,f
[docs]
def make_PSD(traces,window='cosine',fs=1,ADC2unit=1):
"""
Function for making the power spectral density of traces with windowing. It will also convert the trace from ADC to a specified unit.
Parameters
----------
traces: 2D array with the traces
window: see https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window
fs: sampling frequency
ADC2unit: float with the conversion value from ADC units to a wanted unit like mV
Returns
-------
2D array with the periodogram
1D array with the corresponding frequencies
"""
traces = np.array(traces)
f,pxx = scp.signal.welch(traces*ADC2unit,
fs=fs,
window=window,
nperseg=traces.shape[-1]-1,
scaling='density',
)
return pxx,f
[docs]
def make_waterfall(traces,
index:list=None,
window='cosine',
fs:float=None,
ADC2unit:float=None,
useSpacing:bool=False,
logx:bool=False,
logy:bool=False,
logz:bool=False,
xlab:str=None,
ylab:str=None,
unit:str='mV',
fillwith:str='nan',
colorscale:str='Turbo',
zrange=None,
showplot=True):
"""
Function for making the waterfall of the PSD of provided traces
Parameters
----------
traces : 2D array with the traces
index : 1D array with the quantity with respect to which order the traces. If None the provided traces are not ordered
window : specifying windowing type for PSD. See make_window function
fs : sampling frequency
ADC2unit : conversion quantity from the unit of the traces to a specified unit
useSpacing : Switch for enabling correct spacing between the various PSD in the final waterfall. The spacing is based on the index array
logx : switch for enabling logscale on x axis
logy : switch for enabling logscale on y axis
logz : switch for enabling logscale on color-axis
xlab : label for the xaxis
ylab : label for the yaxis
unit : label of the unit of the traces after being multiplied with ADC2unit
fillwith : if useSpacing is enabled specifies how to fill the space between PSDs on the waterfall. If 'min' a value smaller of the minimum in the waterfall is used, if 'nan' the spaces are filled with NaNs
colorscale : string specifying the colorscale
zrange: if not None the PSD is saturated inside range
showplot : switch for enabling the showing of the plot
Returns
-------
2D array with the waterfall
1D array with the frequencies
1D array with the indexes (which if useSpacing is enabled contained the 'missing' indexes of the empty spaces)
go.Figure instance with the plot of the produced waterfall
"""
traces= np.array(traces)
hovertemplate=''
if (index is None or not useSpacing) and xlab is None:
xlab='Index'
hovertemplate ='<b>idx</b>:%{x}'
elif xlab is not None:
hovertemplate =f'<b>{xlab[:4]}</b>:'+'%{x}'
if index is not None:
traces,index = mts.sortmultivec(traces,index,lead=1)
traces=np.array(traces)
if index is None or not useSpacing:
index = np.arange(traces.shape[0])
hovertemplate ='<b>idx</b>:%{x}'
hovertemplate+='<br>'
waterfall,freqs = make_PSD(traces=traces,
window=window,
fs=fs if fs is not None else 1,
ADC2unit=ADC2unit if ADC2unit is not None else 1)
index2plot=index
fillfactor = np.nan
if fillwith == 'min':
fillfactor =np.min(waterfall)*(0.1 if logz else 0.5)
if useSpacing:
step = np.min(np.diff(index))
waterfalltmp=[]
oldstep=-1
index2plot=[]
for i,psd in enumerate(waterfall):
if i==0:
waterfalltmp.append(psd)
oldstep =index[i]
index2plot.append(index[i])
else:
steps_elapsed= int((index[i]-oldstep)/step)
if steps_elapsed>=2:
for j in range(steps_elapsed-1):
waterfalltmp.append(fillfactor*np.ones(psd.shape[-1]))
index2plot.append(j*step + oldstep)
oldstep=index[i]
waterfalltmp.append(psd)
index2plot.append(index[i])
waterfall=waterfalltmp
if ylab is None:
hovertemplate+='<b>Freq</b>:'
ylab = 'Frequency '
if fs is None:
hovertemplate+='%{y} idx'
ylab+='[Frequency Idx]'
else:
hovertemplate+='%{y:.2f} Hz'
ylab+='[Hz]'
hovertemplate+='<br><b>PSD</b>:10^%{z:.2f}'
zlab = 'PSD<br>'
if ADC2unit is None:
zlab+='[ADC^<sup>2</sup>/'
else:
zlab+=f'[{unit}<sup>2</sup>/'
if fs is None:
zlab +='Frequency Idx]'
else:
zlab +='Hz]'
toplot=np.transpose(waterfall)
freqs2plot=freqs
zmin=np.min(toplot)
zmax=np.max(toplot)
if zrange is not None:
zmax=max(zrange)
zmin=min(zrange)
toplot = np.where(toplot<zmax,toplot,zmax)
toplot = np.where(toplot>=zmin,toplot,zmin)
if logz:
toplot = toplot[1:]#Remove DC
freqs2plot=freqs2plot[1:]
color_vals,color_names = gt.makelogaxis(np.min(toplot),np.max(toplot))
toplot=np.log10(toplot)
zmin = min(color_vals)
zmax = max(color_vals)
fig=go.Figure()
fig.add_trace(go.Heatmap(y=freqs2plot,x=index2plot,z=toplot,
colorscale=colorscale,
colorbar_title={'text':zlab},
zmin=zmin,
zmax=zmax,
hovertemplate=hovertemplate,
showlegend=False,
name='',
))
fig.update_xaxes(type='log' if logx else 'linear',title=xlab)
fig.update_yaxes(type='log' if logy else 'linear',title=ylab)
if not useSpacing:
tickvals = np.linspace(min(index),max(index),10)
if logx:
tickvals=np.logspace(np.log10(tickvals[0]),np.log10(tickvals[-1]),len(tickvals))
fig.update_xaxes(tickvals=tickvals,ticktext=[' ']*len(tickvals))
if logz:
fig['data'][0]['colorbar']['tickvals']=color_vals
fig['data'][0]['colorbar']['ticktext']=color_names
if showplot:
fig.show()
return waterfall,freqs,index2plot,fig