Commonly used methods fail to detect known propagation speeds of simulated signals from time-longitude (Hovmöller) diagrams

This work examines the accuracy and validity of two variants of Radon transform and two variants of the TwoDimensional Fast Fourier Transform (2D FFT) that have been previously used for estimating the propagation speed of oceanic signals such as Sea Surface Height Anomalies (SSHA) derived from satellite borne altimeters based on time-longitude (Hovmöller) diagrams. The examination employs numerically simulated signals made up of 20 or 50 modes where one, 10 randomly selected, mode has a larger amplitude than the uniform amplitude of the other modes. Since the dominant input mode is ab-initio known, we can clearly define "success" in detecting its phase/propagation speed. We show that all previously employed variants fail to detect the phase speed of the dominant input mode when its amplitude is smaller than 5 times the amplitude of the other modes and that they successfully detect the phase speed of the dominant input mode only when its amplitude is at least 10 times the amplitude of the other modes. This requirement is an unrealistic limitation on oceanic 15 observations such as SSHA. In addition, three of the variant methods “detect” a dominant mode even when all modes have the exact same amplitude. The accuracy with which the four methods identify a dominant input mode decreases with the increase in the number of modes in the signal. Our findings are relevant to the reliability of phase speed estimates of SSHA observations and the reported “too fast” phase speed of baroclinic Rossby waves in the ocean.

De-Leon and Paldor (2017a) examined the accuracy of various methods in estimating the phase speed of waves by applying these methods to an artificially generated signal made of 3 sine functions (modes) with known phase speeds and amplitudes, compounded by large-amplitude random white noise. All methods have successfully filtered out the high amplitude white-noise from the 3-harmonic signal and accurately detected the main mode (some of them also detected the secondary modes). However, such a signal is too synthetic/ideal and cannot be compared to real oceanic observations that 5 include tens, if not hundreds, of modes with different frequencies and propagation speeds and not just 3 modes compounded by white noise.
In this short study we simulate "oceanic observations" and examine whether the methods detect a single dominant propagation speed out of many (20 or 50) speeds. In Sect. 2 we provide details on the generation of "observed"  signal, the methods for evaluating the propagation speed of the signal and the tests we apply to the signals. The results are shown in 10 Sect. 3 and discussed in Sect. 4.

