Application of DS-DFT to the Fine Spectral Estimation of High-Noise Signals

Abstract

This paper presented an extended double-subsegment discrete Fourier transform (DS-DFT) algorithm as a tool for the fine spectral estimation of high-noise environments, which was previously effective in low-noise scenarios, and as such, its application to the analysis of noisy signals observed by a satellite-based interferometer was investigated. The observation of satellite-borne Doppler asymmetric spatial heterodyne spectroscopy (DASH) was first simulated to obtain the signals of low- and high-noise levels; then, a practical criterion to classify noise levels in the interferograms based on the DS-DFT results was introduced and validated by calculating the SNR. For high-noise signals, DS-DFT remains robust by employing phase differences and amplitude ratios for fine frequency estimation.

1. Introduction

Fourier transform technology is widely used in satellite remote sensing and data processing. It is a very effective method to obtain atmospheric information. For example, the global high-resolution thermospheric imaging Michelson interferometer (MIGHTI), installed on the ICON (Ionospheric Connection) Explorer Mission satellite launched by NASA in October 2019 [1], can be regarded as a combination of a stepping Michelson interferometer and an asymmetrical spatial heterodyne spectrometer (SHS). It can measure horizontal wind speed and temperature in the daytime and nighttime thermosphere (90–300 km).

MIGHTI measures the wind field by observing the Doppler frequency shift of visible to near-infrared radiation, which is emitted by particles moving within a neutral atmosphere. In its limb mode, the observation path passes through the atmosphere at different heights, and the MIGHTI-measured interferogram is essentially the integrated signal of the atmospheric radiation and absorption along the path. An inversion method based on an “onion peeling” algorithm has been established to extract the information of the atmospheric wind field at different heights from the observed or simulated interferogram [2,3,4]. In the inversion process (referred to as the Englert–Harding algorithm hereafter), the key issue is to obtain the atmospheric radiation spectrum from the observed interferogram and accurately estimate its Doppler frequency shift. A discrete Fourier transform (DFT) is often used to achieve this. However, there are two important issues in the Doppler-shifted frequency estimation:

(1)

DFT alone cannot accurately estimate a small Doppler shift frequency. It is well known that DFT spectral resolution depends on the sampling density of a signal sequence. For a single-frequency signal sequence of N sampling points, the normalized spectral resolution is 1/N, and the DFT-derived signal frequency is represented by a set of integers $k_p$ multiplying 1/N. If the signal frequency $f_0$ happens to have a small fraction, $-0.5\le \delta \le 0.5$, the DFT algorithm alone cannot accurately estimate the fractional part $\delta$. In the MIGHTI case, the fractional part can be introduced by a small Doppler shift related to the velocity of the light source. For example, the Doppler shift phases caused by a wind speed of 1 m/s are 2.10 mrad and 1.80 mrad (corresponding to binned horizontal pixels of 6.68 × 10⁻³ and 2.48 × 10⁻³) for MIGHTI-observed 557 nm and 630 nm emission lines, respectively.

(2)

Noise contamination influences the anti-noise robustness and estimation accuracy of a frequency estimator. In addition to various random noises, the interferograms obtained over a wide range of limb brightness apparently have an unexpected, systematic shift in fringe phase depending on the signal strength on the detector, which is determined to be characteristic of the instrument rather than geophysical. This shift is strongest for the dimmest signals and tends toward a zero shift at the brightest signals. In order to improve the accuracy of the estimated frequency and to obtain the fractional part $\delta$, two steps are usually adopted: rough estimation (estimation of integer $k_p$) and fine estimation (estimation of frequency offset $\delta$). The rough estimation can be determined from the frequency corresponding to the peak amplitude of the N-point DFT spectrum. There are two kinds of methods for fine estimation: direct and iterative [5]. The direct method uses the intensity values of the DFT frequency grids near the peak values to design the interpolation functions. Different interpolation schemes produce different direct estimation algorithms. The Englert–Harding algorithm can be considered as a direct estimation algorithm [3,4]. The calculation scheme of direct estimation is relatively simple, but in general, its accuracy is highly correlated with the frequency offset $\delta$. The iterative method can reduce the dependence of the accuracy on the frequency offset $\delta$. The iterations are performed by taking the deterministic estimate of the previous iteration as a rough estimate of the next iteration until convergence reaches [6,7,8]. Although in most cases only two iterations can achieve convergence to meet the accuracy requirements, the calculation is relatively complicated. In addition to the high computational complexity of the iteration method, it also cannot achieve phase information retrieval. The binary Fourier transform proposed by Huang and Xia [9] is a direct algorithm with a relatively simple calculation process but a high precision of frequency estimation results. The basic idea is to divide the input N-point sample into two subsegments with an equal length of data points, M = N/2 and perform a DFT calculation on them, respectively. From the phase difference of their DFT spectral amplitude peaks, the signal frequency $f_0$ and its fraction $\delta$ are estimated according to the shift property of the Fourier transform. But their algorithm only considers low-noise cases.

The purpose of this paper is to extend the binary Fourier transform algorithm of Huang and Xia [9] to high-noise cases, in order to quickly and accurately retrieve the frequency and phase information of a single tone from noisy interferograms observed by satellite-borne Doppler asymmetric spatial heterodyne spectroscopy (DASH), which is similar to the MIGHTI instrument. The DASH interferograms are simulated (Section 2) and used as input to develop an algorithm for an accurate spectral estimation of the signal with high-noise levels (Section 3). The simulation process flowchart and parameters are given in Appendix A, and the detailed derivation of the DS-DFT algorithm for the high-noise case is in Appendix B. On this basis, an operational criterion (Section 4) is proposed to determine the noise level and compare it with the signal-to-noise ratio calculated directly from the analog signals to verify its usability. Our discussions and conclusions are presented in Section 5.

2. Data of High Noise Level

In this study, DASH interferograms of the atmospheric airglow radiation observed from a limb-viewing satellite were used as the input signal. A detailed discussion of the instrument structure, measurement principles, and characteristics of the DASH-observed interferograms can be found in the references [2,3,4]. Here, only an outline is presented to show how to generate input signal data for this study.

2.1. Basic Equations of the Spatial Heterodyne Spectrometer

The DASH interferogram can be described as follows [3]:

$$I(x)={\displaystyle {\int }_{-\infty }^{\infty }A(\kappa ,x)}\beta (\kappa )\{{\mathrm{e}}^{i[\kappa x+\Phi (x,\kappa )]}+{\mathrm{e}}^{-i[\kappa x+\Phi (x,\kappa )]}\}d\kappa +\epsilon$$

(1)

where $I\left(x\right)$is the airglow radiation intensity recorded by each pixel bin in a row of the CCD (charge coupled device) detector (corresponding to one height level of the atmosphere), and $x$ is the position of the pixel bin in the row. Suppose $x=0$ is the midpoint of the row, and the positive integer $n=256$ is the number of pixel bins from the midpoint, each pixel bin contains two CCD pixels of the size $\mathsf{\Delta }x=13.5 \mathsf{\mu }\mathrm{m}$. $\mathsf{\kappa }$ is the heterodyne fringe frequency (Fizeau frequency) corresponding to the airglow emission line of wave number $\sigma$ and is defined by the Littrow wave number ${\sigma }_0$ and grating angle of the spectrometer ${\theta }_L$ as $\mathsf{\kappa }=4\left(\sigma -{\sigma }_0\right)\tan {\theta }_L$. The fringe frequency $\kappa =0$ corresponds to the Littrow wave number. In Appendix A, the flowchart of the simulation and inversion processes is presented in Figure A1, and the related simulation parameters are listed in

Table A1. Parameters for the SHS interferogram simulation.

