Update to change wavelength_calibartion into a library
This commit is contained in:
parent
e30cd203a0
commit
f3e62bfca3
1 changed files with 55 additions and 109 deletions
|
|
@ -1,30 +1,10 @@
|
||||||
import numpy as np
|
import numpy as np
|
||||||
import scipy.signal as sig
|
import scipy.signal as sig
|
||||||
import matplotlib.backends.backend_qtagg as qt
|
|
||||||
import matplotlib.figure as fig
|
|
||||||
import utils
|
|
||||||
|
|
||||||
calib_top = np.load( 'asset/calib_top.npy' )
|
|
||||||
calib_down = np.load( 'asset/calib_down.npy' )
|
|
||||||
|
|
||||||
peaks_reference_sorted = np.loadtxt( 'calibration/OHP_Hg.calib' )
|
|
||||||
peaks_reference = np.sort( peaks_reference_sorted )
|
|
||||||
|
|
||||||
mean_calib_top = np.mean( calib_top , axis = 0 )
|
|
||||||
mean_calib_down = np.mean( calib_down , axis = 0 )
|
|
||||||
|
|
||||||
peaks_calib_top = np.array(
|
|
||||||
utils.retrieve_peaks( mean_calib_top )
|
|
||||||
).astype( int )
|
|
||||||
peaks_calib_down = np.array(
|
|
||||||
utils.retrieve_peaks( mean_calib_down )
|
|
||||||
).astype( int )
|
|
||||||
|
|
||||||
"""
|
|
||||||
Signal without peaks
|
|
||||||
"""
|
|
||||||
|
|
||||||
def remove_peaks( signal , peaks ):
|
def remove_peaks( signal , peaks ):
|
||||||
|
"""
|
||||||
|
remove peaks from a signal
|
||||||
|
"""
|
||||||
peakless_signal = signal.copy()
|
peakless_signal = signal.copy()
|
||||||
for peak in peaks:
|
for peak in peaks:
|
||||||
first = peak
|
first = peak
|
||||||
|
|
@ -41,27 +21,30 @@ def remove_peaks( signal , peaks ):
|
||||||
peakless_signal[ first : peak ] = peakless_signal[ peak ] * np.ones( peak - first )
|
peakless_signal[ first : peak ] = peakless_signal[ peak ] * np.ones( peak - first )
|
||||||
|
|
||||||
return sig.medfilt( peakless_signal , 111 )
|
return sig.medfilt( peakless_signal , 111 )
|
||||||
peakless_calib_top = remove_peaks( mean_calib_top , peaks_calib_top )
|
|
||||||
|
|
||||||
def get_extremities( signal ):
|
def get_extremities( signal , peaks , temp ):
|
||||||
"""
|
"""
|
||||||
It's possible to have an idea of the part the spectrum begins with the
|
It's possible to have an idea of the part the spectrum begins with the
|
||||||
small large peak at the start of it, the peak at 3810 A should be inside
|
small large peak at the start of it, the peak at 3810 A should be inside
|
||||||
it. The same goes for the area at the end.
|
it. The same goes for the area at the end.
|
||||||
ONLY TRUE FOR THE CURRENT CALIBRATION (Hg)
|
ONLY TRUE FOR THE CURRENT CALIBRATION (Hg)
|
||||||
"""
|
"""
|
||||||
argmin = np.argmin(
|
peaks_inside , i = [] , 0
|
||||||
signal[ : len( signal ) // 2 ]
|
argmin = len( signal ) // 2
|
||||||
)
|
while len( peaks_inside ) != 1:
|
||||||
argmax = np.argmax( peakless_calib_top[ : argmin ] )
|
if i > 50 or len( signal[ : argmin ] ) == 0:
|
||||||
peaks_inside = np.where(
|
raise Exception( 'unknown plage, cannot autocalibrate' )
|
||||||
np.logical_and(
|
argmin = np.argmin(
|
||||||
argmax < peaks_calib_top,
|
signal[ : argmin ]
|
||||||
argmin > peaks_calib_top,
|
|
||||||
)
|
)
|
||||||
)[0]
|
argmax = np.argmax( signal[ : argmin ] )
|
||||||
if len( peaks_inside ) < 1:
|
peaks_inside = np.where(
|
||||||
raise Exception( 'unknown plage, cannot autocalibrate' )
|
np.logical_and(
|
||||||
|
argmax < peaks,
|
||||||
|
argmin > peaks,
|
||||||
|
)
|
||||||
|
)[0]
|
||||||
|
i += 1
|
||||||
first_peak = peaks_inside[0]
|
first_peak = peaks_inside[0]
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
|
@ -70,85 +53,48 @@ def get_extremities( signal ):
|
||||||
ONLY TRUE FOR THE CURRENT CALIBRATION (Hg)
|
ONLY TRUE FOR THE CURRENT CALIBRATION (Hg)
|
||||||
"""
|
"""
|
||||||
argmin_1 = np.argmin(
|
argmin_1 = np.argmin(
|
||||||
peakless_calib_top[ len( peakless_calib_top ) // 2 : ]
|
signal[
|
||||||
) + len( peakless_calib_top ) // 2
|
len( signal ) // 2 :
|
||||||
|
- int(
|
||||||
|
0.1 * len( signal ) // 2
|
||||||
|
)
|
||||||
|
] # not at the end
|
||||||
|
) + len( signal ) // 2
|
||||||
peaks_inside = np.where(
|
peaks_inside = np.where(
|
||||||
argmin_1 < peaks_calib_top
|
argmin_1 < peaks,
|
||||||
)[0]
|
)[0]
|
||||||
if len( peaks_inside ) < 1:
|
if len( peaks_inside ) < 1:
|
||||||
raise Exception( 'unknown plage, cannot autocalibrate' )
|
raise Exception( 'unknown plage, cannot autocalibrate' )
|
||||||
return ( first_peak , peaks_inside[0] )
|
return ( first_peak , peaks_inside[0] )
|
||||||
|
|
||||||
first_peak_ref = 3 # hard-coded
|
def only_keep_calib( peaks_data , peaks_calib ):
|
||||||
last_peak_ref = 20 # hard-coded
|
"""
|
||||||
first_peak_cal , last_peak_cal = get_extremities( peakless_calib_top )
|
only keep data peaks corresponding to calibration
|
||||||
|
"""
|
||||||
|
diff_calib = ( peaks_calib[ 1 : ] - peaks_calib[ : -1 ] ).astype( float )
|
||||||
|
diff_calib -= np.min( diff_calib )
|
||||||
|
diff_calib /= np.max( diff_calib )
|
||||||
|
|
||||||
"""
|
diff_data = ( peaks_data[ 1 : ] - peaks_data[ : -1 ] ).astype( float )
|
||||||
delete if too much peak in calib
|
diff_data -= np.min( diff_data )
|
||||||
"""
|
diff_data /= np.max( diff_data )
|
||||||
|
|
||||||
peaks_calib_top = peaks_calib_top[ first_peak_cal : last_peak_cal ]
|
peaks = [ -1 ]
|
||||||
peaks_reference = peaks_reference[ first_peak_ref : last_peak_ref ]
|
for i in range( len( diff_calib ) ):
|
||||||
|
good , sum_ = -1 , 0
|
||||||
|
for j in range( peaks[ - 1 ] + 1 , len( diff_data ) ):
|
||||||
|
sum_ += diff_data[ j ]
|
||||||
|
if sum_ - diff_calib[ i ] > 0.002:
|
||||||
|
print( sum_ - diff_calib[ i ] )
|
||||||
|
raise Exception( 'reference peak not found' )
|
||||||
|
if sum_ - diff_calib[ i ] > - 0.002:
|
||||||
|
good = j
|
||||||
|
break
|
||||||
|
if good == -1:
|
||||||
|
raise Exception( 'reference peak not found and not exceeded' )
|
||||||
|
peaks.append( good )
|
||||||
|
peaks.append( peaks[-1] + 1 ) # append the last peak
|
||||||
|
|
||||||
diff_ref = ( peaks_reference[ 1 : ] - peaks_reference[ : -1 ] ).astype( float )
|
return np.array(
|
||||||
diff_ref -= np.min( diff_ref )
|
[ peaks_data[ i ] for i in peaks[ 1 : ] ] # remove the first -1 value
|
||||||
diff_ref /= np.max( diff_ref )
|
).astype( int )
|
||||||
diff_cal = ( peaks_calib_top[ 1 : ] - peaks_calib_top[ : -1 ] ).astype( float )
|
|
||||||
diff_cal -= np.min( diff_cal )
|
|
||||||
diff_cal /= np.max( diff_cal )
|
|
||||||
|
|
||||||
peaks = [ -1 ]
|
|
||||||
for i in range( len( diff_ref ) ):
|
|
||||||
good , sum_ = -1 , 0
|
|
||||||
for j in range( peaks[ - 1 ] + 1 , len( diff_cal ) ):
|
|
||||||
sum_ += diff_cal[ j ]
|
|
||||||
if sum_ - diff_ref[ i ] > 0.002:
|
|
||||||
print( sum_ - diff_ref[ i ] )
|
|
||||||
raise Exception( 'reference peak not found' )
|
|
||||||
if sum_ - diff_ref[ i ] > - 0.002:
|
|
||||||
good = j
|
|
||||||
break
|
|
||||||
if good == -1:
|
|
||||||
raise Exception( 'reference peak not found and not exceeded' )
|
|
||||||
peaks.append( good )
|
|
||||||
peaks.append( peaks[-1] + 1 )
|
|
||||||
|
|
||||||
peaks_calib_top = np.array(
|
|
||||||
[ peaks_calib_top[ i ] for i in peaks[ 1 : ] ]
|
|
||||||
).astype( int )
|
|
||||||
|
|
||||||
"""
|
|
||||||
Calibration
|
|
||||||
"""
|
|
||||||
|
|
||||||
polyval_wavelength = np.polyfit(
|
|
||||||
peaks_calib_top,
|
|
||||||
peaks_reference,
|
|
||||||
1 ,
|
|
||||||
)
|
|
||||||
|
|
||||||
wavelength = np.polyval(
|
|
||||||
polyval_wavelength,
|
|
||||||
np.arange(
|
|
||||||
0 ,
|
|
||||||
len( mean_calib_top ),
|
|
||||||
1 ,
|
|
||||||
)
|
|
||||||
)
|
|
||||||
|
|
||||||
manager = qt.FigureManager(
|
|
||||||
qt.FigureCanvas(
|
|
||||||
fig.Figure(
|
|
||||||
figsize = ( 10 , 5 )
|
|
||||||
),
|
|
||||||
),
|
|
||||||
0,
|
|
||||||
)
|
|
||||||
manager.canvas.figure.gca().plot(
|
|
||||||
peaks_reference,
|
|
||||||
peaks_reference - np.polyval( polyval_wavelength , peaks_calib_top ),
|
|
||||||
'o-' ,
|
|
||||||
)
|
|
||||||
manager.show()
|
|
||||||
manager.start_main_loop()
|
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue