From bf40d1db67d853fb54a107c27717060c042134cd Mon Sep 17 00:00:00 2001 From: linarphy Date: Tue, 23 May 2023 14:35:21 +0200 Subject: [PATCH] Add auto calibration based on peakless signal --- wavelength_calibration.py | 180 +++++++++++++++++++++++--------------- 1 file changed, 109 insertions(+), 71 deletions(-) diff --git a/wavelength_calibration.py b/wavelength_calibration.py index c526d46..3e954c6 100644 --- a/wavelength_calibration.py +++ b/wavelength_calibration.py @@ -13,20 +13,109 @@ reference = np.load( 'asset/reference.npy' ) mean_calib_top = np.mean( calib_top , axis = 0 ) mean_calib_down = np.mean( calib_down , axis = 0 ) -peaks_calib_top = utils.retrieve_peaks( mean_calib_top ) -peaks_reference = utils.retrieve_peaks( reference[1] , window_size = 1 , max_window_size = 1 ) - -polyval = np.polyfit( peaks_calib_top , peaks_reference[ 2 : ] , 1 ) - -peaks_values = [ reference[ 1 , i ] for i in np.array( peaks_reference ).astype( int ) ] -sorting = np.argsort( peaks_values )[ :: -1 ] -wavelength = [ reference[ 0 , i ] for i in np.array( [ peaks_reference[ j ] for j in sorting ] ).astype( int ) ] - -polyval_wavelength = np.polyfit( peaks_reference , np.sort( wavelength ) , 1 ) - -wavelength = np.polyval( polyval_wavelength , np.polyval( polyval , np.arange( 0 , len( mean_calib_top ) , 1 ) ) ) +peaks_calib_top = np.array( + utils.retrieve_peaks( mean_calib_top ) +).astype( int ) +peaks_reference = np.array( + utils.retrieve_peaks( + reference[1] , + window_size = 1 , + max_window_size = 1, + ) +).astype( int ) """ +Signal without peaks +""" + +peakless_calib_top = mean_calib_top.copy() +for peak in peaks_calib_top: + first = peak + peak , old = first - 2 , first - 1 + while peakless_calib_top[ peak ] <= peakless_calib_top[ old ]: + old = peak + peak -= 1 + peakless_calib_top[ peak : first ] = peakless_calib_top[ peak ] * np.ones( first - peak ) + + peak , old = first + 2 , first + 1 + while peakless_calib_top[ peak ] <= peakless_calib_top[ old ]: + old = peak + peak += 1 + peakless_calib_top[ first : peak ] = peakless_calib_top[ peak ] * np.ones( peak - first ) + +peakless_calib_top = sig.medfilt( peakless_calib_top , 111 ) + +""" +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 +it. The same goes for the area at the end. +ONLY TRUE FOR THE CURRENT CALIBRATION (Hg) +""" +argmin = np.argmin( + peakless_calib_top[ : len( peakless_calib_top ) // 2 ] +) +argmax = np.argmax( peakless_calib_top[ : argmin ] ) +peaks_inside = np.where( + np.logical_and( + argmax < peaks_calib_top, + argmin > peaks_calib_top, + ) +)[0] +if len( peaks_inside ) < 1: + raise Error( 'unknown plage, cannot autocalibrate' ) +first_peak_cal = peaks_inside[0] +first_peak_ref = 3 # hard-coded + +""" +The next peak after the minimum at the end of the spectrum should be +5079 A. +ONLY TRUE FOR THE CURRENT CALIBRATION (Hg) +""" +argmin_1 = np.argmin( + peakless_calib_top[ len( peakless_calib_top ) // 2 : ] +) + len( peakless_calib_top ) // 2 + +peaks_inside = np.where( + argmin_1 < peaks_calib_top +)[0] +if len( peaks_inside ) < 1: + raise Error( 'unknown plage, cannot autocalibrate' ) +last_peak_cal = peaks_inside[0] +last_peak_ref = 20 # hard-coded + +print( first_peak_cal , last_peak_cal ) +polyval = np.polyfit( + peaks_calib_top[ first_peak_cal : last_peak_cal ], + peaks_reference[ first_peak_ref : last_peak_ref ], + 1 , +) # We suppose there is as much calib peak than ref peak for now + +peaks_values = [ reference[ 1 , i ] for i in peaks_reference ] +sorting = np.argsort( peaks_values )[ :: -1 ] +wavelength = np.array( [ + reference[ 0 , i ] for i in np.array( [ + peaks_reference[ j ] for j in sorting + ] ).astype( int ) +] ) + +polyval_wavelength = np.polyfit( + peaks_reference, + np.sort( wavelength ), + 1, +) + +wavelength = np.polyval( + polyval_wavelength, + np.polyval( + polyval, + np.arange( + 0 , + len( mean_calib_top ), + 1 , + ) + ) +) + manager = qt.FigureManager( qt.FigureCanvas( fig.Figure( @@ -35,73 +124,22 @@ manager = qt.FigureManager( ), 0, ) - -manager.canvas.figure.add_subplot( 2 , 1 , 1 , xlim = ( 3800 , 4000 ) , xmargin = 0 ).plot( +manager.canvas.figure.gca().plot( wavelength , mean_calib_top, ) -manager.canvas.figure.axes[0].set_title( 'raw data' ) - -manager.canvas.figure.add_subplot( 2 , 1 , 2 , xlim = ( 3800 , 4000 ) , xmargin = 0 ).plot( - reference[0], - reference[1], -) -manager.canvas.figure.axes[1].set_title( 'reference' ) - -manager.show() -manager.start_main_loop() -""" -for peak in np.array( peaks_calib_top ).astype( int ): - first = peak - peak , old = first - 2 , first - 1 - while mean_calib_top[ peak ] <= mean_calib_top[ old ]: - old = peak - peak -= 1 - mean_calib_top[ peak : first ] = mean_calib_top[ peak ] * np.ones( first - peak ) - - peak , old = first + 2 , first + 1 - while mean_calib_top[ peak ] <= mean_calib_top[ old ]: - old = peak - peak += 1 - mean_calib_top[ first : peak ] = mean_calib_top[ peak ] * np.ones( peak - first ) - -manager = qt.FigureManager( - qt.FigureCanvas( - fig.Figure( - figsize = ( 10 , 5 ) - ), - ), - 0, -) manager.canvas.figure.gca().plot( - mean_calib_top -) -manager.canvas.figure.gca().plot( - sig.medfilt( mean_calib_top , 71 ) -) -manager.show() -manager.start_main_loop() - -manager = qt.FigureManager( - qt.FigureCanvas( - fig.Figure(), - ), - 1, -) -signal_calib = sig.convolve( - sig.medfilt( mean_calib_top , 71 ), - np.ones( 50 ), - 'same' -) -manager.canvas.figure.gca().plot( - signal + wavelength , + peakless_calib_top, ) manager.canvas.figure.gca().vlines( [ - np.argmin( signal[ 25 : -25 ] ) + wavelength[ argmin ] , + wavelength[ argmax ] , + wavelength[ argmin_1 ], ] , - np.min( signal[ 25 : -25 ] ), - np.max( signal[ 25 : -25 ] ), + np.min( peakless_calib_top ), + np.max( peakless_calib_top ), color = 'red' , ) manager.show()