Heterodyne fringe frequency (Fizeau frequency) $\kappa =4\left(\sigma -{\sigma }_0\right)\tan {\theta }_L$$({\mathrm{c}\mathrm{m}}^{-1}$)	$\mathsf{\sigma }$$({\mathrm{c}\mathrm{m}}^{-1}$): Wave number ${10}^7/557.735 \mathrm{n}\mathrm{m}$ ${10}^7/630.031 \mathrm{n}\mathrm{m}$ ${\sigma }_0 ({\mathrm{c}\mathrm{m}}^{-1}$): Littrow wave number ${10}^7/554.670 \mathrm{n}\mathrm{m}$ ${10}^7/633.910 \mathrm{n}\mathrm{m}$ ${\theta }_L$ = 8.2°: Littrow angle
Spectral resolution of interferogram $\mathrm{d}\mathsf{\sigma }={\displaystyle \frac{1}{8·x_{\mathit{max}}·\tan {\theta }_L}}$$({\mathrm{c}\mathrm{m}}^{-1}$)	$x$ (cm): Sample position $dx$ = 13.5 (μm): Sample interval (CCD pixel size) $dx={\displaystyle \frac{1}{4N·\mathrm{d}\mathsf{\sigma }·\tan {\theta }_L}}$ N = 150: Sample number (column no. of a row) Sample Range: [−$x_{max},$+$x_{max}$] $x_{max}={\displaystyle \frac{N·dx}{2}}$(μm) or $x_{max}={\displaystyle \frac{1}{2}}W\cos {\theta }_L$ Width of grating (cm): $W={\displaystyle \frac{2x_{max}}{\cos {\theta }_L}}$
Maximum wave number ${\sigma }_{max} ({\mathrm{c}\mathrm{m}}^{-1}$) Minimum wave number ${\sigma }_{min} ({\mathrm{c}\mathrm{m}}^{-1}$)	${\sigma }_{max}={\sigma }_0$ $W={\displaystyle \frac{N}{8\left({\sigma }_{max}-{\sigma }_{min}\right)\sin {\theta }_L}}$

and

Table A2. Parameters for the noise simulation.

Parameter Name	Specifications
Readout noise ${\sigma }_{read}={\displaystyle \frac{2C_{avg}}{A_a·\mathsf{\Omega }·d\sigma ·D^{*}·{\tau }_{tot}·\sqrt{A_{d}T}}}$ ${\sigma }_{read}=4.2$ electrons/pixel/read out	$C_{avg}=1.5259$: Average $C\left(\kappa \right)$ $A_a$ = $60\times 21 {\mathrm{m}\mathrm{m}}^2:$ Aperture area $\mathsf{\Omega }=\pi /\left(4f^2+1\right)=$ 0.1307 (rad): FOV of SHS f = 2.4: F number $D^{*}$ = ${2\times 10}^{10}(\mathrm{W}/\mathrm{c}\mathrm{m}/\sqrt{\mathrm{H}\mathrm{z}}$): Detection rate $A_d=13.5\times 13.5({\mathsf{\mu }\mathrm{m}}^2$): Area of CCD pixel T = 114.9 (frames/s): Frame rate $d\sigma$$({\mathrm{c}\mathrm{m}}^{-1}$): Spectral resolution of CCD detector
Conversion factor of spectral radiance to electron numbers $C\left(\kappa \right)={\displaystyle \frac{A_d·\mathsf{\Omega }·t_{int}·QE}{(hc/\sigma )·d\sigma }}\sqrt{{\displaystyle \frac{N}{2}}}$	$A_d,$$\mathsf{\Omega }$, $d\sigma$: See above $t_{int}=0.5\left(\mathrm{m}\mathrm{s}\right)$: Integration time $N$ = 150: Sample number of interferogram $\sigma \left({\mathrm{c}\mathrm{m}}^{-1}\right):$ Doppler-shifted wave numbers $h$ Planck constant $c$ Light speed $QE={0.70\times e}^{-f\left(\sigma \right)}$: Quantum efficiency $f\left(\sigma \right)={\displaystyle \frac{1-\sigma }{2{\sigma }_{max}}}$
Total transmission function of SHS: ${\tau }_{tot}={\tau }_1·{\tau }_2·{\displaystyle \frac{{\eta }_A+{\eta }_B}{2}}$	Transmission function of entrance and exit optics ${\tau }_1={0.85\times e}^{-f_1\left(\sigma \right)}$  $f_1={\displaystyle \frac{{\sigma +\sigma }_{max}-{\sigma }_{min}}{{\left(2\ast {\sigma }_{max}-{\sigma }_{min}\right)}^2}}$ ${\tau }_2={0.80\times e}^{-f_2\left(\sigma \right)}$  $f_2={\displaystyle \frac{{\sigma +\sigma }_{max}-{\sigma }_{min}}{{\left(8\ast {\sigma }_{max}-{\sigma }_{min}\right)}^2}}$ Transmission function of grating: ${\eta }_A=0.80+0.01\ast \cos \left(0.2\ast \pi \ast \sigma \right)$ ${\eta }_B=0.85+0.01\ast \cos \left(0.2\ast \pi \ast \sigma \right)$
Quantization noise ${\sigma }_{quant}={\displaystyle \frac{L_{max}}{\sqrt{12}\left(2^M-1\right)}}{C}_{avg}$	M = 17: CCD bit number $L_{max}=\mathrm{s}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \mathrm{l}\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{l}/C_{avg}$: Maximum radiance detected by CCD, converted from its saturation level $={32\times 10}^6$(electrons)
Dark current ${\sigma }_{dark}=0.02 \mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{n}\mathrm{s}/\mathrm{s}/\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}$	Each physical pixel has an individual dark current simulated by a random number generator

.

In Equation (1), $\beta \left(\kappa \right)$ and $A(\kappa ,x)$ is the spectral density and the line shape function of the line, respectively, and $\mathsf{\Phi }(x,\kappa )$ is the known phase determined by the instrument parameters. $A(\kappa ,x)\beta \left(\kappa \right)$ can be expressed as the integral of the line of sight (LOS) corresponding to the CCD row (through different atmospheric altitude layers) [10]:

$$A(\kappa ,x)\beta (\kappa )=\frac{1}{2}{\displaystyle {\int }_0^{\infty }L(\kappa ,h)}{B}_0(\kappa ,h)\mathrm{d}s$$

(2)

where $B_0(\kappa ,h)$ is the spectral density (volume emissivity) of the Fizeau wave number $\kappa$ at height h, and $L(\kappa ,h)$ is the line function of the spectral line at height h, which is assumed to be a Gaussian function. Note that the wave number $\sigma$ consists of two parts: the intrinsic wave number of the airglow radiation for the emitting particles at rest (corresponding Fizeau frequency) and the Doppler frequency shift for the particles moving with the neutral atmosphere (corresponding Fizeau frequency). Therefore, the atmospheric wind speed at different altitudes h is different, thus the Doppler frequency shift is a function of altitude.

In a realistic measurement, the observed interferogram also has measurement noise $\epsilon$, including instrument measurement noise and background radiation noise. In our simulation, the background radiation noise is not considered. The instrument noise only includes readout noise ${\sigma }_{read}$, quantization noise ${\sigma }_{quant}$, and dark current ${\sigma }_{dark}$. The simulation parameters of the three noises are listed in Table A2 in Appendix A. The total noise standard deviation can be obtained as follows:

$${\sigma }_{tot}=\sqrt{{\sigma }_{read}^2+{\sigma }_{quant}^2+{\sigma }_{dark}^2}$$

(3)

2.2. Simulated DASH Interferogram

According to Equations (1) and (2) and the simulation parameters given in Appendix A, the DASH interferogram was simulated and is presented in Figure 1. In the simulation, the 557.7 nm airglow radiation was considered. Its spectral density $B_0(\kappa ,h)$ was taken from the VER (volume emission rate) model [11], and the atmospheric temperature and horizontal wind were from the MSIS (mass-spectrometer-incoherent-scatter) model [12] and the HWM (horizontal wind) model [13]. The three instrumental measurement noises were contained in the simulated interferogram.