Generating the simulated "observations"
The  signals (where  is any oceanic variable) used here are generated numerically by summing up N purely propagating sine functions (modes, hereafter) of the form sin(kx-t) where k is the zonal wavenumber, x is longitude, t is time and =kC 15 is the frequency (where C is the zonal propagation speed). The number of participating modes, N, is taken to be either 20 or 50 and N-1 of these modes have an amplitude of 1 while the amplitude of the additional N th , randomly chosen, mode is either larger than or equal to 1. The sum of all N modes constitutes the -signal which is analyzed by the methods described in 2.2.
The "spatial domain" (x) is chosen between "longitudes" 70-130° on a Cartesian grid with a 1/4° resolution, and the "time" (t) duration is 20 years (1044 weeks) with temporal resolution of one datum per week, similar to publicly distributed 20 products by e.g. Aviso. The values of the propagation speeds, C, are uniformly distributed in the -18 to 0 cm s -1 range (i.e. each of the 20 or 50 modes is assigned a different propagation speed within that range), which is typical for baroclinic Rossby waves in the ocean (see e.g. Fig. 7 of Killworth et al., 1997;Barron et al., 2009). The values of the frequencies, , are selected randomly so that the period, 2/, falls in the range between 5 and 200 weeks while the values of the zonal wavenumbers, k, equal /C. The resulting signal was low-pass filtered by applying a 5-week-running-average at each grid 25 point to eliminate short term variations.
The  signal made up of the filtered signal (i.e. the sum of N pure sine waves) at a given latitude, is plotted as a function of longitude and time (Hovmöller diagram). When a single dominant mode exists that has a certain propagation speed the pattern on this diagram is a straight line whose slope is the inverse of the dominant propagation speed (since the abscissa is longitude and the ordinate is time). An example of a time-longitude diagram of an artificial signal is shown in Fig. 1a (for speed of the dominant input mode. The challenge is to estimate the dominant speeds using different methods and examine their success in detecting the propagation speed of the known dominant input mode. The methods examined here are detailed in the next subsection.

Methods for estimating the "observed" propagation speed
Four variants of methods have been employed for identifying the preferred direction of the same-amplitude contours on the 5 Hovmöller diagram, each variant relates a certain measure of the power/intensity in a mode to its propagation speed. The first method is the Radon transform, used by e.g. Chelton and Schlax (1996); Chelton et al. (2003) and Tulloch et al. (2009) for analysing satellite observations of the ocean. In this method, one calculates the sum of the amplitudes along lines inclined at an angle  and displaced a distance s from the origin. Then, the sum of squares of the values of these sums along all lines having the same angle is calculated and the angle at which this sum of squares is maximal, is the best estimate for the 10 orientation of the lines on the image. The dominant propagation speed of the signal is then proportional to the tangent of this angle of maximum sum-of-squares. In the second method, which is a variant of the Radon transform, the variance of the amplitudes is calculated along every angle  instead of the sum of amplitudes. This method was applied e.g. by Polito and Liu (2003 to local auto-correlation of segments of Hovmöller diagram) and Barron et al. (2009). Another, independent, method commonly used (e.g. Zang and Wunsch, 1999;Osychny and Cornillon, 2004) is the two-dimensional Fast Fourier 15 Transform (2D FFT). An example of (ω, k) diagram, obtained by applying 2D FFT to the signal of Fig. 1a is shown in Fig. 1b. Here we use two variants of the 2D FFT method: In the first variant one sweeps over the 2D FFT spectra to find the direction in (ω, k) plane with maximum "energy" (this is the third method) while in the second variant one finds the maximal amplitude of the 2D FFT (i.e. one of the bright points in Fig. 1b) and calculates the ratio ω/k where ω and k are the frequency and zonal wavenumber of the maximal amplitude (this is the fourth method). Detailed description of these methods and the 20 interpretation of observed signals are found in De-Leon and Paldor (2017a).
Each of the four variants of the methods can yield an estimate of the dominant propagation speed of a signal, based on the extremum of a graph that relates the calculated measure of a mode's intensity to its propagation speed. An estimation of the propagation speed based on a local extremum of this graph is accepted when this extremum is narrow and isolated compared to other local extrema. When the normalized amplitude (the term normalized amplitudes is used for the calculated 25 amplitudes divided by the maximal amplitude in the domain) of a distinct peak is 1 while the (normalized) amplitudes of all other peaks are smaller than 0.8, this mode is considered the dominant mode. In the variance method, where the extrema are minima, the dominant mode is accepted when its amplitude is 0.0 while the normalized amplitudes of all other minima are larger than 0.2. The results shown in Fig. 2 demonstrate the emergence of a "rejected" peak that does not differ significantly from other peaks (2d) and "accepted" peaks whose intensities are significantly larger than those of other peaks (2a, 2c) or 30 "accepted" trough that is significantly smaller than the other troughs (2b).

Examining the accuracy of dominant mode detection
Two types of tests are applied to these "observations" to examine the accuracy of the various methods in assessing the existence of a dominant propagation speed (i.e. mode) in the signal.
The first is a true-positive/false-negative test in which there is a dominant mode in the "observed" signal and a given method indicates whether this dominant mode exists (true-positive; TP) or not (false-negative; FN). In our case, one of the sine functions (this is the dominant input mode) is chosen randomly and its amplitude is set to be larger than 1. We check for 5 each method if (at all) it identifies a dominant mode and if so, if it matches the propagation speed of this (larger amplitude) input mode. We divide the interval of propagation speeds into N (=20 or 50) bins of equi-distant values and the determination of the success of the methods in identifying a dominant mode is as follows: if the dominant mode found by the method falls in the expected bin of the dominant input modewe score it by 1 ("TP"). If it is found in one of its next neighbors, it is scored by 1/2. A score of 0 ("FN") is assigned when the method cannot find any dominant mode and when it 10 detects a dominant mode more than 1 bin away from the correct bin. The second is a false-positive/true-negative test in which no dominant mode exists in the "observed" signal and a given 15 method indicates that there exists a dominant mode (false-positive; FP) or not (true-negative; TN). In our case, this is done by generating a signal in which all modes have identical amplitudes (=1) and checking whether a method erroneously detects a certain propagation speed as dominant. If a dominant mode is detected, we score it by 1 ("FP"), if a peak is detected but is too wide we score it by 1/2 and if there is no dominant mode (i.e. no distinct peak or more than one peak) we score it by 0 ("TN"). We repeat this procedure 50 times for each of N=20 or N=50, sum up the scores and calculate the percentage of 20 erroneous detection of dominant mode in the signals by FP/(FP+TN)*100.

Results
An example of false determination of the dominant mode is shown in Fig. 2 for the signal shown in Fig. 1a. Figure Figure 2b shows the (normalized) mean of variances as a function of C (black markers), and here, too, the calculated propagation speed (the curve's minimum point) does not agree with the dominant input mode's 30 speed. Figure 2c shows the (normalized) distribution of the sum-of-squares of the spectral coefficients (2D FFT-amplitudes) along different ω/k lines (sweeping) as a function of C (black markers) where a distinct peak exists but it's located far from the dominant input mode's propagation speed. Figure 2d shows the 20 highest (normalized) 2D FFT amplitudes of the (ω, k) diagram of Fig. 1b. There are many peaks with no clear single maximum and the dominant input mode has one of the lowest amplitudes. For this signal, none of the methods identified correctly the dominant input mode.
The statistics of success of each method in detecting the dominant input mode for the TP/FN test is shown in Fig. 3 for dominant input mode's amplitude of 2.5, 5 and 10. The results for amplitudes of 1.5 and 2 are not shown as the success rate 5 of all methods at these amplitudes is between 10% and 30% only. In each case we repeated the procedure 50 times (i.e. for 50 signals) for N=20 and then for N=50, summed the scores and calculated the percentage of success by TP/(TP+FN)*100.
The conclusions from these results are: 1) in order to identify the dominant input mode with more than 70% certainty, its amplitude should be larger than 5. 2) No method has clear advantage over the other methods. 3) Clearly, as N increases the dominant mode's amplitude has to increase, too, for a successful identification (so as to ensure that the ratio between the 10 dominant mode's amplitude and the sum of all amplitudes is similar for different values of N since, e.g., 2.5/20>2.5/50). 4) The amplitude at which successful detection occurs decreases with the decrease in the ranges of propagation speeds and periods. Thus, for propagation speeds in the range of -10 to -2 cm s -1 and periods between 15 and 100 weeks, an amplitude of 5 is successfully detected in over 90% of the cases, while an amplitude of 2 yields poor results (results not shown).
The statistics of erroneous detection of a dominant mode (the FP/TN test where the amplitudes of all input modes equal 15 1 i.e. there is no dominant input mode) is shown in Fig. 4 for each of the methods (here again we have generated 50 signals for N=20 and for N=50 and calculated the percentage of erroneous detection by FP/(FP+TN)*100). The 2D FFT maxima method is the only method for which the percentage of error is smaller than 20% while other 3 methods err in at least 50% of the cases (i.e. they identify a single clear peak in one of the propagation speeds). As N increases this erroneous detection percentage decreases slightly in the 2 Radon variants but increases slightly in the 2D FFT sweeping method. 20

