The Working Horse

The Schottky signal is of noisy nature owing to the indeterministic arrival times of the charged particles at the detector. Its randomness always makes human feel difficult to pictorially perceive and quantitatively describe. Fortunately with the Wiener-Khinchin theorem [1,2], one is able to attain the Power Spectral Density (PSD) of the Schottky signal in the frequency domain, which is deterministic and thus easy to interpret. The "internal structure" of the signal has never been revealed in a more straightforward way.

In reality, the computation of the PSD may not be viable, since the ensemble mean is needed, in contrast to only one instance of the ensemble is available in most cases. Hence the Spectral Density Esimation (SDE) becomes the central topic for the practical applications of the Wiener-Khinchin theorem. Beginning with the concept of the Periodogram proposed by A. Schuster more than a century ago [3], a number of different methods has emerged ever since. The representative ones include replacing the ensemble mean with the time average under the assumption of ergodicity, or applying a set of orthonormal windows to the same sequence to improve the time resolution of the estimation.

The former evolves to two similar average schemes: Bartlett's [4] and Welch's [5], while the latter mainly leads to the Thomson's Multitaper methods [6]. One of the superiorities of the multitaper method over the time-average method, among many others such as the non-subjective way to define the frequency resolution and the possibility to trade the estimation bias for the variance, is manifested in that the multitaper can inherently embrace the statistical F-test for the line peak detection [7]. Besides those, all methods possess a well-defined parameter called Degrees of Freedom, which is crucial to draw the Confidence Band of the SDE at a given Confidence Level [7].

Table of Contents

• Class Structure
  • fft_1d
  • spectrum
  • periodogram_1d
  • multitaper_1d
  • adaptive_multitaper_1d
  • time_average_1d
  • fft_2d
  • spectrogram
  • periodogram_2d
  • multitaper_2d
  • adaptive_multitaper_2d
  • time_average_2d
  • confidence_band
  • plot_1d
  • plot_2d
• References

Class Structure

The class Processing developed in the script processing.py tries to address the aforementioned SDE computations with the primary aim at an actual data file as the input. It is an inheritance of the class Preprocessing, which is detailed on a separate page: Loading Files. The instantiation of the class DPSS, which is described on another page: DPSS, is also indispensable for the multitaper estimations.

Inside its core, the Fast Fourier Transform (FFT) is heavily used for the efficient execution of the code. The implementation of FFT is CPU based, by the "fastest" library FFTW. By default, all the CPU threads are utilized to distribute the computational burden. This can be changed by tweaking the class attribute n_thread.

Processing(file_str)

Instantiate the class.

• arguments:
  • file_str: str
    Name of the data file. Two file formats (wv & tiq) are supported, which are distinguished from the file extensions. *.wvh and *.wvd both indicate a wv file.

fft_1d

fft_1d(window_length=1000, n_offset=0, padding_ratio=0, window=None, beta=None)

Compute the Fourier transform of an, optionally, zero-padded and windowed time series.

• arguments:
  • window_length: int, default 1000
    Length of the time series for one FFT frame.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 0
    Ratio of the zero-padded frame length to the original series length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the series length.
  • window: str, default None
    Tapering window to be applied onto the time series. Valid options are ["bartlett", "blackman", "hamming", "hanning", "kaiser"]. If not specified, a boxcar window, or equivalently no window, is acknowledged.
  • beta: float, default None
    Only if the kaiser window is specified, a complementary argument beta is also expected to define the window shape.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. If the frame length is odd, the 0 frequency is in the middle of the array; if even, the negative Nyquist frequency is in the beginning.
  • spectrum: 1d numpy array
    A complex array of Fourier transformed values in units of volt. The array is re-ordered according to the correspondence with the frequencies.

spectrum

spectrum(window_length=1000, n_offset=0, padding_ratio=2, window=None, beta=None)

Compute the amplitude of the Fourier transform of an, optionally, zero-padded and windowed time series.