The left and right panels of Figure 1 correspond to the two fields of view (FOVs), whose center lines of sight were directed at an angle of 45° and 135°, respectively, with respect to the satellite motion direction. The horizontal direction of each FOV (along the satellite orbit, in the CCD column) is the $\mathrm{x}$ direction of the interferogram with a resolution of two pixels. The vertical direction of each FOV (the CCD row) is the altitude of the earth’s atmosphere, and at the center of each row $\mathrm{x}=0$. In the vertical direction, the pixel bin used for the 557.7 nm observation included 16 pixels, corresponding to an altitude resolution of about 2 km. The observed altitude area is about 80–120 km. The upper limit of the LOS integral is 120 km, above which the radiation intensity is close to zero and no simulation is performed. The atmosphere layers of different altitudes seen by the DASH were projected onto different CCD rows. Fizeau interference fringes of the same height appear simultaneously on the same CCD row. In the image, the intensity variation between pixels is characterized by a vertical stripe pattern. This feature is mainly controlled by the vertical gradient of the airglow emission, as well as the altitude variations in the satellite and earth rotations, and the atmospheric wind and temperature fields. In the simulations, the atmosphere was assumed to be horizontally uniform, and changes in the instrument visibility and phase between pixels in a row are ignored. However, the intensity fluctuation between pixels in a given row is significantly larger. This is because different pixel bins in the same row actually see Fizeau interference fringes with different phases. Despite the phase modulation, approximately 19 intensity peaks in each row can be seen, corresponding to the Fizeau fringe frequency ${\mathsf{\kappa }}_0=13.773 {\mathrm{c}\mathrm{m}}^{-1}$ for the 557.7 nm line. The height distribution of the interference fringe intensity is similar to that of the volume emission rate given by the VER model (not shown). The peak values of the radiation intensity are distributed in the narrow height region of 95 ± 5 km, and the maximum value is about 5.25 KR. Also, it should be noted that the vertical coordinates of the two FOVs in Figure 1 are different. The same row of CCD in the two FOVs observes different atmospheric heights, with a difference of about 2 km. This difference is caused by the curvature effect of the spherical atmosphere in the different viewing regions of the two FOVs.

2.3. Fourier Transform of the DASH Interferogram

The spaceborne DASH measurements aim to obtain the height profile of the atmospheric wind speed from the observed interferogram $I\left(x\right)$ in Equations (1) and (2). This is achieved by solving or inverting the integral equations. The first step of the Englert–Harding inversion algorithm is to obtain the spectrum of the airglow emission through the Fourier transform (FFT) of the observed $I\left(x\right)$. If the interferogram sequence $I\left(x\right)$ is assumed to be a single frequency ${\kappa }_0$, its FFT spectrum $F\left[I\left(x\right)\right]$ has two peaks of amplitude around $\pm {\kappa }_0$. This spectrum is called a bilateral spectrum, and the corresponding input interferogram is called a bilateral interferogram. A single-side spectrum with only one peak amplitude at $+{\kappa }_0$ can be obtained by applying the δ-function $\mathsf{\delta }\left(\kappa -{\kappa }_0\right)$ to the bilateral spectrum $F\left[I\left(x\right)\right]$, i.e., $F_{\delta }=F\left[I\left(x\right)\right]·\delta (\kappa -{\kappa }_0)$, and a single-side interferogram corresponding to a single frequency $+{\kappa }_0$ can be obtained by an inverse Fourier transform $F^{-1}{(F}_{\delta })$ [3]. The δ-function was appropriated by a Gaussian function with a suitable half-width in this study.

Figure 2 and Figure 3 show the interferograms and spectra for the simulated observations of 557 nm airglow by the DASH FOV1 at 111.9 km. For comparison, Figure 2 does include any measurement noise, i.e., $\epsilon =0$, but Figure 3 contains the readout noise, quantization noise, and dark current. Both figures display the double-sided interferogram $I\left(x\right)$ (upper left panel) simulated from Equations (1) and (2), the double-sided spectrum generated by the Fourier transform $F\left[I\left(x\right)\right]$ (middle- and lower-left panels for the amplitude and phase spectra, respectively), the single-sided spectrum selected with the δ-function $F_{\delta }=F\left[I\left(x\right)\right]·\delta (\kappa -{\kappa }_0)$ (middle- and lower-right panels for the amplitude and phase spectra, respectively), and the single-sided interferogram $F^{-1}{(F}_{\delta })$ (upper-right panel) derived from the inverse Fourier transform. The simulated (or measured) double-sided interferograms $I\left(x\right)$ (black lines) in both figures were bias-removed and apodized with the Haming window (red lines) before the Fourier transform. The $I\left(x\right)$ in Figure 2 and Figure 3 were collected from a row of CCD bins in the simulated FOV1 interferogram shown in Figure 1.

In the absence of instrumental noise (Figure 2), the bilateral spectrum $F\left[I\left(x\right)\right]$ has two symmetric amplitude peaks about the origin (κ = 0), while the amplitudes at the other frequencies (wave number) are close to zero (middle-left panel), indicating the absence of other frequency components or white noise characterized by a continuous spectrum. The phase angles of the spectrum are expressed in continuous arguments and show a piecewise linear change with the Fizeau frequency κ (lower left panel). The Gaussian function was applied to the bilateral spectrum $F\left[I\left(x\right)\right]$, and the unilateral spectrum $F_{\mathsf{\delta }}$ (middle-right panel) was obtained with only one amplitude peak and the same phase angles as those of the bilateral spectrum (lower-right panel), since the Gaussian function has the phase response of a zero-phase shift. The unilateral interferogram $F^{-1}{(F}_{\delta })$ obtained from the inverse Fourier transform of the unilateral spectrum is complex, with real (black lines) and imaginary (red lines) parts (top-right panel). Due to the apodization effect of the Haming window, the unilateral interferogram and apodized bilateral interferogram rapidly decay at both ends of the sampling window.

In the presence of the measurement noises (Figure 3), the bilateral and unilateral interferograms, as well as their amplitude and phase spectra, are completely different from those of the noise-free case. The amplitude spectrum of $F\left[I\left(x\right)\right]$ (middle-left panel) generally shows frequency components of non-zero amplitude from a continuous spectrum of white noise, although there are still two symmetric peaks about the origin (κ = 0). The phases of the bilateral spectrum (lower-left panel) oscillate strongly with the Fizeau frequency κ due to the noises but still show a linear trend approximately. The unilateral spectrum ${\mathrm{F}}_{\delta }$ selected by the Gaussian function from the bilateral spectrum (middle-right panel) has only one peak amplitude but shows approximately symmetric frequency components distributed around the central frequency (wave number) of the peak and extended to zero frequency. Due to the phase response of the zero-phase shift for the Gaussian function, the unilateral phase spectrum (lower-right panel) is the same as the bilateral one, preserving the phase oscillations caused by the noises. The unilateral interferogram $F^{-1}{(F}_{\delta })$ (upper-right panel) obtained from the inverse Fourier transform of the unilateral spectrum is also significantly different from that of the noise-free case.

2.4. Amplitudes and Phases of the Unilateral Spectrum