Discussion
None of the methods can identify a dominant input mode unless its amplitude is significantly larger than the others (by a factor larger than 5! in the present study's ranges of propagation speeds and periods) and most of them (except the 2D FFT maxima) erroneously detect a dominant mode when there is no such input mode. Though the 2D FFT maxima method does not falsely detect dominant mode when it does not exist, its performance in detecting dominant input mode when it exists is 25 not satisfactory. For realistic signals of the ocean we don't know that there is a dominant mode with sufficiently large amplitude so none of the methods is reliable for estimating the propagation speed of e.g. Rossby waves.
When the values of  and k are chosen in a range corresponding to the resolution limit and the Nyquist frequency, the success of 2D FFT in identifying the dominant input mode's propagation speed increases significantly, compared to the case where the values of C are determined in advance,  is chosen randomly in a particular range and k is set as a result of /C so 30 aliasing can occur. For that reason, even if the signal includes only one mode (i.e. one C), and both  and k are chosen in the latter manner, there can be a wrong identification by the 2D FFT (but the Radon transform identifies it correctly). In the ocean we don't know ab-initio which wave numbers and frequencies exist so we cannot filter them out of the signal; hence aliasing can occur, and the percentage of success in detecting the real dominant mode is expected to decrease further.
The erroneous identification of the Radon and variance methods can be partially attributed to the non-linear relation between the angle  and the propagation speed, C, which is proportional to tan( ) so the equi-distant values of C are converted to  -values that are very close to one another. Figure 5 shows the distribution of the sum of squares of the Radon 5 transform versus  for the signal shown in Fig. 1a (while the distribution of the sum of squares of the Radon transform versus C for that signal is shown in Fig. 2a). It is clear from this figure that the peak is located in the vicinity of  values corresponding to many C values. However, the performance of the Radon variants improves with the increase of the dominant input mode's amplitude, so the non-linear relation between the angle and propagation speed is not the only reason for the mismatch. 10 As N (the number of modes) increases (it is impossible to establish a-priori a bound on the number of modes in the ocean), the dominant input mode's amplitude should be larger in order to be separately identified from modes with similar characteristics. In addition, the width of the bins becomes narrower as N increases so fewer results can be evaluated as success. Also, in narrower ranges of propagation speeds and periods and for the same number of bins, all methods correctly detect the dominant input mode at lower threshold amplitudes (results not shown). 15 The weakness of the methods in identifying the dominant mode points to the difficulty in comparison between theories and observations of baroclinic Rossby waves in the ocean and this difficulty might explain the lack of "continuity" of propagation speed estimates between adjacent latitudes in one (or more) methods. It can also explain why a validation of the higher order trapped wave theory (where the  term is treated consistently) has been confirmed by observations only in the Indian Ocean south of Australia (De-Leon and Paldor, 2017b) and not in other parts of the world ocean. In contrast to the 20 harmonic theory whose propagation speed estimates are always slower than the observed speeds, the trapped wave theory differs from observations sporadically (see Fig 5 of De-Leon and Paldor, 2017b).