#!/usr/bin/env python3
#
'''
.. module:: Bfit_Utils
:synopsis: Algorithms used for baseline optimal filtering
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import os
import numpy
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from tesfdmtools.methods.IQrotate import IQphase,IQphrot
from tesfdmtools.utils.cu import cu
from tesfdmtools.methods.gaussfit import gaussfit
[docs]def optfilt(hdf,channel,freq,indx,base1,base2,noisspec,debug=False,freqcutoff=None,\
rotate=False,bsec=0.05,absolute=False,prlen=None,avpulsepos=None):
'''
compute optimal filter for pulses
Args:
* `hdf` = HDF5 input file object
* `channel` = channel number being processed
* `freq` = frequency number (pixel) being processed
* `indx` = index of selected events, to be processed
* `base1` = baseline level at start of record
* `base2` = baseline level at end of record
* `noisspec` = noise spectrum
Kwargs:
* `debug` = if True, plot various fit parameters [default: False]
* `freqcutoff` = cutoff frequency (Hz) for optimal filtering [default: None]
* `rotate` = rotation angle (radians) for I/Q pulse record
* `bsec` = section of record to take for background (only used for 'rotate' or 'absolute') [default: 0.05]
* `absolute` = use sqrt(I^2+Q^2) signal [default: False]
* `prlen` = length of record to process (None = entire record) [default: None]
* `avpulsepos` = average position of pulse maximum; used for limited `prlen`. [default: None]
Returns:
* `optf` = optimal filter
* `norm` = optimal filter scaling factor
'''
print("Optimal filter, channel=",channel,' freq=',freq)
aphase=0.0
if type(rotate) is bool:
if rotate:
aphase=getphase(hdf,channel,freq,indx,debug=debug)
print("Optimal fit phase rotation: ",aphase)
else:
if type(rotate) is float:
aphase=rotate
print("Optimal fit phase rotation: ",aphase)
sirecord=hdf.channel[channel].freq[freq][indx[0]][:,0] # read first record to get record parameters
recp1=0 # use entire record
recp2=sirecord.size
if ( prlen is not None ) and ( apulsepos is not None ): # limit record length ?
p1=int(apulsepos-100-bsec*sirecord.size)
if p1 > 0:
recp1=p1
p2=recp1+prlen
if p2 < sirecord.size:
recp2=p2
print("Optimal fit record length limited to: ",recp1,recp2)
sirecord=sirecord[recp1:recp2]
bl=int(bsec*sirecord.size) # section for background level
samplerate=float(cu(hdf.channel[channel].freq[freq].attrs['sample_rate'])) # get sample rate
if freqcutoff is None:
slen=sirecord.size//2 # max samplingfrequency
else:
slen=int(float(freqcutoff)/float(samplerate)*sirecord.size) # frequency cutoff
xax=numpy.arange(sirecord.size,dtype=float) # make x-axsis for record
iax=numpy.arange(sirecord.size,dtype=int)
bax=numpy.concatenate((iax[0:bl],iax[-bl:])) # baseline section
fpulses=numpy.zeros((slen-1),dtype=float) # initialize array to store pulse power spectra
for i,irec in enumerate(indx): # go through all record in selection list
bas1=base1[irec]
bas2=base2[irec]
record=hdf.channel[channel].freq[freq][irec]
irecord=record[:,0]
qrecord=record[:,1]
if aphase != 0.0:
irecord,qrecord=IQphrot(irecord,qrecord,aphase)
if record[0,0] < 0: # wrong rotation, rotate by pi
irecord=-record[:,0]
qrecord=-record[:,1]
if absolute:
irecord=numpy.sqrt(irecord.astype(float)**2+qrecord.astype(float)**2)
irecord=irecord[recp1:recp2] # limit record length
if absolute or ( aphase != 0.0 ):
pp=numpy.polyfit(bax,irecord[bax],1)
cirecord=irecord-(pp[0]*xax+pp[1])
else:
cirecord=irecord-(bas1+xax/float(sirecord.size-1)*(bas2-bas1)) # retrieve I-signal and subtract background
fpulses=fpulses+numpy.absolute(numpy.fft.fft(cirecord)[1:slen]) # compute power spectrum
if ( i % 1000 ) == 0:
print("[optfilt] process record: ",irec)
print("number of pulses processed: ",indx.size)
apulse=fpulses/float(indx.size) # compute average power spectrum for all selected pulses
nns=(noisspec**2)[1:slen]
norm=numpy.sum(apulse**2/nns) # average pulse normalization
optf=apulse/nns # optimal filter
if debug:
fax=numpy.arange(1,slen)/float(slen)*samplerate
f, ax = plt.subplots(1)
ax.plot(fax,optf,'b-')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power')
ax.set_title('Optimal filter')
plt.show()
plt.close('all')
return (optf,norm)
[docs]def nfilt(hdf,channel,freq,ptp,ofilter,norm,debug=False,absolute=False,\
rotate=False,freqcutoff=None,prlen=None,threshold=1000):
'''
Perform optimal filter for baseline records
Args:
* `hdf` = HDF5 input file object
* `channel` = channel number being processed
* `freq` = frequency number (pixel) being processed
* `ptp` = delta difference within records
* `ofilter` = the optimal filter
* `norm` = optimal filter scaling factor
Kwargs:
* `debug` = if True, plot various fit parameters [default: False]
* `absolute` = use sqrt(I^2+Q^2) signal [default: False]
* `rotate` = rotation angle (radians) for I/Q pulse record
* `freqcutoff` = cutoff frequency (Hz) for optimal filtering [default: None]
* `prlen` = length of record to process (None = entire record) [default: None]
* `threshold` = treshold for event detection
Returns:
* `ofits` = list with optimal fits for noise records
* `indx` = index of records with the optimal fits
'''
print("Optimal filter, channel=",channel,' freq=',freq)
aphase=0.0
if type(rotate) is bool:
if rotate:
aphase=getphase(hdf,channel,freq,indx,debug=debug)
print("Optimal fit phase rotation: ",aphase)
else:
if type(rotate) is float:
aphase=rotate
print("Optimal fit phase rotation: ",aphase)
sirecord=hdf.channel[channel].freq[freq][0][:,0] # read first record to get record parameters
recp1=0 # use entire record
recp2=sirecord.size
if ( prlen is not None ) and ( apulsepos is not None ): # limit record length ?
p1=int(apulsepos-100-bsec*sirecord.size)
if p1 > 0:
recp1=p1
p2=recp1+prlen
if p2 < sirecord.size:
recp2=p2
print("Optimal fit record length limited to: ",recp1,recp2)
sirecord=sirecord[recp1:recp2]
samplerate=float(cu(hdf.channel[channel].freq[freq].attrs['sample_rate'])) # get sample rate
if freqcutoff is None:
slen=sirecord.size//2 # max samplingfrequency
else:
slen=int(float(freqcutoff)/float(samplerate)*sirecord.size) # frequency cutoff
xax=numpy.arange(sirecord.size,dtype=float) # make x-axsis for record
iax=numpy.arange(sirecord.size,dtype=int)
ptph,ptpbs=numpy.histogram(ptp,range=[0,threshold],bins=100) # histogram of records deltas
ptpb=(ptpbs[0:-1]+ptpbs[1:])/2.0
hmax=ptpb[numpy.argmax(ptph)]*3.0
thr=numpy.min(numpy.array([threshold,hmax])) # threshold for noiserecord
indx,=numpy.where(ptp < thr)
print("Used baseline-noise threshold: ",thr," (threshold,hmax: ",threshold,hmax,")")
ofits=numpy.zeros(indx.size) # buffer for optimal fits
for i,irec in enumerate(indx): # go through all records in selection list
record=hdf.channel[channel].freq[freq][irec]
irecord=record[:,0]
qrecord=record[:,1]
if aphase != 0.0:
irecord,qrecord=IQphrot(irecord,qrecord,aphase)
if record[0,0] < 0: # wrong rotation, rotate by pi
irecord=-record[:,0]
qrecord=-record[:,1]
if absolute:
irecord=numpy.sqrt(irecord.astype(float)**2+qrecord.astype(float)**2)
irecord=irecord[recp1:recp2] # limit record length
blp=numpy.polyfit(xax,irecord,1)
irecord=irecord-(blp[0]*xax+blp[1]) # subtract baseline
pwsp=numpy.absolute(numpy.fft.fft(irecord)[1:slen]) # compute power spectrum
ofits[i]=numpy.sum(pwsp*ofilter)/norm # optimal filter
return (ofits,indx)
[docs]def fitgauss(pltdata,iplt,hdf,channel,freq,spectrum,bins,debug=False):
'''
Fit gaussian function to the optimal fit spectrum
Args:
* `pltdata` = dictionary for the fit results
* `iplt` = number of the plot to process
* `hdf` = HDF5 input file object
* `channel` = channel number being processed
* `freq` = frequency number (pixel) being processed
* `spectrum` = spectrum of the pulse fit parameters
* `bins` = centers of the spectrum bins
Kwargs:
* `debug` = if True, plot debug information [default: False]
Returns:
* `afit` = average fitted position
* `fitpars` = fit parameters
* fitpars[0] = Gaussian scale
* fitpars[1] = Gaussian position
* fitpars[2] = Gaussian position width
* fitpars[3-5] = Errors for fitpars[0-2]
'''
fffx,=numpy.where( hdf.channel[channel].freq_conf['pixel_index'] == freq )
freqindx=fffx[0]
khz=hdf.channel[channel].freq_conf['freq'][freqindx]
fid="ch=%d, px=%d %7.2f (kHz)" % (channel,freq,khz) # plot title
pars,perr,fit=gaussfit(bins.astype(float),spectrum.astype(float))
afit=pars[1] # average fit in (eV)
yrr=numpy.sqrt(spectrum) # 1-sigma errors
yerr=numpy.where(yrr <= 0.0,1.0,yrr) # for zero counts insert one-count error
fwhmtxt="%7.7s %12.4g +/- %12.4g (eV)" % ("Av.Fit: ",afit,pars[2]) # average fit + error
txt=' '.join(fwhmtxt.split())+'\n'
txt=txt+( "%d counts\n" % (numpy.sum(spectrum)) )
print("Av.Fit: %12.4g +/- %10.4g (eV)" % (afit,pars[2])) # average fit with 1-sigma error
pltdata[iplt]={'xax':bins,'spectrum':spectrum,'yerr':yerr,'gauss':fit,\
'ptitle':fid,'ptxt':txt} # spectral fit data
fitpars=numpy.concatenate((pars,perr))
return (afit,fitpars)