According to Equations (1) and (2), the amplitude $A\left({\kappa }_0,x\right)\beta \left({\kappa }_0\right)$ and total phase ${\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ can be obtained from the unilateral interferogram of the Fizeau frequency $+{\kappa }_0$ and expressed as [3]:

$$\begin{array}{c}\hfill {[A({\kappa }_0,x)\beta ({\kappa }_0)]}^2={\{\mathrm{Re}[I_{\kappa 0}(x)]\}}^2+{\{\mathrm{Im}[I_{\kappa 0}(x)]\}}^2\\ \hfill {\Phi }_{tot}(x,{\kappa }_0)={\kappa }_{0}x+\Phi (x,{\kappa }_0)={\tan }^{-1}\left\{{\displaystyle \frac{\mathrm{Im}[I_{\kappa 0}(x)]}{\mathrm{Re}[I_{\kappa 0}(x)]}}\right\}\end{array}$$

(4)

Here, for the sake of simplicity, the unilateral interferogram $F^{-1}{(F}_{\delta })$ is denoted as $I_{\kappa 0}\left(x\right)$. According to Equation (4), the difference between the total phase ${\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ at any x and the total phase ${\mathsf{\Phi }}_{tot}\left(0,{\kappa }_0\right)$ at x = 0, i.e., $∆{\mathsf{\Phi }}_{tot}\left(x\right)={\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)-{\mathsf{\Phi }}_{tot}\left(0,{\kappa }_0\right)$, should be equal to the phase $2\pi {\kappa }_{0}x$ generated by the Fizeau frequency (where the phase angle is expressed in rad). However, in the presence of the Doppler shift in the Fizeau frequency induced by the atmospheric wind, the phase difference between $∆{\mathsf{\Phi }}_{tot}\left(x\right)$ and $2\mathsf{\pi }{\kappa }_{0}x$, i.e., $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x\right)=∆{\mathsf{\Phi }}_{tot}\left(x\right)-2\mathsf{\pi }{\kappa }_{0}x$, is not zero. Thus, the Englert–Harding algorithm focuses on the behavior of

$$\delta \begin{array}{cc}{\Phi }_{tot}(\mathrm{x},{\kappa }_0)=[{\Phi }_{tot}(x,{\kappa }_0)-{\Phi }_{tot}(0,{\kappa }_0)]-2\pi {\kappa }_{0}x={\phi }_0({\kappa }_0)+D_{\kappa }({\kappa }_0)\cdot x& \end{array}$$

(5)

where the phase difference $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x\right)$is assumed to be a linear function of $x$. The slope $D_{\kappa }\left({\kappa }_0\right)$ is related to the Doppler shift in the Fizeau frequency due to the atmospheric wind motion, and the intercept ${\phi }_0\left({\kappa }_0\right)$ is introduced to account for phase errors in the calibration of the zero-wind phase, phase distortion, and various noise sources in the optical system. Ideally, in the absence of these calibration errors, ${\phi }_0\left({\kappa }_0\right)$ should be zero, and Figure 4 and Figure 5 show the amplitudes $⌈A\left({\kappa }_0,x\right)\beta \left({\kappa }_0\right)⌉$ and phase differences $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ computed by Equations (3) and (5) for the unilateral interferograms without and with the measurement noises (Figure 2 and Figure 3), respectively.

In the noise-free case (Figure 4), the amplitudes $⌈A\left({\kappa }_0,x\right)\beta \left({\kappa }_0\right)⌉$ (left panel) are seen to decay rapidly with the increase in absolute x values due to the apodization effect of the Haming window. Since there is no noise, the amplitudes vary smoothly with x. Note that the amplitudes are asymmetric with respect to the origin x = 0 and do not reach their maximum at x = 0. This is because the phase ${\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ in Equation (1) was assumed not to be zero at x = 0 but includes the phase distortion term $\mathsf{\Theta }\left(x,\kappa \right)=4\pi \sigma \left(\mathsf{\Delta }d\right)$. This term is introduced by the asymmetry of the two optical arms of the interferometer with their optical path difference $\mathsf{\Delta }d=2.5 \mathrm{c}\mathrm{m}$ and the wave number $\sigma$ of the airglow radiation. If the instrument phase distortion term is not considered, the phase ${\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)=0$ at $x=0$, and the total phase ${\mathsf{\Phi }}_{tot}\left(0,{\kappa }_0\right)=0$, the amplitudes should be symmetric with respect to the origin (x = 0), and their maximal value at x = 0.

For the noise-free case, the phase differences $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ in Equation (5) are presented in the right panel of Figure 4. The total phases ${\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ were computed from Equation (4) for the unilateral interferogram (top-right panel of Figure 2) and expressed by continuous arguments.${\mathsf{\Phi }}_{tot}\left(0,{\kappa }_0\right)$ is the total phase at x = 0, and $2\mathsf{\pi }{\kappa }_{0}x$ denotes the phases generated by the Fizeau frequency ${\kappa }_0$. It can be seen that the phase differences $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ (black lines) vary nearly linearly with x, except for both ends of the x interval, where the absolute values of $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ change rapidly due to the apodization effect of the Haming window. The phase differences $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ were fitted by a linear function $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)={\phi }_0\left({\kappa }_0\right)+D_{\kappa }\left({\kappa }_0\right)·x$, which is shown by the red dotted line. For the fitting, about 10% of data points at both ends of the x interval were omitted. The fitting coefficients are ${\phi }_0\left({\kappa }_0\right)=-1.487\times {10}^{-3} \mathrm{r}\mathrm{a}\mathrm{d}$ and $D_{\kappa }\left({\kappa }_0\right)=0.941081 {\mathrm{c}\mathrm{m}}^{-1}$. ${\phi }_0\left({\kappa }_0\right)$ is an additional term introduced to eliminate the calibration errors of the zero-wind phase, phase distortion, and noise effects. In the ideal state without calibration error and noise, ${\phi }_0\left({\kappa }_0\right)$ should be zero.$D_{\kappa }\left({\kappa }_0\right)$ is associated with the Doppler shift in the Fizeau frequency and is caused by the atmospheric wind motions. The mean of the Doppler phase $D_{\kappa }\left({\kappa }_0\right)·x=-8.0865\times {10}^{-4} \mathrm{r}\mathrm{a}\mathrm{d}$, and its standard deviation is $3.511\times {10}^{-1} \mathrm{r}\mathrm{a}\mathrm{d}$.

In the presence of the measurement noises (Figure 5), the amplitudes of the unilateral interferogram $⌈A\left({\kappa }_0,x\right)\beta \left({\kappa }_0\right)⌉$ (left panel) increase by a maximum of nearly 50% (60 W/m² vs. 40 W/m²) and show large oscillations with x, even with two almost-equal peaks. The phase differences $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ (right panel) also have significant oscillations due to the influence of the noises, though they still show an approximately linear change with x. These features are quite different from those of the noiseless case (left panel of Figure 4). For $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right),$the fitting coefficient ${\phi }_0\left({\kappa }_0\right)=0.023558 \mathrm{r}\mathrm{a}\mathrm{d},$ which was larger than the noise-free one ($-1.487\times {10}^{-3} \mathrm{r}\mathrm{a}\mathrm{d}$) by 0.025 rad, and $D_{\kappa }\left({\kappa }_0\right)=1.28930 {\mathrm{c}\mathrm{m}}^{-1}$, which increased by 37% in comparison with the noise-free case ($0.941081 {\mathrm{c}\mathrm{m}}^{-1}$). The mean $1.5996\times {10}^{-2} \mathrm{r}\mathrm{a}\mathrm{d}$ of the Doppler phase$D_{\kappa }\left({\kappa }_0\right)·x$ increased by almost two orders of magnitude, and the standard deviation $4.9251\times {10}^{-1} \mathrm{r}\mathrm{a}\mathrm{d}$ increased by 40% compared to the mean of $8.0865\times {10}^{-4} \mathrm{r}\mathrm{a}\mathrm{d}$ and the standard deviation of $3.511\times {10}^{-1} \mathrm{r}\mathrm{a}\mathrm{d}$ in the noise-free case, respectively.

The above results show that the existence of strong measurement noises introduces large deviations in the estimation of the LOS-integrated amplitudes $⌈A\left({\kappa }_0,x\right)\beta \left({\kappa }_0\right)⌉$ and the Doppler phase shift $D_{\kappa }\left({\kappa }_0\right)$.

3. Fourier Transform of Double Subsegments

The binary Fourier transform proposed by Huang and Xia [9] is a direct algorithm with a relatively simple calculation process but a high precision in frequency estimation results. In this section, the key points of Huang and Xia’s DS Fourier transform algorithm are summarized and extended to the case of high noise.

3.1. DFT of N-Point Samples and Their Double Subsegments