• arguments:
  • window_length: int, default 1000
    Length of the time series for one FFT frame.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original series length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the series length.
  • window: str, default None
    Tapering window to be applied onto the time series. Valid options are ["bartlett", "blackman", "hamming", "hanning", "kaiser"]. If not specified, a boxcar window, or equivalently no window, is acknowledged.
  • beta: float, default None
    Only if the kaiser window is specified, a complementary argument beta is also expected to define the window shape.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. If the frame length is odd, the 0 frequency is in the middle of the array; if even, the negative Nyquist frequency is in the beginning.
  • spectrum: 1d numpy array
    A real array of Fourier transform amplitudes in units of V/kHz. With such a voltage-density expression, the total power is the numerical integral, instead of summation, of the squared spectrum. The array is re-ordered according to the correspondence with the frequencies.

periodogram_1d

periodogram_1d(window_length=1000, n_offset=0, padding_ratio=2, window=None, beta=None)

Estimate the PSD of an, optionally, zero-padded time series by using the periodogram.

• arguments:
  • window_length: int, default 1000
    Length of the time series for one FFT frame.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original series length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the series length.
  • window: str, default None
    Tapering window to be applied onto the time series. Valid options are ["bartlett", "blackman", "hamming", "hanning", "kaiser"]. If not specified, a boxcar window, or equivalently no window, is acknowledged.
  • beta: float, default None
    Only if the kaiser window is specified, a complementary argument beta is also expected to define the window shape.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. If the frame length is odd, the 0 frequency is in the middle of the array; if even, the negative Nyquist frequency is in the beginning.
  • spectrum: 1d numpy array
    A real array of the periodogram estimate of the PSD in units of V²/kHz. The array is re-ordered according to the correspondence with the frequencies.
  • n_dof: int
    Degrees of freedom. By definition, always return 2.

multitaper_1d

multitaper_1d(window_length=1000, n_offset=0, padding_ratio=2, half_bandwidth=3, n_taper=4, f_test=False)

Estimate the PSD of an, optionally, zero-padded time series with the basic multitaper method.

• arguments:
  • window_length: int, default 1000
    Length of the time series for one FFT frame.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original series length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the series length.
  • half_bandwidth: float, default 3
    Half of the bandwidth of the DPSS tapers in units of frequency grid.
  • n_taper: int, default 4
    Number of DPSS tapers. To avoid spectral leakage, it is advisalbe to keep n_taper less than 2*half_bandwidth.
  • f_test: bool, default False
    If True, also compute F-test-relevant statistics.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. If the frame length is odd, the 0 frequency is in the middle of the array; if even, the negative Nyquist frequency is in the beginning.
  • spectrum: 1d numpy array
    A real array of the basic multitaper estimate of the PSD in units of V²/kHz. The array is re-ordered according to the correspondence with the frequencies.
  • n_dof: int
    Degrees of freedom. By definition, return 2*n_taper.
  • statistics: 1d numpy array, optional
    A real array of the F-test statistics whose elements correspond to the SDE in the spectrum. Only returned if f_test is True.
  • critical_value: float, optional
    A real critical value beyond which the statistics shall be considered to be significant. Only returned if f_test is True.

adaptive_multitaper_1d

adaptive_multitaper_1d(window_length=1000, n_offset=0, padding_ratio=2, half_bandwidth=3, n_taper=4, f_test=False)

Estimate the PSD of an, optionally, zero-padded time series with the adaptive multitaper method.

• arguments:
  • window_length: int, default 1000
    Length of the time series for one FFT frame.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original series length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the series length.
  • half_bandwidth: float, default 3
    Half of the bandwidth of the DPSS tapers in units of frequency grid.
  • n_taper: int, default 4
    Number of DPSS tapers. To avoid spectral leakage, it is advisalbe to keep n_taper less than 2*half_bandwidth.
  • f_test: bool, default False
    If True, also compute F-test-relevant statistics.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. If the frame length is odd, the 0 frequency is in the middle of the array; if even, the negative Nyquist frequency is in the beginning.
  • spectrum: 1d numpy array
    A real array of the adaptive multitaper estimate of the PSD in units of V²/kHz. The array is re-ordered according to the correspondence with the frequencies.
  • n_dof: 1d numpy array
    A real array of degrees of freedom whose elements correspond to the SDE in the spectrum.
  • statistics: 1d numpy array, optional
    A real array of the F-test statistics whose elements correspond to the SDE in the spectrum. Only returned if f_test is True.
  • critical_value: 1d numpy array, optional
    A real array of critical values beyond which the corresponding statistics shall be considered to be significant. Only returned if f_test is True.

