from typing import List,Dict,Union,Tuple
import math
import numpy as np
import scipy.optimize
from show_config.show_configuration import Parameters
from .level0imagearray import Level0_ImageCollection, Level0_Image, ImageStats
#------------------------------------------------------------------------------
# Level1_Spectra
#------------------------------------------------------------------------------
[docs]class Level1_Spectra:
"""
Instance Attributes:
.. py:attribute:: mjd
The UTC time of the image represented as a modified julian date. Float.
.. py:attribute:: exposure_time
The exposure time in micro-seconds seconds. Float.
.. py:attribute:: sensor_temp
The detector temperature in Celsius. Float.
.. py:attribute:: sensor_pcb_temp
The detector printed circuit board temperature in Celsius. Float.
.. py:attribute:: top_intferom_temp
The temperature in Celsius of the top of the interferometer. Float.
.. py:attribute:: bottom_intferom_temp
The temperature in Celsius of the bottom of the interferometer. Float
.. py:attribute:: optobox_temp
The temperature in Celsius of the optics box. Float
.. py:attribute:: Q7_temp
The temperature in Celsius of Q7. Float
.. py:attribute:: high_gain
The detector gain setting. The value is `True` for high gain and `False` for low gain.
.. py:attribute:: comment
An optional comment string for each image. The comment is typically supplied by the operator during data acquisition.
.. py:attribute:: spectrum
The Level 1 spectral image expressed as a numpy 2-D float array of dimension (H,M) where H is the number of height rows and M is the number
of transform frequencies, typically half the number of interfrogram pixels.
.. py:attribute:: error
The error on the Level 1 spectrum . The error field may be None meaning no error value is available. If it is not None then it will be
a numpy 2-D float array of dimension (H,M) where H is the number of height rows and M is the number of transform fequencies.
It will be the the same size as the spectrum attribute.
.. py:attribute:: ifgram_bounds
The sub-window bounds used to select the useful detector area from the original interferogram. This is a four elelemnt tuple (x0,x1,y0,y1)
"""
def __init__(self):
self.mjd = 0.0
self.exposure_time = 0.0
self.sensor_temp = 0.0
self.setpoint_temp = 0.0
self.sensor_pcb_temp = 0.0
self.top_intferom_temp = 0.0
self.bottom_intferom_temp = 0.0
self.optobox_temp = 0.0
self.Q7_temp = 0.0
# self.tec_on = 0.0
self.high_gain = 0.0
self.comment = None
self.spectrum = None
self.error = None
self.ifgram_bounds = []
#------------------------------------------------------------------------------
# class L0Algorithms:
#------------------------------------------------------------------------------
[docs]class L0Algorithms:
def __init__(self, parameters: Parameters ) -> None:
self.params = parameters # type: Parameters
#------------------------------------------------------------------------------
# L0Algorithms::rms_signal_error_from_detector_specs
#------------------------------------------------------------------------------
[docs] def rms_signal_error_from_detector_specs(self, dnsignal: np.ndarray ) ->np.ndarray:
"""
Fetch the theoretical error on the signal read from the detector in DN. Considers Poisson counting statistics
on total number of electrons an dthe detector readout noise.
:param dnsignal: The signal readout from the detector. Note that you may have to subtract any DC bias
ont the detector before calling this routine
:return: The calculated error in DN.
"""
Edn = self.params.config['level0']['electrons_per_DN'] # Get electrons per DN
Rde = self.params.config['level0']['detector_readout_noise'] # get readout noise in electrons
Ne = dnsignal*Edn # get the number of electrons in signal, (this is the square of the error!)
Netot = np.sqrt( Ne + Rde*Rde ) # get the total electron error, Ne poisson error^2 + readout error^2
dnerror = Netot/Edn # convert total electron error to DN
return dnerror # return the dnerror
#------------------------------------------------------------------------------
# L0Algorithms::removedc_bias_from_apodized_interferogram
#------------------------------------------------------------------------------
[docs] def apodized_interferogram_with_zero_bias(self, x: np.ndarray, y : np.ndarray ) -> np.ndarray:
"""
Apply an apodization function, *y*, to the each height row of the interferogram, *x*. Apply
a dc offset to the interferogram, *x*, such that the average of
the product of the interferogram, *x*, and the apodization function, *y*, is zero. This eliminates the large
zero order component bleeding into nearby frequencies after we perform the FFT.
Note that removing DC bias from the interferogram before apodization is applied results in a slight non-zero
dc component that makes a quite large spike in the zero harmonic which bleeds into neighbouring frequencies. This technique
removes that bias and helps avoid bleed through of the zero harmonic in the fft spectrum.
:param x: Original interferogram, a 2d-array of size (H,M). We need to average over M
:param y: a 2d array of size (H,M). We need to average over M
:return: the apodized integrforogram with zero bias for each height row, 2-d array (H,M)
"""
H,M = x.shape
xy = x*y
xyav = np.average( xy, 1) # Get the average signal of the apodized signal at each height, array (H)
yav = np.average( y, 1) # Get the average value of the apodization function at each height, array (H)
b = xyav/yav # Get the DC correction to subtrcat from original signal to zero bias apodized signal, array (H)
bav = np.tile(b, (M, 1)).transpose() # make the DC correction into 2-D array
xy = (x-bav)*y
return xy
#------------------------------------------------------------------------------
# L0Algorithms::rms_frequency_error_from_rms_spatial_error
#------------------------------------------------------------------------------
[docs] def rms_frequency_error_from_rms_spatial_error(self, spatial_error_array: np.ndarray, include_apodization=True, real_component_only= False ) -> Tuple[np.ndarray, np.ndarray]:
"""
Calculates the theoretical root mean square error in Fourier transform space given the root mean square error in interferogram space. The algorithm is applied to
each height row independently. For each height row only one value, the RMS error is returned.
:param spatialerrorarray: 2D array (H,M) of errors for M points at each height (H) row.
:param include_apodization: Default True. If True then divide theoretrical RMS error by a factor of 2
:param real_component_only: Default False. If True then divide theoretical RMS errror by sqrt(2)
:return: the theoretical root mean square error at each height row. An array (H)
"""
H,N = spatial_error_array.shape
Ex = np.sqrt( np.sum( spatial_error_array*spatial_error_array , 1)/N ) # Get the rms error in spatial space, Ex(H)
factor = 0.5*math.sqrt(N) if include_apodization else math.sqrt(N) # get the factor for the fourier transform RMS error
if real_component_only: factor = factor/math.sqrt(2) # See if we are getting real+imag or only real component
Ef = factor*Ex # and get the RMS error in transform space.
return Ef, Ex
#------------------------------------------------------------------------------
# L0Algorithms::interferogram_subwindow
#------------------------------------------------------------------------------
[docs] def interferogram_subwindow(self, level0: np.ndarray, userdefinederror : np.ndarray=None) -> Tuple[np.ndarray, np.ndarray, Tuple[int, int, int, int]]:
"""
Fetch the interferogram image and error arrays sub-windowed into the useful detector area.
:param level0: The incoming :class:`~showapi.level0.level0imagearray.Level0_Image` record.
:param userdefinederror: Default None. If not None then derive the error on the record's image from this array rathe rthan the *error* field stored in the record
:return:
A three element tuple containing, (i) the sub-windows image, (ii) the sub-windowe error and (iii) the detector bounding region as (x0,x1,y0,y1).
"""
x0 = self.params.fft_window_start
x1 = self.params.fft_window_end
y0 = self.params.height_window_start
y1 = self.params.height_window_end
ifgram = np.copy( level0[ y0:(y1+1), x0:(x1+1)]) # Select the subset of pixels used for the interferogram
e = userdefinederror # select the inteferogram error source
error = np.copy( e[y0:(y1+1), x0:(x1+1)]) if e is not None else None # sub-window the inteferogram error if its not none
return ifgram, error, (x0,x1,y0,y1)
#------------------------------------------------------------------------------
# L0Algorithms::subwindow_array
#------------------------------------------------------------------------------
[docs] def subwindow_array(self, image: np.ndarray) -> np.ndarray:
"""
Fetch this instruments sub-window from a given image.
:param image: An incoming 2-D array which will be sub-windowed. It is assumed it is in the same shape as the detector
:return: The sub-windowed image as a 2-D array
"""
x0 = self.params.fft_window_start
x1 = self.params.fft_window_end
y0 = self.params.height_window_start
y1 = self.params.height_window_end
ifgram = np.copy( image[ y0:(y1+1), x0:(x1+1)]) # Select the subset of pixels used for the interferogram
return ifgram
#------------------------------------------------------------------------------
# L0Algorithms::interferogram_to_spectrum
#------------------------------------------------------------------------------
[docs] def interferogram_to_spectrum( self, level0 : np.ndarray, userdefinederror : np.ndarray = None )-> Tuple[ np.ndarray, np.ndarray, np.ndarray, np.ndarray] :
"""
Converts the interferogram to spectra by applying an apodization (Hanning) window and taking the FFT at each
selected height level. A sub-window (see below) is selected to create clean interferograms of dimension (H,M).
The clean interferogram is apodized and fourier transformed to make a spectral 2-D image of size (H, M/2+1), ie we remove
frequencies above the nyquist limit. The new spectral image is returned as a :class:`~showapi.level0.l0algorithms.Level1_Spectra` object.
The code selects a sub-window as defined by the instrument parameters fft_window_start,
fft_window_end, height_window_start and height_window_end. This sub-window is meant to eliminate all the useless
edge areas of the detector. There are no requirements on the size of the sub-window, for example it does not have to
be a power of 2 etc. Users should eliminate bad and marginal regions near the edges of the detector with the sub-window
as this is much better leaving them in where they ultimately corrupt the entire signal.
The code removes a DC bias from the original intefreogram such that the zero order harmonic of the FFT is zero. There
is a subtlety in this process as we remove a DCbias from the originl signal so the average of the product of the
orginal signal times the apodization function is zero. This is not quite the same as ensuring the average of the
original signal is zero. If this correction is not made then we see the zero harmonic get quite large and it bleeds out
into neighbouring spectral pixels during the FFT.
:param level0:
A level 0 record, :class:`~showapi.level0.level0imagearray.Level0_Image`. This contains both the image which is transformed and a header. Both header and transormed image are copied to the
level 1 spectral object, :class:`~showapi.level0.l0algorithms.Level1_Spectra`. By default an error estimate of the transform is made from the error estimate in the Level 0 record. The error estimate can
be overridden using a user defined error estimate
:param userdefinederror:
Default = None. A user defined 2-d array specifying the error on the Level 0 interferogram image. If provided this array will be used in the error analysis and propagation instead of
the Level 0 records error field. 2-d image.
:return:
A three element tuple containing (i) the desired :class:`~showapi.level0.l0algorithms.Level1_Spectra` object, (ii) the windowed and apodized interferogram, array(H,M) and (iii) the windowed but not apodized interferogram, array (H,M)
"""
ifgram, iferror, bounds = self.interferogram_subwindow(level0, userdefinederror=userdefinederror) # Select the subset of pixels used for the interferogram, and get the bounding region indices
H,M = ifgram.shape # get the shape of the sub-window
han = np.tile(np.hanning(M), (H, 1)) # Select the apodization function for the FFT
ifgramwh = self.apodized_interferogram_with_zero_bias( ifgram, han ) # apply apodization and zero bias the average of the product
spectra = np.fft.rfft( ifgramwh ) # Compute the real FFT. The real fft does not compute the negative frequencies
if iferror is not None:
EFrms = np.sqrt( np.sum( iferror*iferror, 1 )/2) # Get the root mean square error in fourier transform space for each height, array (H)
HH,FF = spectra.shape # get the shape of the spectra, we want F as its about half of M
error = np.tile( EFrms, (FF,1)).transpose() # Now tile the root mean square error across the F spectral elements at each height
else:
error = None
return spectra, error, ifgramwh, ifgram
# ------------------------------------------------------------------------------
# _cosinefunc( param, M, s):
# ------------------------------------------------------------------------------
@staticmethod
def _cosinefunc(param, M, s):
"""
Evaluates the difference between a cosine function modulated by gaussian envelope and measured spectrum ``s``.
:param M:
:param s:
:return:
"""
A = param[0]
w = param[1]
phi = param[2]
K = param[3]
x0 = param[4]
sigma = param[5]
x = np.arange(0, M)
sigma2 = sigma * sigma
x2 = np.square( x-x0)
g = np.exp(x2/(-2 * sigma2))
y = (A * np.cos(w * x + phi)*g + K) - s
return y
# ------------------------------------------------------------------------------
# _cosinefunc_jacobian(param, M, s):
# ------------------------------------------------------------------------------
@staticmethod
def _cosinefunc_jacobian(param, M, s):
"""
Generates a cosine function with Gaussian envelope used to fit to the Krypton
fringe data
:param M: A tuple of parameters
:param s: The measured spectrum
:return: The diffreence between modelled function and measured spectrum
"""
A = param[0]
w = param[1]
phi = param[2]
K = param[3]
x0 = param[4]
sigma = param[5]
x = np.arange(0, M)
J = np.zeros([M, 6])
sigma2 = sigma*sigma
x1 = (x-x0)
x2 = np.square(x1)
g = A*np.exp( x2/(-2*sigma2) )
y = np.cos(w * x + phi)*g
z = -np.sin(w * x + phi)*g
J[:, 0] = y/A
J[:, 1] = z*x
J[:, 2] = z
J[:, 3] = np.zeros([M]) + 1
J[:, 4] = y*x1/(sigma2)
J[:, 5] = y*x2/(sigma2*sigma)
return J
# ------------------------------------------------------------------------------
# def fit_cosine_with_exponential(self, algo: L0Algorithms, krdata: Level0_Image):
# ------------------------------------------------------------------------------
[docs] def fit_cosine_with_gaussian_envelope(self, krdata_subwindow: np.ndarray ):
"""
Fits a cosine with a gaussian enevelope to each row of a sub-windowed interferogram using a least squares approach
:param krdata_subwindow:
A 2-D interferogram sub-window of size (H,M) where H is the number of height rows and M is the number of interferogram bins.
The cosine and gaussian enevelope is fitted to each row of the sub-window.
:return:
a 2-D array of size (6,H) where H is the number of heightrows and 6 is the number of fitted parameters.
Theory: The code fits the following function to the measured interferogram,
.. math::
y = A\\cos (\\omega x + \\phi)e^{-\\frac{(x-x_0)^2}{2\\sigma^2} } + K
where :math:`x` is the inteferogram bin number and we fit for the following 6 variables,
1. :math:`A`
2. :math:`\\omega`
3. :math:`\\phi`
4. :math:`K`
5. :math:`x_0`
6. :math:`\\sigma`
The least squares algorithm uses the analytic differentials to find the best fit parameters. Let
.. math::
z = A\\sin (\\omega x + \\phi)e^{-\\frac{(x-x_0)^2}{2\\sigma^2}}
then the differentials are given by
.. math::
\\frac{\\partial y}{\\partial A} &= \\frac{y-K}{A} \\\\
\\frac{\\partial y}{\\partial \\omega} &= -zx \\\\
\\frac{\\partial y}{\\partial \\phi} &= -z \\\\
\\frac{\\partial y}{\\partial K} &= 1 \\\\
\\frac{\\partial y}{\\partial x_0} &= (y-K)\\frac{(x-x_0)}{\\sigma^2} \\\\
\\frac{\\partial y}{\\partial \\sigma} &= (y-K)\\frac{(x-x_0)^2}{\\sigma^3} \\\\
"""
image = krdata_subwindow
H, M = image.shape # Get the Number of heights and number of interferogram pixel
han = np.tile(np.hanning(M), (H, 1)) # Select the apodization function for the FFT
ifgramwh = self.apodized_interferogram_with_zero_bias(image, han) # apply apodization and zero bias the average of the product
spectra = np.fft.rfft(ifgramwh) # and get the real Fourier transform
x = np.arange(0, M, 1) # get the values of x for each interferogram bin
fitimage = np.zeros( image.shape )
fitresults = np.zeros( [6,H] )
for ih in range(H): # For each valid height
s = image[ih, :] # get the Krypton data for this height
ffts = spectra[ih, :] # Get the spectrum for this height
imax = np.argmax(np.absolute(ffts)) # Find the maximum frequency in the spectrum. This is probably Krypton
w = 2 * math.pi * imax / M # Find the nominal frequency of this component.
phi = np.angle(ffts[imax]) # Get the phase of this compoennts
A = 4 * np.absolute(ffts[imax]) / M # and the magnirude of the cosine
K = np.average(s) # Guess the nominal background signal
xx0 = np.argmax(s) # Centre of gaussian envelope is on the maxim interferogram signal
sigma = M/2 # Guess the standard deviation is half the width
x0 = [A, w, phi, K, xx0, sigma ] # Set up the least squares vector
results = scipy.optimize.least_squares(L0Algorithms._cosinefunc, x0, jac=L0Algorithms._cosinefunc_jacobian, args=[M, s]) # and so the least squares fit.
x = results.x # get the results from the fit.
fitimage[ih,:] = results.fun
#print(results.message) # NEED TO CHECK THAT TEH FIT WORKED.
#print('A initial %15.7f ---> %15.7f' % (x0[0], x[0]))
#print('w initial %15.7f ---> %15.7f' % (x0[1], x[1]))
#print('Phi initial %15.7f ---> %15.7f' % (x0[2], x[2]))
#print('K initial %15.7f ---> %15.7f' % (x0[3], x[3]))
#print('x0 initial %15.7f ---> %15.7f' % (x0[4], x[4]))
#print('sigma initial %15.7f ---> %15.7f' % (x0[5], x[5]))
#res = { 'A': x[0], 'w' : x[1], 'phi':x[2], 'K' : x[3], 'x0':x[4], 'sigma':x[5] }
fitresults[:,ih] = x
fitimage += image
return fitresults, fitimage