A wave signal $r\left[n\mathsf{\Delta }\right]$ of single frequency $f_0$, amplitude A, and phase angle ${\theta }_0$ is sampled at n points ($0\le n\le N-1$) with sampling interval Δ. The frequency $f_0$ is usually estimated by the DFT of the entire N-point sequence. In order to improve the accuracy of the estimated frequency and to obtain the fractional part $\delta$, Huang and Xia [9] proposed an algorithm to divide the whole sequence into two subsegments with the same length and perform DFT against each of them, respectively. The frequency $f_0$is estimated from the phase difference between the two spectral amplitude maxima.

Without loss of generality, N is assumed to be even, that is N = 2 M. The sample sequence $r\left[n\mathsf{\Delta }\right](0\le n\le 2M-1)$ is divided into two subsegments $\left\{r_0\left[n\mathsf{\Delta }\right]\right\}$ and $\left\{r_1\left[n\mathsf{\Delta }\right]\right\}$. The integer frequencies $k_0^{*}$ and $k_1^{*}$ are determined from the subsegments $\left\{r_0\left[n\mathsf{\Delta }\right]\right\}$ and $\left\{r_1\left[n\mathsf{\Delta }\right]\right\}$, respectively, and it is assumed $k_0^{*}\ne k_1^{*}$. The corresponding frequency offsets are ${\epsilon }_0$ and ${\epsilon }_1$ ($-0.5\le \left[{\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1\right]\le 0.5$). So, the signal frequency$f_0$ can be written as follows:

$$f_0={\displaystyle \frac{({\mathrm{k}}_0^{*}+{\epsilon }_0)}{M\Delta }}={\displaystyle \frac{({\mathrm{k}}_1^{*}+{\epsilon }_1)}{M\Delta }}\begin{array}{cc},& {\mathrm{k}}_0^{*},\end{array}{\mathrm{k}}_1^{*}\in Z^{+}$$

(6)

As shown in Appendix B in detail, the DFT spectra of the two subsegments $\left\{r_0\left[n\mathsf{\Delta }\right]\right\}$ and $\left\{r_1\left[n\mathsf{\Delta }\right]\right\}$ can be written as follows:

$$\begin{array}{l}{R}_0(k)=A{\displaystyle \frac{\sin [\pi {\epsilon }_0]}{\sin [\pi {\epsilon }_0/M]}}{e}^{j[{\theta }_0+(M-1){\epsilon }_0\pi /M]}\begin{array}{cc},& 0\le k\le M-1\end{array}\\ R_1(k)=A{\displaystyle \frac{\sin [\pi {\epsilon }_1]}{\sin [\pi {\epsilon }_1/M]}}{e}^{j[{\theta }_0+2\pi f_{0}M\Delta +(M-1){\epsilon }_1\pi /M]}\begin{array}{cc},& 0\le k\le M-1\end{array}\end{array}$$

(7)

Equations (6) and (7) are Equations (A9) and (A10) in Appendix B.

3.2. Phases of DS-DFT Spectra

The frequency offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ are associated with the phase terms of Equation (7). Let $k=k_0^{*}$ and $k=k_1^{*}$ be the frequencies determined from the amplitude peaks of the two DS-DFT spectra, respectively. As mentioned previously, the normalized frequencies $k_0$ and $k_1$ derived from the DFT analysis are integer numbers. The phase terms of Equation (7) can be written as

$$\begin{array}{l}{\Phi }_0[k_0^{*}]={\theta }_0+(M-1){\epsilon }_0\pi /M+2\pi k_0\\ {\Phi }_1[k_1^{*}]={\theta }_0+2\pi f_{0}M\Delta +(M-1){\epsilon }_1\pi /M+2\pi k_1\end{array}$$

(8)

where $k_0,k_1\in Z^{+}$. The phase angles ${\mathsf{\Phi }}_0\left[k_0^{*}\right]$ and ${\mathsf{\Phi }}_1\left[k_1^{*}\right]$ include a periodicity of $2\pi$. The phase difference between ${\mathsf{\Phi }}_0\left[k_0^{*}\right]$ and ${\mathsf{\Phi }}_1\left[k_1^{*}\right]$ can be expressed as

$$\begin{array}{l}\Delta \Phi ={\Phi }_1[k_1^{*}]-{\Phi }_0[k_0^{*}]\\ =2\pi f_{0}M\Delta +\left(M-1\right){\displaystyle \frac{{\epsilon }_1-{\epsilon }_0}{M}}\pi +2\pi \left(k_1-k_0\right)\\ =2\pi f_{0}M\Delta -\left(M-1\right){\displaystyle \frac{k_1^{*}-k_0^{*}}{M}}\pi +2\pi \left(k_1-k_0\right)\end{array}$$

(9)

By using Equation (6), the frequency offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ can be written as

$$\begin{array}{l}2\pi {\epsilon }_0=[\Delta \Phi +(M-1){\displaystyle \frac{\Delta k^{*}}{M}}\pi ]-2\pi k_0^{*}-2\pi (k_1-k_0)\\ 2\pi {\epsilon }_1=[\Delta \Phi +(M-1){\displaystyle \frac{\Delta k^{*}}{M}}\pi ]-2\pi k_1^{*}-2\pi (k_1-k_0)\end{array}$$

(10)

where $\mathsf{\Delta }{k}^{*}=k_1^{*}-k_0^{*}$. This term is introduced by our assumption that $k_0^{*}\ne k_1^{*}$ and is the only difference from Equation (9) of Huang and Xia [9], where the $\mathsf{\Delta }{k}^{*}$ term vanishes due to the assumption that $k_0^{*}=k_1^{*}$. The ambiguity of $2\pi (k_1-k_0)$ due to the multi-valued periodicity can be eliminated by two different approaches as discussed below.

3.2.1. Principal Value of Argument: DS-H Approach

The phase angles in Equation (10) are taken as their principal values to eliminate the multi-valued uncertainty and to uniquely solve the fractional frequencies ${\epsilon }_0$ and ${\epsilon }_1$ in the interval of [−π, 2π]. This approach was adopted by Huang and Xia [9] and is referred to as the DS-H approach hereafter. Based on this approach, the frequency estimates of $f_0$ can be obtained from Equations (6) and (10) and written as

$$\begin{array}{l}M\Delta *{\widehat{f}}_0=k_0^{*}+{\epsilon }_0=k_0^{*}+[\Delta \Phi /2\pi +(M-1)\frac{\Delta k^{*}}{2M}]\\ M\Delta *{\widehat{f}}_1=k_1^{*}+{\epsilon }_1=k_1^{*}+[\Delta \Phi /2\pi +(M-1)\frac{\Delta k^{*}}{2M}]\end{array}$$

(11)

where $\widehat{f_0}$ and $\widehat{f_1}$ are the frequency estimates of $f_0$ derived from the two DS-DFT spectra. They can be assumed to be equal, since the wave signal $r\left[n\mathsf{\Delta }\right]$ is assumed to have only a single frequency of $f_0={\displaystyle \frac{k_p+\delta }{N}}$with an integer $k_p$ and $-0.5\le \delta <0.5$. By combining the two estimates in Equation (11) with $\widehat{f_0}=\widehat{f_1}=\widehat{f}$, the frequency estimates can be expressed as

$$\begin{array}{l}M\Delta *{\widehat{f}}_1+M\Delta *{\widehat{f}}_0=k_0^{*}+k_1^{*}+2[\Delta \Phi /2\pi +(M-1)\frac{\Delta k^{*}}{2M}]\\ M\Delta *\widehat{f}=(k_0^{*}+k_1^{*})/2+[\Delta \Phi /2\pi +(M-1)\frac{\Delta k^{*}}{2M}]\end{array}$$

(12)

3.2.2. Continuous Argument: DS-W Approach

As discussed in Section 2, for the bilateral and unilateral spectra of DASH interferograms, the spectral phase angles generally vary linearly with the fringe frequency κ. This feature is based on the continuous expression of the arguments. Thus, there is no need to remove the multi-valued uncertainty of the phases. With continuous arguments, Equation (12) can be expressed as follows:

$$M\Delta *\widehat{f}=(k_0^{*}+k_1^{*})/2+[\Delta \Phi /2\pi +(M-1)\frac{\Delta k^{*}}{2M}]-(k_1-k_0)$$

(13)

This expression is referred to as the DS-W approach hereafter. It is possible to introduce phase jumps in the continuous arguments, if the difference between two adjacent angles exceeds 3.5 rad in absolute value. In this case, the two adjacent phases are re-processed to ensure that their absolute difference never exceeds 3.5 rad and to correct the phase angles defined in the interval [0, 2π] to a continuous fashion.

3.3. Amplitudes of DS-DFT Spectra

The frequency offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ are also associated with the amplitude terms of Equation (7), i.e.,

$$\begin{array}{l}{A}_0=A{\displaystyle \frac{\sin [\pi {\epsilon }_0]}{\sin [\pi {\epsilon }_0/M]}}\\ A_1=A{\displaystyle \frac{\sin [\pi {\epsilon }_1]}{\sin [\pi {\epsilon }_1/M]}}\end{array}$$

(14)

For simplicity, let $x=\pi {\epsilon }_0/M$ and $x=\pi {\epsilon }_0$, or $x=\pi {\epsilon }_1/M$ and $x=\pi {\epsilon }_1$. By using Taylor’s expansion $\mathrm{s}\mathrm{i}\mathrm{n}x=\sum _{n=0}^{\infty }{(-1)}^n{\displaystyle \frac{x^{2n+1}}{\left(2n+1\right)!}}$ and keeping only the terms of n = 0 and 1, Equation (14) can be approximated by

$$\begin{array}{l}{A}_0=AM[1-{(\pi {\epsilon }_0)}^2/6]\\ A_1=AM[1-{(\pi {\epsilon }_1)}^2/6]\end{array}$$

(15)

By taking the ratio $A_0/A_1$ from Equation (15) and neglecting the fourth-order term, an approximation can be obtained in the form of

$${(\pi {\epsilon }_0)}^2-{(\pi {\epsilon }_1)}^2=6(1-A_0/A_1)$$

(16)

3.4. Frequency Offsets of DS-DFT Spectra

As shown in Equation (7), the frequency offsets ${\epsilon }_0$ and ${\epsilon }_1$ are related to the phases and amplitudes of the DS-DFT spectra and are explicitly expressed by Equations (10) and (16) from which ${\epsilon }_0 \mathrm{and} {\epsilon }_1$ can be solved.

3.4.1. Frequency Offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ in the Case of ${\mathrm{k}}_0^{\mathrm{*}}={\mathrm{k}}_1^{\mathrm{*}}$

In the case of $k_0^{*}=k_1^{*}$, Equations (10) and (16) can be written as

$$\begin{array}{l}2\pi {\epsilon }_0+2\pi {\epsilon }_1=2\Delta \Phi -2\pi (k_1^{*}+k_0^{*})=a\\ {(2\pi {\epsilon }_0)}^2-{(2\pi {\epsilon }_1)}^2=24(1-A_0/A_1)=b\end{array}$$

(17)

It is not difficult to obtain

$$\begin{array}{c}2\mathsf{\pi }{\epsilon }_0={\displaystyle \frac{a+b/a}{2}}\\ 2\mathsf{\pi }{\epsilon }_1={\displaystyle \frac{a-b/a}{2}}\end{array}$$

(18)

3.4.2. Frequency Offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ in the Case of ${\mathrm{k}}_0^{\mathrm{*}}\ne {\mathrm{k}}_1^{\mathrm{*}}$

When $k_0^{*}\ne k_1^{*}$, Equations (10) and (16) can be written as

$$\begin{array}{l}2\pi {\epsilon }_1-2\pi {\epsilon }_0=-2\pi \Delta k^{*}=a\\ {(2\pi {\epsilon }_0)}^2-{(2\pi {\epsilon }_1)}^2=24(1-A_0/A_1)=b\end{array}$$

(19)

It is easy to obtain

$$\begin{array}{c}2\mathsf{\pi }{\epsilon }_0=-{\displaystyle \frac{a-b/a}{2}}\\ 2\mathsf{\pi }{\epsilon }_1={\displaystyle \frac{a+b/a}{2}}\end{array}$$

(20)

3.4.3. Frequency Offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\mathsf{\epsilon }}_1$ in the Case of ${\mathrm{A}}_0\ne {\mathrm{A}}_1$