time_average_1d

time_average_1d(window_length=1000, n_offset=0, padding_ratio=2, n_average=10, estimator='p', **kwargs)

Estimate the PSD of a time series by evenly sectioning it into several segments and averaging their SDEs with, optionally, zero-padding. The SDE can be computed with the periodogram, the basic multitaper, or the adaptive multitaper method. Usually the Bartlett's averaging scheme is adopted, except for the periodogram estimate with windows other than the boxcar. In such cases, the Welch's averaging scheme with half overlapping portions is used instead. For the adaptive multitaper method, the average is weighted by the degrees of freedom of each SDE.

• arguments:
  • window_length: int, default 1000
    Length of the series segment for one FFT frame.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original segment length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the segment length.
  • n_average: int, default 10
    Number of consecutive in time SDEs for average. For Welch's averaging scheme, the actual number is 2*n_average-1 by counting the overlapped segments.
  • estimator: str, default 'p'
    Estimating method to be adopted for the SDE. Valid options are ['p', 'm', 'a'], which stand for Periodogram, basic Multitaper, and Adaptive multitaper, respectively.
  • **kwargs:
    Some complementary arguments to be passed on to the backstage method depending on the estimator. If estimator='p', window, and possibly, beta should be specified. If estimator='m' | 'a', half_bandwidth and n_taper should be specified.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. If the frame length is odd, the 0 frequency is in the middle of the array; if even, the negative Nyquist frequency is in the beginning.
  • spectrum: 1d numpy array
    A real array of the averaged SDE in units of V²/kHz. The array is re-ordered according to the correspondence with the frequencies.
  • n_dof: int | 1d numpy array
    Degrees of freedom depending on the estimator. If estimator='p' | 'm', n_dof is a real number. If estimator='a', n_dof is a real array whose elements correspond to the SDE in the spectrum.

fft_2d

fft_2d(window_length=1000, n_frame=100, n_offset=0, padding_ratio=0, window=None, beta=None)

Compute the Fourier transforms of the segments of an evenly sectioned, and optionally, zero-padded and windowed time series.

• arguments:
  • window_length: int, default 1000
    Length of the series segment for one FFT frame.
  • n_frame: int, default 100
    Number of the FFT frames. It should, in principle, be positive. However, the actual n_frame may depend on the available IQ pairs left in the data block after deducting the n_offset. If n_frame is negative, all the available IQ pairs, after trimming off the incomplete segment, will be selected.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 0
    Ratio of the zero-padded frame length to the original segment length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the segment length.
  • window: str, default None
    Tapering window to be applied onto the segments. Valid options are ["bartlett", "blackman", "hamming", "hanning", "kaiser"]. If not specified, a boxcar window, or equivalently no window, is acknowledged.
  • beta: float, default None
    Only if the kaiser window is specified, a complementary argument beta is also expected to define the window shape.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. The length of the frequencies is one point longer than the frame length, such that the frequencies delimits the spectrogram in frequency. If the frame length is odd, the negative Nyquist frequency is in the beginning of the array; if even, the negative Nyquist frequency is the mean of the beginning two elements in the array.
  • times: 1d numpy array
    A real array of uniformaly distanced time grids in units of second. It is relative to the beginning of the data block. The length of the times is larger than the actual n_frame by one, such that the times delimits the spectrogram in time.
  • spectrogram: 2d numpy array
    A complex array of Fourier transformed values in units of volt. The array is re-ordered in frequency according to the correspondence with the frequencies, and stacked consecutively by frames.

spectrogram

spectrogram(window_length=1000, n_frame=100, n_offset=0, padding_ratio=2, window=None, beta=None)