As shown in Equations (18) and (20), the frequency offsets ${\epsilon }_0$ and ${\epsilon }_1$ are related to the coefficients $(a\pm b/a)$ in the cases of both $k_0^{*}=k_1^{*}$ and $k_0^{*}\ne k_1^{*}$, where the coefficient $a$ is related to the phase difference ${\mathsf{\Delta }\mathsf{\Phi }=\mathsf{\Phi }}_1\left[k_1^{*}\right]-{\mathsf{\Phi }}_0\left[k_0^{*}\right]$ and the coefficient $b$ to the amplitude ratio $(1-A_0/A_1)$ under the first-order approximation of Equation (16). When $A_0=A_1$, and thus $b=0$, the frequency offsets ${\epsilon }_0$ and ${\epsilon }_1$ are determined only by the phase difference, and ${\epsilon }_0={\epsilon }_1$. These are the low-noise cases discussed by Huang and Xia [9]. In contrast, when $A_0\ne A_1$, and thus $b\ne 0$ in high-noise cases, the frequency offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ are also dependent on the difference between $A_0$ and $A_1$. The larger the difference is, the more influence on ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$. If

$$a^2\le \left|b\right|$$

(21)

The frequency offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ are mainly determined by the difference between the $A_0$ and $A_1$ of the double-segment DFT. So, Equation (21) provides a basis for judging the difference between $A_0$ and $A_1$.

3.5. Application of the DS-DFT Algorithm to the DASH Interferograms

As application examples of the DS-DFT algorithm, the DASH interferograms without and with the measurement noises were taken as the input signal sequences. The noisy interferograms are shown in Figure 1; the noiseless ones are not shown. As mentioned in Section 2, the Fizeau interference fringes corresponding to the same atmospheric height appeared simultaneously on the same CCD row (in the horizontal direction of Figure 1) and were taken as an input signal sequence to perform the DS-DFT calculations in a manner of row by row. The altitude profiles were obtained for the phases ${\mathsf{\Phi }}_1$ and ${\mathsf{\Phi }}_2$, the amplitudes $A_1$ and $A_2$, and the frequency offsets ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$, as well as other quantities. Here and hereafter, the subscripts (1, 2) denote the two subsegments, respectively, replacing the (0,1) previously used in the fashion of Huang and Xia [9]. The calculations were performed for the DASH interferograms of the 557 nm and 630 nm airglow emissions observed by two FOVs. Figure 6 and Figure 7 show the calculation results of the DS-DFT algorithm for the DASH interferograms of the 557 nm emission observed by FOV1 in the absence and presence of the measurement noises. The results for the 557 nm emission in FOV2 and the 630 nm emissions in both FOVs showed similar features discussed below and are not presented.

The first panels from the right of Figure 6 and Figure 7 show the phases ${\mathsf{\Phi }}_1$ (black), ${\mathsf{\Phi }}_2$ (red), and their difference ${\mathsf{\Delta }\mathsf{\Phi }=\mathsf{\Phi }}_2-{\mathsf{\Phi }}_1$ (blue), computed from Equations (8) and (9). The phases ${\mathsf{\Phi }}_1\left(k_0^{*}\right)$ and ${\mathsf{\Phi }}_2\left(k_1^{*}\right)$ obtained directly from the DS-DFT spectral amplitude maxima (called DFT phases) are denoted by solid lines. Also, the phases ${\mathsf{\Phi }}_1\left(k\right)$ and ${\mathsf{\Phi }}_2\left(k\right)$ are fitted by a linear function $\mathsf{\Phi }\left(k\right)=a+b·k$, respectively. The phases ${\mathsf{\Phi }}_1\left(k_0^{*}\right)$ and ${\mathsf{\Phi }}_2\left(k_1^{*}\right)$ are calculated from the linear functions at $k=k_0^{*}$ and $k_1^{*}$ (called fitting phases), respectively, and indicated by dotted lines (called fitting phases). The DFT phases are seen to be quite close to the fitting phases, indicating a linear trend of the phase spectrum $\mathsf{\Phi }\left(k\right)$ with frequency $k$. This lays the foundation for the estimation of the Doppler phase shift either from the DFT phases in the Huang–Xia algorithm or from the fitting phases in the Englert–Harding algorithm.