Compute the amplitudes of the Fourier transforms of the segments of an evenly sectioned, and optionally, zero-padded and windowed time series.

• arguments:
  • window_length: int, default 1000
    Length of the series segment for one FFT frame.
  • n_frame: int, default 100
    Number of the FFT frames. It should, in principle, be positive. However, the actual n_frame may depend on the available IQ pairs left in the data block after deducting the n_offset. If n_frame is negative, all the available IQ pairs, after trimming off the incomplete segment, will be selected.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original segment length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the segment length.
  • window: str, default None
    Tapering window to be applied onto the segments. Valid options are ["bartlett", "blackman", "hamming", "hanning", "kaiser"]. If not specified, a boxcar window, or equivalently no window, is acknowledged.
  • beta: float, default None
    Only if the kaiser window is specified, a complementary argument beta is also expected to define the window shape.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. The length of the frequencies is one point longer than the frame length, such that the frequencies delimits the spectrogram in frequency. If the frame length is odd, the negative Nyquist frequency is in the beginning of the array; if even, the negative Nyquist frequency is the mean of the beginning two elements in the array.
  • times: 1d numpy array
    A real array of uniformaly distanced time grids in units of second. It is relative to the beginning of the data block. The length of the times is larger than the actual n_frame by one, such that the times delimits the spectrogram in time.
  • spectrogram: 2d numpy array
    A real array of Fourier transform amplitudes in units of V/kHz. With such a voltage-density expression, the total power in a frame is the numerical integral, instead of summation, of the squared spectrogram. The array is re-ordered in frequency according to the correspondence with the frequencies, and stacked consecutively by frames.

periodogram_2d

periodogram_2d(window_length=1000, n_frame=100, n_offset=0, padding_ratio=2, window=None, beta=None)

Estimate the PSDs of the segments of an evenly sectioned, and optionally, zero-padded time series by using the periodogram.

• arguments:
  • window_length: int, default 1000
    Length of the series segment for one FFT frame.
  • n_frame: int, default 100
    Number of the FFT frames. It should, in principle, be positive. However, the actual n_frame may depend on the available IQ pairs left in the data block after deducting the n_offset. If n_frame is negative, all the available IQ pairs, after trimming off the incomplete segment, will be selected.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original segment length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the segment length.
  • window: str, default None
    Tapering window to be applied onto the segments. Valid options are ["bartlett", "blackman", "hamming", "hanning", "kaiser"]. If not specified, a boxcar window, or equivalently no window, is acknowledged.
  • beta: float, default None
    Only if the kaiser window is specified, a complementary argument beta is also expected to define the window shape.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. The length of the frequencies is one point longer than the frame length, such that the frequencies delimits the spectrogram in frequency. If the frame length is odd, the negative Nyquist frequency is in the beginning of the array; if even, the negative Nyquist frequency is the mean of the beginning two elements in the array.
  • times: 1d numpy array
    A real array of uniformaly distanced time grids in units of second. It is relative to the beginning of the data block. The length of the times is larger than the actual n_frame by one, such that the times delimits the spectrogram in time.
  • spectrogram: 2d numpy array
    A real array of the periodogram estimate of the PSD in units of V²/kHz. The array is re-ordered in frequency according to the correspondence with the frequencies, and stacked consecutively by frames.
  • n_dof: int
    Degrees of freedom. By definition, always return 2.

multitaper_2d

multitaper_2d(window_length=1000, n_frame=100, n_offset=0, padding_ratio=2, half_bandwidth=3, n_taper=4)

Estimate the PSDs of the segments of an evenly sectioned, and optionally, zero-padded time series with the basic multitaper method.

• arguments:
  • window_length: int, default 1000
    Length of the series segment for one FFT frame.
  • n_frame: int, default 100
    Number of the FFT frames. It should, in principle, be positive. However, the actual n_frame may depend on the available IQ pairs left in the data block after deducting the n_offset. If n_frame is negative, all the available IQ pairs, after trimming off the incomplete segment, will be selected.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original segment length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the segment length.
  • half_bandwidth: float, default 3
    Half of the bandwidth of the DPSS tapers in units of frequency grid.
  • n_taper: int, default 4
    Number of DPSS tapers. To avoid spectral leakage, it is advisalbe to keep n_taper less than 2*half_bandwidth.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. The length of the frequencies is one point longer than the frame length, such that the frequencies delimits the spectrogram in frequency. If the frame length is odd, the negative Nyquist frequency is in the beginning of the array; if even, the negative Nyquist frequency is the mean of the beginning two elements in the array.
  • times: 1d numpy array
    A real array of uniformaly distanced time grids in units of second. It is relative to the beginning of the data block. The length of the times is larger than the actual n_frame by one, such that the times delimits the spectrogram in time.
  • spectrogram: 2d numpy array
    A real array of the basic multitaper estimate of the PSD in units of V²/kHz. The array is re-ordered in frequency according to the correspondence with the frequencies, and stacked consecutively by frames.
  • n_dof: int
    Degrees of freedom. By definition, return 2*n_taper.

adaptive_multitaper_2d

adaptive_multitaper_2d(window_length=1000, n_frame=100, n_offset=0, padding_ratio=2, half_bandwidth=3, n_taper=4)

Estimate the PSDs of the segments of an evenly sectioned, and optionally, zero-padded time series with the adaptive multitaper method.

• arguments:
  • window_length: int, default 1000
    Length of the series segment for one FFT frame.
  • n_frame: int, default 100
    Number of the FFT frames. It should, in principle, be positive. However, the actual n_frame may depend on the available IQ pairs left in the data block after deducting the n_offset. If n_frame is negative, all the available IQ pairs, after trimming off the incomplete segment, will be selected.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original segment length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the segment length.
  • half_bandwidth: float, default 3
    Half of the bandwidth of the DPSS tapers in units of frequency grid.
  • n_taper: int, default 4
    Number of DPSS tapers. To avoid spectral leakage, it is advisalbe to keep n_taper less than 2*half_bandwidth.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. The length of the frequencies is one point longer than the frame length, such that the frequencies delimits the spectrogram in frequency. If the frame length is odd, the negative Nyquist frequency is in the beginning of the array; if even, the negative Nyquist frequency is the mean of the beginning two elements in the array.
  • times: 1d numpy array
    A real array of uniformaly distanced time grids in units of second. It is relative to the beginning of the data block. The length of the times is larger than the actual n_frame by one, such that the times delimits the spectrogram in time.
  • spectrogram: 2d numpy array
    A real array of the adaptive multitaper estimate of the PSD in units of V²/kHz. The array is re-ordered in frequency according to the correspondence with the frequencies, and stacked consecutively by frames.
  • n_dof: 2d numpy array
    A real array of degrees of freedom whose elements correspond to the SDE in the spectrogram.

time_average_2d

time_average_2d(window_length=1000, n_frame=100, n_offset=0, padding_ratio=2, n_average=10, estimator='p', **kwargs)

Estimate the PSDs of the segments of an evenly sectioned time series by further sectioning each segment into sub-segments and averaging their SDEs with, optionally, zero-padding. The SDE can be computed with the periodogram, the basic multitaper, or the adaptive multitaper method. Usually the Bartlett's averaging scheme is adopted, except for the periodogram estimate with windows other than the boxcar. In such cases, the Welch's averaging scheme with half overlapping portions is used instead. For the adaptive multitaper method, the average is weighted by the degrees of freedom of each SDE.