As shown in Equations (8) and (9), phases ${\mathsf{\Phi }}_1$,${\mathsf{\Phi }}_2$, and their difference $\mathsf{\Delta }\mathsf{\Phi }$ are associated with the signal frequency $f_0$ and frequency offsets${ \epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$. In the case discussed here, the signal frequency $f_0$ and frequency offsets contain the Doppler shift due to the relative motion of the airglow emitting particles with respect to the detector, including the motion induced by the atmospheric wind field. Therefore, the height variations of ${\mathsf{\Phi }}_1$, ${\mathsf{\Phi }}_2$, and $\mathsf{\Delta }\mathsf{\Phi }$ mainly reflect the altitudinal change in the relative motion between the emitting particles and the detectors. The conversion factor of the Doppler phase shift to the relative velocity of motion is $c/\left(2\pi k_0\right)\approx 1.21\times {10}^5{\mathrm{m}\mathrm{s}}^{-1}{\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}}^{-1}$, where the speed of light is $c\approx 3.0\times {10}^{10} \mathrm{c}\mathrm{m}·{\mathrm{s}}^{-1}$, and Fizeau frequency $k_0=13.773 {\mathrm{c}\mathrm{m}}^{-1}$. The Doppler phase shift is very small. Even assuming that the relative velocity of the emitting particle and detector is 8000 m/s (including the velocity of the satellite), the Doppler phase shift is only about $1.15 \mathrm{m}\mathrm{r}\mathrm{a}\mathrm{d}={0.066}^{°}$. In the noiseless case (Figure 6), ${\mathsf{\Phi }}_1$,${\mathsf{\Phi }}_2$, and $\mathsf{\Delta }\mathsf{\Phi }$ show neither a significant height change nor this small Doppler phase shift. In contrast, in the noisy case (Figure 7), the altitude behavior of ${\mathsf{\Phi }}_1$,${\mathsf{\Phi }}_2$, and $\mathsf{\Delta }\mathsf{\Phi }$ is quite different, showing significant height variations. This is thought to be caused by taking the continuous arguments between the phase angles of two adjacent heights. When the absolute difference of the two phases exceeds 3.5 rad, the phases are reprocessed to ensure that the absolute difference never exceeds 3.5 rad. As described in Section 2, the phase errors caused by the measurement noises are estimated at $\pm 0.02 \mathrm{r}\mathrm{a}\mathrm{d}=\pm 20 \mathrm{m}\mathrm{r}\mathrm{a}\mathrm{d}=\pm 1.15°$. This error may result in the different values of the continuous arguments. Meanwhile, ${\mathsf{\Phi }}_1$,${\mathsf{\Phi }}_2$, and $\mathsf{\Delta }\mathsf{\Phi }$ in Figure 7 also show obvious height changes. There are significant differences between the regions below and above 105 km, with larger height changes below 105 km and smaller ones above. This corresponds to the rapid reduction in airglow emission rate above 105 km and the obvious change in the atmospheric wind field near 105 km in the considered case (not shown).