• arguments:
  • window_length: int, default 1000
    Length of the sub-segment for one FFT frame.
  • n_frame: int, default 100
    Number of the FFT frames after average. It should, in principle, be positive. However, the actual n_frame may depend on the available IQ pairs left in the data block after deducting the n_offset. If n_frame is negative, all the available IQ pairs, after trimming off the insufficient sub-segments for one average, will be selected.
  • n_offset: int, default 0
    Offset from the beginning of the data block in units of IQ pairs.
  • padding_ratio: float, default 2
    Ratio of the zero-padded frame length to the original sub-segment length. The actual frame length will, if not radix-2 compatible, be inflated to the next power of 2. Mostly, the padding_ratio should be no less than 1. Otherwise, the zero-padding is suppressed, such that the frame length equals exactly to the sub-segment length.
  • n_average: int, default 10
    Number of sub-segments for one average. For Welch's averaging scheme, the actual number is 2*n_average-1 by counting the overlapped sub-segments.
  • estimator: str, default 'p'
    Estimating method to be adopted for the SDE. Valid options are ['p', 'm', 'a'], which stand for Periodogram, basic Multitaper, and Adaptive multitaper, respectively.
  • **kwargs:
    Some complementary arguments to be passed on to the backstage method depending on the estimator. If estimator='p', window, and possibly, beta should be specified. If estimator='m' | 'a', half_bandwidth and n_taper should be specified.
• returns:
  • frequencies: 1d numpy array
    A real array of uniformly distanced frequency grids in units of kHz. It is calculated based on the sampling rate, ranged between the negative Nyquist frequency and the positive one, and sorted in an ascending order. The length of the frequencies is one point longer than the frame length, such that the frequencies delimits the spectrogram in frequency. If the frame length is odd, the negative Nyquist frequency is in the beginning of the array; if even, the negative Nyquist frequency is the mean of the beginning two elements in the array.
  • times: 1d numpy array
    A real array of uniformaly distanced time grids in units of second. It is relative to the beginning of the data block. The length of the times is larger than the actual n_frame by one, such that the times delimits the spectrogram in time.
  • spectrogram: 2d numpy array
    A real array of averaged SDEs in units of V²/kHz. The array is re-ordered in frequency according to the correspondence with the frequencies, and stacked consecutively by averaged frames.
  • n_dof: int | 2d numpy array
    Degrees of freedom depending on the estimator. If estimator='p' | 'm', n_dof is a real number. If estimator='a', n_dof is a real array whose elements correspond to the SDEs in the spectrogram.

confidence_band

confidence_band(sde, level, n_dof)

Compute the confidence band of the SDE.

• arguments:
  • sde: float | numpy array
    Spectral density estimate. It can be a real number or a real array.
  • level: float
    Confidence level. It should be positive and less than 1.
  • n_dof: float | numpy array
    Degrees of freedom. It can be a real number or a real array, as long as it is of the same shape as the sde.
• returns:
  • upper_bound: float | numpy array
    The upper bound of the confidence band. It is of the same shape as the sde.
  • lower_bound: float | numpy array
    The lower bound of the confidence band. It is of the same shape as the sde.

plot_1d

plot_1d(frequencies, spectrum)

Plot the one dimensional SDE.

• arguments:
  • frequencies: 1d numpy array
    A real array of frequency grids in units of kHz.
  • spectrum: 1d numpy array
    A real array of the SDE in units of V²/kHz.
• returns:
  • None

plot_2d

plot_2d(frequencies, times, spectrogram)

Plot the two dimensional SDEs on a color-coded map.

• arguments:
  • frequencies: 1d numpy array
    A real array of frequency grids in units of kHz.
  • times: 1d numpy array
    A real array of time grids in units of second.
  • spectrogram: 2d numpy array
    A real array of SDEs in units of V²/kHz.
• returns:
  • None

References

1. Wiener, N. "Generalized harmonic analysis". Acta Math. 55, 117 (1930).
2. Khintchine, A. "Korrelationstheorie der stationären stochastischen Prozesse". Math. Ann. 109, 604 (1934).
3. Schuster, A. "The periodogram of magnetic declination as obtained from the records of the Greenwich Observatory during the years 1871—1895". Cambridge, Trans. Phil. Soc., 18, 107 (1899).
4. Bartlett, M. S. "Smoothing Periodograms from Time-Series with Continuous Spectra". Nature 161, 686 (1948).
5. Welch, P. "The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms". IEEE Trans. Audio Electroacoust. 15, 7073 (1967).
6. Thomson, D. J. "Spectrum estimation and harmonic analysis". Proc. IEEE 70, 1055 (1982).
7. Percival, D. B. & Walden, A. T. "Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques". Cambridge University Press, (1993).