The second panels from the right of Figure 6 and Figure 7 show the DS-DFT amplitudes $A_1$ (black), $A_2$ (red), and their ratios $A_1$/$A_2$ (blue). They are associated with the frequency shifts ${\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ ($-0.5\le \left[{\epsilon }_0 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1\right]\le 0.5)$, as shown in Equation (14). In the absence of the measurement noises (Figure 6), $A_1$ and $A_2$ at all heights are equal, and their height variations mainly reflect the changes in the airglow emissivity at different heights. The ratio $A_1$/$A_2$ is identical to 1 and has no height variation. In the noisy case (Figure 7), the signal-to-noise ratios (SNRs) are different heights. In the 85–105 km region of high signal intensity and high SNR level, the DS-DFT amplitudes $A_1$ and $A_2$ are nearly equal, and their ratios $A_1$/$A_2$ are close to 1. However, in the region above 105 km where the signal intensity is weak and the SNR level is low, the ratios $A_1$/$A_2$ are significantly different from 1, with deviations of about 3–12%, and show large height changes.

The second panels from the left of Figure 6 and Figure 7 show the frequency offsets ${\epsilon }_1$ (black) and ${ϵ}_2$ (red) and the coefficients $a/2\pi$ (blue) and $\left(b/a\right)/2\pi$ (green), while the first panels of Figure 6 and Figure 7 display the Fizeau fringe frequencies $k_0$ (blue), $k_1$ (black), and $k_2$ (red) roughly estimated from the entire sampling sequence and from the DS-DFT spectra, as well as the fine estimation of the fringe frequencies $\widehat{f}=k_{es}$ (green) calculated from Equation (6). In the noiseless case (Figure 6), $k_1$ and $k_2 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ and ${\epsilon }_2$ from Equation (18) do not show significant height changes. In contrast, in the noisy case (Figure 7), $k_1$ and $k_2 \mathrm{a}\mathrm{n}\mathrm{d} {\epsilon }_1$ and ${\epsilon }_2$ from Equation (20) have no appreciable height change below 105 km but exhibit quite obvious height variations above 105 km. This corresponds to the relatively high SNR below 105 km and the low SNR above 105 km due to a sharp reduction in airglow emissivity. The fine fringe frequency $\widehat{f}=k_{es}$ is the key parameter of this study. At individual heights of the high SNR regions (the noiseless case in Figure 6 or below 105 km of the noisy case in Figure 7), when the DS-DFT estimated Fizeau frequency $k_1=k_2$ (black and red dots coincide), the fine frequency $\widehat{f}=k_{es}$ (green) estimated from Equation (6) is between $k_1=k_2$ and $k_0$(blue) estimated from the entire sampling sequence, implying that $\widehat{f}=k_{es}$ provides an accurate estimation in comparison with $k_0$ as shown by Huang and Xia [9]. However, in the low SNR heights above 105 km in Figure 7, the DS-DFT estimated Fizeau frequency $k_1\ne k_2$, and the fine Fizeau frequency $\widehat{f}=k_{es}$(green) from Equation (6) tends to the Fizeau frequency ${\kappa }_0$(blue) estimated from the entire sampling sequence, suggesting that $\widehat{f}=k_{es}$ is unlikely better than $k_0$ in the low SNR case. Nevertheless, the DS-DFT algorithm still can provide a criterion of noise level determined from the observed interferograms, as discussed below.

4. DS-DFT Criterion of Noise Level for the DASH Interferogram

4.1. DS-DFT Criterion of Noise Level

According to the discussion in Section 3.5, the frequencies $k_1^{*}$ and $k_2^{*}$ and the amplitudes $A_1$ and $A_2$ determined at the DS-DFT spectral amplitude maxima can be used to qualitatively characterize the noise level of the DASH interferograms. Three empirical scenarios can be obtained:

The two frequencies and two amplitudes are equal, respectively, namely, $k_1^{*}=k_2^{*}$ and $A_1=A_2$. This represents the case of the low-noise level, including the noise-free case.

The two frequencies are equal, but the two amplitudes are not, namely, $k_1^{*}=k_2^{*}$ and $A_1\ne A_2$. This represents the case of the moderate-noise level.

The two frequencies are not equal, i.e., $k_1^{*}\ne k_2^{*}$, usually implying that the two amplitudes are also not the same, $A_1\ne A_2$. This represents the case of the high-noise level.

For convenience, the cases of the moderate- and high-noise levels are combined together and referred to as the case of the high-noise level.

4.2. SNR Calculation

The effectiveness of the DS-DFT criterion can be verified by comparison with the SNR calculated directly from the simulated interferograms. The SNR is usually defined as the ratio between the effective powers of signal and noise in the form of [14]

$$SNR=10\cdot {\log }_{10}{\displaystyle \frac{P_S}{P_n}}(\mathrm{db})$$

(22)

However, the signal and noise power $P_S$ and $P_n$ is usually difficult to calculate and can be approximately converted into the ratio between the effective amplitudes of signal and noise, $I_S$ and $I_n$:

$$SNR=20\cdot {\log }_{10}{\displaystyle \frac{I_S}{I_n}}(\mathrm{db})$$

(23)

where ${P_S\propto I}_S^2$ and ${P_n\propto I}_n^2$ are assumed. In order to obtain the effective amplitude $I_S$ and $I_n$, the simulated interferogram without the measurement noise (not shown) is firstly calculated, and the interference intensity $I_S^0$ of each pixel (or pixel bin) is taken as the effective amplitude of the signal. Let the interference intensity of each pixel (or pixel bin) in the noisy interferogram (Figure 1) be $I_S$, the effective amplitude of the noise can be estimated from $I_n=\left|I_S-I_S^0\right|$. Thus, the SNR can be expressed as

$$SNR=20\cdot {\log }_{10}{\displaystyle \frac{I_S^0}{\left|I_S-I_S^0\right|}}(\mathrm{db})$$

(24)

4.3. Comparison of the DS-DFT Criterion with the SNR Calculation

Figure 8 shows the SNRs of each pixel bin calculated according to Equation (24) (the left panels) and the CCD-column-averaged SNR of each row (i.e., averaged SNR of each height) (the right panels) for the FOV1 simulated interferogram of 557 nm airglow emission (Figure 1). An SNR greater than, equal to, and less than zero means that the amplitude of signal is greater than, equal to, and less than that of the noise, correspondingly. In both figures, it can be seen that the SNRs below 105–107 km are greater than zero, indicating a low-noise level, while the SNRs above 105–107 km are close to or less than zero, indicating a high-noise level.

The noise level at each height is also determined by the DS-DFT criterion. The high- ($k_1^{*}\ne k_2^{*}$), moderate- ($k_1^{*}=k_2^{*}$, but $A_1\ne A_2$), and low- ($k_1^{*}=k_2^{*}$ and $A_1=A_2$) noise levels are indicated by the pink diamond, pink square, and green circle, respectively. It can be seen that the noise levels determined from the DS-DFT criterion are quite consistent with those calculated from the SNR. In practice, it is very difficult, even impossible, to separate the information of signal and noise (i.e., the effective power $P_S$ and $P_n$ or the effective amplitude $I_S$ and $I_n$) from realistic measurements and to calculate the SNR from Equations (22)–(24). Therefore, the DS-DFT criterion provides an operational algorithm to determine noise level directly from observational data.

4.4. Application of the DS-DFT Criterion at Different Noise Levels

In order to investigate the performance of the DS-DFT criterion at different noise levels, a series of simulated DASH interferograms was generated with random noises of magnitudes at 0 to 80% and an increment of 10% of the simulated emission intensities. Figure 9 shows the cases of the 10%, 20%, and 30% noise levels. For all three cases, the SNRs below 105–107 km are greater than zero, indicating a low-noise level, while the SNRs above 105–107 km are close to or less than zero, indicating a high-noise level. In contrast, the DS-DFT criterion indicates the existence of high-noise points below 105–107 km for the 30% case. When noise levels increase up to 80%, more and more high-noise points appear below 105–107 km, which are still determined as low-noise levels by the SNR ratio, due to the large emission intensities at those altitudes. This implies that the DS-DFT criterion is more sensitive to the change in noise levels.

5. Discussion and Conclusions

The interferograms observed by the spaceborne DASH instrument contain signals from atmospheric airglow emissions and measurement noises from readout noise, quantization noise, and dark current, as well as other sources (Figure 1). Usually, a discrete Fourier transform (DFT) is applied to the measured interferograms, and the derived amplitude and phase spectra are used to extract the information of the intensity and frequency of the airglow emission, including the very small Doppler shift frequency induced by the motion of the emitting particles with respect to the detector [3].

This study pays special attention to the phase spectrum extracted from the DASH interferograms of high noise. Various digital filtering schemes with a frequency response of zero-phase shift could be adopted to suppress the influence of noise but conserve the phase relations unchanged. The low-pass Gaussian filter and bidirectional Butterworth finite impulse response (FIR) filter were tested (not shown here) [15,16]. It was found that the low-pass filters can greatly suppress the large random oscillations in the wave amplitudes $⌈A\left({\kappa }_0,x\right)\beta \left({\kappa }_0\right)⌉$ and phases $\mathsf{\delta }{\mathsf{\Phi }}_{tot}\left(x,{\kappa }_0\right)$ of the unilateral interferogram and enhance the linear feature of the phase spectrum. However, these types of low-pass filters have an inherent deficiency near the zero frequency, where the low-frequency components associated with noise are still perceptible, and their phases remain nearly unchanged after filtering. The derived Doppler phase shift $D_{\kappa }\left({\kappa }_0\right)$ may still have a large error.

In order to compensate for the deficiency of the Butterworth FIR Filter, the warped filter (WFIR) can be used [17]. However, WFIR has a nonlinear frequency scale, which requires a high resolution at low frequency, complex structure, and time-consuming operation. Thus, instead of applying the WFIR filter, this paper sought other simple and fast methods for processing the interferograms of high noise. The two-subsegment Fourier transform algorithm proposed by Huang and Xia [9] for low-noise signals is extended to the high-noise case, and fine estimation of the frequency spectrum is obtained (§3). Then, a criterion for determining the noise level of the measured signal (the observed interferogram) is proposed and compared with the usual calculation results of the SNR (§4) to provide an operational method for processing the observed signals.

The noise levels determined from the DS-DFT criterion are quite consistent with those directly calculated from Equations (22)–(24) for the SNR of the simulated interferograms, implying the validity of the DS-DFT criterion. Since it is very difficult, even impossible, to extract the information of signal and noise and directly calculate the SNR from realistic measurements, the DS-DFT criterion provides an effective way to determine noise level directly from observational data.

The consistency between the proposed DS-DFT criterion and the SNR calculation suggests the validity and feasibility of the criterion. Noise contamination influences the anti-noise robustness and estimation accuracy of a frequency estimator. The signal-to-noise ratio (SNR) reflects the anti-noise robustness. Below an SNR threshold, the estimator’s mean squared error (MSE) will exhibit a rapid increase, potentially resulting in an inaccurate estimation [14]. The MSE is also used to evaluate the frequency estimation accuracy arising from noise contamination. Physically, the key indicator MSE indicates the spectral leakage surrounding a spectral peak. For a periodic signal, the minimum MSE or minimum spectral leakage corresponds to the case of the signal energy concentrating on the single spectral peak, resulting in the estimator attaining optimal anti-noise performance [5].

However, it must be noted that our experimental results have only shown the basic function of the DS-DFT algorithm; more work is required to fully prove the effectiveness of the innovation and to verify its robustness and accuracy under different environmental conditions. For example, an analysis of statistical error and confidence intervals should be supplemented to enhance the persuasiveness of the data results. The use of a statistical test (such as the k-fold approach) should be conducted in order to assess the stability of the results shown in the paper. A quantitative metric to emphasize the improvement of our method in fine frequency estimation and noise-level classification should be set first, and then a more detailed comparison analysis of DS-DFT performance under different geophysical conditions, noise levels, and characteristics (such as colored noise and Gaussian noise), as well as a comparison with other high-noise spectrum estimation methods, is required. The future development of these challenge directions will further verify the applicability and robustness of the method and clearly show its improvement over existing methods.