Local Frequencies in Superoscillatory Phenomena

Abstract

Superoscillations correspond to a non-linear phenomenon theoretically addressed by Aharonov in 1991. The resulting waves or functions have the particularity of being of limited bandwidth and contain faster amplitude variations than that corresponding to the fastest components obtained applying the Fourier transform. Also, the amplitude developed in the region where it occurs is small, since it decreases exponentially. These characteristics prevent its determination using the Fourier transform since it is not a stationary phenomenon. With this perspective, we have tested other methods for determining these features, such as wavelet transforms and Hilbert-Huang transform. Wavelet transforms can capture both low- and high-frequency components of the signal. The Hilbert-Huang transform allows the decomposing of a signal into the so-called intrinsic mode functions (IMF) together with a trend, and obtaining instantaneous frequencies. We also proposed a methodology using Gabor-adaptive windows to perform detection. Finally, filtering results were added using a multiresolution analysis decomposition that allows separating the super-oscillatory part of one and therefore localizes the oscillations in time, that is, local frequencies.

1. Introduction

Superoscillation (SO) is a counterintuitive non-linear phenomena introduced first theoretically by Aharonov in 1991 [1]. In order to simplify the study, we work with some particular bandlimited functions introduced by Berry [2] having variations faster than the corresponding harmonic obtained when applying the Fourier transform. The main properties of these functions are briefly summarized.

The fast variation of the amplitude occurred in a small temporal domain in some cases less than one period. In consequence, to analyze these localized features it is possible to use alternative techniques to Fourier such as wavelets and the empirical mode decomposition (EMD) that allows treating signals expanding the traditional Fourier signal model. In signal processing literature, this “expansion” [3] is described as an adaptive harmonic model. EMD and wavelets are introduced and described likewise results were obtained for SO signals. One of the first calculations used to determine the instantaneous frequency was introduced by Berry but it failed for some SO functions [4]. Although there are current transformations that outperform them, timely modifications have been proposed such as adaptive and dynamical windows as a part of the Gabor transform. This window is designed to adjust its width in agreement with the maxima and minima of the analyzed signal and in consequence, detects features of the signal that correspond with SO.

Additionally, we introduce other discussions such as continuous and discrete wavelet filtering methods to localize SO through multiresolution analysis.

2. First steps with superoscillations

SO was first introduced theoretically by Aharonov in 1991 in the context of weak measurement performed on an ensemble of a quantum system. As a result, he obtained eigenvalues larger than the largest eigenvalue of the observable. In a similar form, there exist some particular functions having variations faster than the corresponding harmonic obtained when the Fourier transform is applied. Additionally, these functions are bandlimited, making this phenomenon more surprising. On the other hand, the fast variation of the amplitude occurred in a small temporal domain in such a way that the wave does not describe a full period in the majority of the cases.

2.1 Mathematical background and Fourier transform

To tackle the study, we begin first considering the expressions given by Berry. He generalized the mathematical properties of SO introducing methods to construct band limited signals containing local oscillations. SO are band-limited functions with the counterintuitive property that they can vary arbitrarily faster than their fastest Fourier component, over arbitrarily long intervals. Modern studies had been made in quantum theory but were anticipated in radar and optics. The mathematical understanding – still being explored – recognizes that functions are extremely small where they superoscillate; this has implications for information theory. Optical vortices, sub-wavelength microscopy, two level systems and related areas of nanoscience are now moving from the theoretical and the demonstrative to the practical. Our work uses a typical SO function formulated by Berry.

2.2 Theoretical considerations

Probably, one of the most studied and simplest is the periodic SO function given by Eq. (1) in which n is a large even integer and a > 1. Its period is nπ, and is band-limited, because when expanded in a Fourier series, the component oscillations are all of the form exp(iknt) with |kn| < 1. Such mathematical expressions could be formulated in terms of cosine and sine functions.

$$\mathit{Re}\left\{F\left(n,a,t\right)\right\}=\mathit{Re}{\left[\mathit{cos}\left(\frac{{\omega }_{0}t}{n}\right)+\mathit{ia}\mathit{sin}\left(\frac{{\omega }_{0}t}{n}\right)\right]}^n$$E1

To clarify the meaning of each term used in Eq. (1), ω₀ is the highest frequency in the Fourier spectrum of the function, ${}^{“}{a}^{”}$ is usually called the degree of SO in the region near t = 0. Finally, n measures the extent of this SO region. It is SO, because for |ω₀n| < √n it can be approximated by exp(iax) as is demonstrated in [1]. Outside the interval |ω₀tn| < n, f(x) the function increases anti-Gaussianly up to its maximum value |f(nπ/2)| = an. We introduce here an example when n = 5. In Figure 1 three harmonics have been represented and ω₀ corresponds to 200000π Hz (that is f₀ = 100 kHz). Additionally, the curve in red corresponds to the SO frequency and in black the resultant wave. It can be appreciated as an overlap between curves. It takes place in the region near t = 0.

Figure 1.

Representation of the super-oscillatory function and its harmonic composition (dots colored curves) for 20, 60, and 100 kHz. The black curve in solid was calculated from Eq. (2) for n = 5 and a = 2.

The harmonic Fourier components previously shown are easily appreciated if Eq. (1) is written in exponential form as:

$$\mathit{Re}\left(F\left[{\omega }_0,a,n\right]\right)=\sum _{j=0}^n{C}_j\left(n,a\right)\mathit{exp}\left(i{\omega }_0\left(1-\frac{2j}{n}\right)\right)$$E2

where

$$C_j\left(n,a\right)=\frac{n!{\left(\frac{1-a}{2}\right)}^j{\left(\frac{a+1}{2}\right)}^{n-j}}{j!\left(n-j\right)!}$$E3

are the coefficients in the summation, where it is clear that if n is odd, there are (n + 1)/2 term representing the harmonic components of the wave. However, the bandwidth of this wave is higher than the highest Fourier component. The Fourier transform of Eq. (2) is

$$\mathit{\text{FT}}\left(\omega ,a,n\right)=\sum _{j=0}^n{C}_j\left(n,a\right)\delta \left(\omega -{\omega }_0\left(1-2j/n\right)\right)$$E4

Eq. (4) can be illustrated following the previous example for n = 5 and $a$ = 2. Here is shown the real part of the SO function in the frequency domain. Evidently the function is, after applying the Fourier analysis, limited to 100 kHz.

This last feature makes Fourier transform inadequate; partial overlapping observed in Figure 2 between the resultant function and the SO wave, suggests the fact that SO is a non-stationary phenomenon.

Figure 2.

Fourier transform for n = 5 and a = 2 of the function given in Eq. (1). The height of each stem corresponds with the magnitude of C_j(n,a) coefficient given in Eq. (4). This function has a discrete spectrum and is bandlimited at 100 kHz. Changes in the sign of coefficients, correspond to phases 0 and π of the components.

At this point, it is necessary a brief review considering that signals with frequencies that do not change with time are called stationary signals able to be analyzed by Fourier series methods or the Fourier transform.

However real-life signals are commonly non-stationary; that is, their frequencies may change with time (i.e., spectrum varies in time). So, in this case spectral components at a certain instant are unknown. Despite this fact, it is possible to determine what spectral components exist at any given interval of time by means of a procedure known as “windowing” of the signal. What is more in agreement with the signal variation the window must change its temporal width to be efficient. This change implies a resolution adjustment of this window (scaling). The previous considerations result in a transform called wavelet, which possesses a variable resolution. The analysis of such signals requires expanding the traditional Fourier signal model. In signal processing literature, this “expansion” is described as an adaptive harmonic model [3]. So, non-stationary signals modify their frequencies in localized time intervals. Considering such signals, it is convenient to revise different definitions of instantaneous frequencies. Likewise, apart from Fourier, there exist other methods to analyze signals with frequencies not present concentrated in some temporal segments, for example, Hilbert Huang transform (H-HT) and wavelet transforms.

Table 1 summarizes the main similarities and differences between techniques for signals time-frequency analysis.

	Fourier	Wavelet	Hilbert-Huang
Basis	A priori	A priori	Adaptive
Frequency	Global uncertainty	Regional. Partial uncertainty	Local. Certainty
Presentation	Energy. Frequency	Energy-Time-Frequency	Energy-Time- Frequency
Non-linear	No	No	Yes
Non-stationary	No	Yes	Yes
Theoretical base	Complete	Complete	Empirical

Table 1.

Comparison between Fourier, Wavelet, and H-HT properties.

3. Definition of instantaneous frequency

An important feature to characterize any oscillatory function is the local frequency. Taking into account local changes of frequency for several types of signals, these definitions become more important. However, there exist various definitions to calculate this value because each of them is more or less efficient in calculating the instantaneous frequency in agreement with the considered signal.

In order to define an instantaneous frequency for x(t), the Hilbert transform Eq. (5) is used, through which we can find the complex conjugate, y(t), of any real-valued function x(t) of L^p class,

$$H\left[x\left(t\right)\right]=\frac{1}{\pi }\mathit{PV}\underset{-\infty }{\overset{\infty }{\int }}\frac{x\left(\tau \right)}{t-\tau }\mathit{d\tau }=\frac{1}{\pi }x\left(t\right)\ast \left(\frac{1}{t}\right)$$E5

in which PV indicates the principal value of the singular integral. An alternative interpretation of this integral is the convolution between x(t) and the function 1/t. Beginning from the Hilbert transform, and computing the real or the imaginary part we have:

$$z\left(t\right)=x\left(t\right)+\mathit{iy}\left(t\right)=a\left(t\right)e^{\mathit{i\theta }\left(t\right)}$$E6

Thus, z(t) Eq. (6) is an analytical function. Such functions allow determining the phase function in a relatively easy way.

Here

$$a\left(t\right)=\sqrt{x{\left(t\right)}^2+y{\left(t\right)}^2}$$E7

and

$$\theta \left(t\right)=\mathit{\text{atan}}\left(\frac{y\left(t\right)}{x\left(t\right)}\right)$$E8

The function $a\left(t\right)$ in Eq. (7) is the instantaneous amplitude and $\theta \left(t\right)$ given by Eq. (8) is the phase function, then it follows that:

$$\omega \left(t\right)=\frac{\mathit{d\theta }\left(t\right)}{\mathit{dt}}$$E9

As the instantaneous frequency is defined through a derivative Eq. (9), it is very local. It can be used to describe the detailed variation of frequency, including the intra-wave frequency variation. As simple as this principle is, the implementation is not at all trivial. To represent the function in terms of a meaningful amplitude and phase, however, requires that the function satisfies certain conditions [5]. In particular, for SO functions beginning from Eq. (6) we introduce expressions of instantaneous frequency. Based on Eq. (1) the SO signal can be expressed in the following way:

$$f_{\mathit{SO}}\left(t\right)=Re\left\{A\left(t\right)\mathit{exp}\left(i\underset{0}{\overset{t}{\int }}\omega \left(t^{\text{'}}\right){\mathit{dt}}^{\text{'}}\right)\right\}$$E10

Here we consider $f_{\mathit{SO}}=Re\left(F\left[n,a,t\right]\right)$.

Berry and Popescu [6] derived the local frequency function ω(t), considering that A(t) presents slow variations. With these considerations, the local frequency function can be written as:

$$\omega \left(t\right)=\mathit{Im}\left\{\frac{\partial }{\partial t}\mathit{log}\left(f_{\mathit{SO}}\left(t\right)\right)\right\}$$E11

for the case of a SO with parameters a and n, the expression is

$$\omega \left(t\right)=\frac{\mathit{Na}}{{\mathit{cos}}^2\left({\omega }_{0}t\right)+a^2{\mathit{sin}}^2\left({\omega }_{0}t\right)}$$E12

Eq. (12) describes the local frequency along the SO function showing clearly the points where the frequency rises. The SO feature is observed in Figure 3 as an interval around t = −1.5 μs as expected. In this case, from Figure 3 it can be seen that the maximum frequency is 10⁶/2π = 200 kHz in agreement with the product af₀.

Figure 3.

Representation of the instantaneous frequency obtained from Eq. (12) corresponding to a SO function for n = 5 and a = 2.

The above local frequency derivation can be generalized to any function comprising a sum of modes. If the sum of modes is defined with a set of amplitudes C_n and frequencies ω_n the local generalized frequency function ω(t) can be expressed as:

$$\omega \left(t\right)=\mathit{Im}\left\{\frac{\partial }{\mathit{\delta t}}\mathit{log}\left(\sum _{n=0}^N{C}_n\mathit{exp}\left(i{\omega }_{n}t\right)\right)\right\}$$E13

The H-HT includes both the empirical mode decomposition (EMD) and the Hilbert spectral analysis (HSA) methods. They were introduced 20 years ago for analyzing data from non-linear and non-stationary processes, in order to gain a deeper insight into the underlying process that generates the data that cannot be handled by Fourier-based analysis. Considering the SO features, they can be treated in a suitable form by this methodology.

Studies by Ref. [7] established that the EMD is a dyadic filter, and it is equivalent to an adaptive wavelet.

Components of the EMD are usually physically meaningful, for the characteristic scales and are defined by the physical data. The sifting process is, in fact, a Reynolds-type decomposition: separating variations from the mean, except that the mean is a local instantaneous mean, so that the different modes are almost orthogonal to each other, except for the non-linearity in the data. Once obtained from the intrinsic mode function components, we can apply the Hilbert transform on each IMF component and compute the instantaneous frequency as the derivative of the phase function. After performing the Hilbert transform on each IMF component, we can express the original data as the real part, in the following way

$$x\left(t\right)=\mathit{Re}\left\{\sum _{j=1}^n{C}_j\left(t\right)\mathit{exp}(i{\int }_0^t{\omega }_j\left(t^{\text{'}}\right){\mathit{dt}}^{\text{'}}\right\}$$E14

The above equation gives both, amplitude and frequency of each component as a function of time. The same data, if expanded in a Fourier representation, would have a constant amplitude and frequency for each component. The contrast between EMD and Fourier decomposition is clear: the IMF represents a generalized Fourier expansion with a time-varying function for amplitude and frequency. It is worth mentioning that for data coming from non-linear and non-stationary processes, the EMD processing shows the advantage of being implemented in the time domain [8]. The method consists of decomposing the data into a collection of intrinsic mode functions (IMFs) to which the Hilbert analysis can be applied. An IMF is any function with the same number of extrema and zero crossings, with its envelopes, as defined by all the local maxima and minima, being symmetric with respect to zero.

In this way, a decomposition of the data into K ${\mathit{IMF}}_i$ modes is achieved, and a residue r_k(t) which can be either a constant, a monotonic mean trend, or a curve having only one extremum. Finally, the original signal can be represented as the sum of all IMFs and a final residual Eq. (1)

$$x\left(t\right)=\sum _{i=1}^n{\mathit{IMF}}_i\left(t\right)+r_K\left(t\right)$$E15

where ${\mathit{IMF}}_i$ is the i-th $\mathit{IMF}$ K the number of $\mathit{\text{IMFs}}$ and $r_k\left(t\right)$ is the final residual. This frequency–time distribution of the amplitude is designated as the Hilbert spectrum. In Figure 4 we present the $\mathit{\text{IMFs}}$ obtained by EMD considering a SO function. In turn, a Hilbert transform was applied to each $\mathit{IMF}$ to explore the range of frequencies involved. It can be observed that in the first mode frequencies greater than 100 kHz (out of the bandwidth of the signal) appear with low amplitude as expected while in the second mode, a frequency of 60 kHz is noticed inside of the bandwidth with high amplitude.

Figure 4.

Hilbert spectrum and intrinsic modes IMF_i(t) for i = 1, 2 (left upper panel). IMF₂ (t) shows components around 65 kHz distributed in intervals in agreement with the low amplitude super-oscillation signal observed in the SO signal (in red). In the lower panel IMF₃(t).

4. The short time Fourier transform and the continuous wavelet transform (CWT)

One class of time-frequency representations widely used in signal processing is based on the use of time windows, that is, smooth and well-localized functions in an interval. The window frames a portion of the signal and allows the Fourier transform to be applied locally. In this way, temporally localized frequency information in the effective domain of the window is re-elevated. By temporarily shifting the window, the signal domain is covered and the complete time-frequency information of the signal is obtained:

$$G\left(\tau ,\omega \right)=\underset{-\infty }{\overset{+\infty }{\int }}x\left(t\right)g\left(t-\tau \right)exp\left(-\mathit{i\omega t}\right)\mathit{dt}$$E16

Eq. (16) models the temporal shifting of the window g(t) over the signal x(t). Assuming that the window g(t) is well-localized in time in an interval centered at t = 0 with length ∆_t and its Fourier transform is localized in frequency in a band centered at $\omega$ = 0 of width ∆_ω, the shifted and modulated windows given by g(t-τ)exp(−iωt) are well-localized elementary functions in time-frequency, in a rectangle centered at the point (τ,ω) of dimension ∆_t∆_ω.

The set of values G(τ,ω) gives us a map with complete information about the signal in the time-frequency domain. This mapping over the time-frequency domain is known to represent an attractive generalization of the Fourier transform called “short Fourier transform”. It can be understood as a localized treatment of the signal by means of sliding bandpass filters of constant bandwidth.

One of the main problems of this approach is how to determine the optimal window width at which the function g(t,τ) describes correctly both amplitude and frequency properties of the signal x(t). The wavelet approach is essentially an adjustable window spectral analysis. The continuous wavelet transform (CWT) is defined as:

$$W\left(a,b,x,\psi \right)=\underset{-\infty }{\overset{+\infty }{\int }}x\left(t\right)\frac{1}{\sqrt{\left|a\right|}}{\psi }^{\ast }\left(\frac{t-b}{a}\right)\mathit{dt}=x^{\ast }\frac{1}{\sqrt{\left|a\right|}}{\psi }^{\ast }\left(\frac{-b}{a}\right)$$E17

in which $\psi$(.) is a wavelet function called the mother wavelet that satisfies certain conditions, where $\mathrm{a}$ and b are the dilation and the translation parameters respectively. From the above equation $W\left(a,b,x,\psi \right)$ can be interpreted in a simple way as is the energy of x(t) of scale a at t = b.

The CWT uses a window with variable width that depends on the frequency interval under analysis. The window is automatically adjusted depending on frequency. The most widely applied wavelets are Gaussian Wavelet, Morlet Wavelet, and Mexican Hat Wavelet [9].

5. The Gabor transform

The implementation of the Gabor transform for signal processing is efficient when it comes to locating and characterizing events with well-defined, non-overlapping, and relatively long-frequency patterns with respect to the analysis window. The window for this transform is a Gaussian window like:

$$g\left(t,\tau \right)=exp\left(-\frac{{\left(t-\tau \right)}^2}{2{\sigma }^2}\right)$$E18

It is known from the uncertainty principle that ∆_t∆_ω ≥ C. The equality is reached by Gaussian windows. In contrast, it is totally unsuitable for detecting details of short duration. So, SO belong to the group of functions with frequency patterns of short duration and are not good candidates to be processed with Gabor transform with fixed windows. Therefore, an alternative to the Gabor transform is to use modulated windows of variable dimension, adjusted to the oscillation frequency of a given function. Due to this temporal fit of the windows that are called dynamic-adaptive windows [4] this method has been shown to be efficient in localizing frequencies and setting an alternative to the calculation given by Eq. (12) for instantaneous frequencies.

Alternative Gabor transforms, using local adaptive windows, are designed considering maxima and minima of the waveform to be analyzed. We present results for the case of a SO waveform signal given by Eq. (1). The corresponding parameters used are n = 5 and $a$ = 2. To detect SO, two kinds of windows were examined: a fixed Gaussian window and a dynamic adaptive Gaussian. Both can be expressed by Eq. (18) but for the adaptive windows the parameter σ must be introduced as a t function. This function was adjusted in agreement with the position of the maximum and minimum of the function given by Eq. (1). As a result, for the case of dynamic windows we obtain the following curve as it is shown in Figure 5. Minimums fall in the regions where SO function show low amplitude as expected.

Figure 5.

σ(t) parameter function for dynamic adaptive windows obtained considering the time distribution of maximum and minima of the SO function given by Eq. (1).

For the fixed windows, the width was taken as 30 ps. As the window shifts over the SO signal, different spectra are calculated. As an example, for a given shift we show its calculation (Figure 6).

Figure 6.

In the first panel it is shown the window function in red shifting over the SO function in black. In the second panel, in detail. The adaptive window. The third panel shows the Fourier transform of the product between the window function and the signal. The fourth panel shows the two side Fourier transform of the signal with its characteristic frequencies at 20, 60 and 100 kHz.

Shift steps depend on the number of maxima and minima obtained from the SO curve but other points can be taken to increase the resolution like the middle point between maxima and minima time interval. Here, we propose a slight variation in the methodology. To determine the instantaneous frequency, all these spectra were represented in a waterfall scheme showing the evolution of FWHM in each, as can be seen in Figure 7. Frequencies at which the amplitude of normalized spectra is one and a half, were determined in each case, to represent a significant value of the instantaneous frequencies present in the SO signal. Other authors take more complex criteria because of the shape of the spectra [4] by determining the maximum values and the surrounding minima and taking it as the instantaneous frequency.

Figure 7.

Evolution of the spectra obtained as a consequence of the shifting over the superoscillatory function of a dynamically adaptive window. Periodically, as the window moves, the spectrum undergoes narrowing and widenings. To appreciate this effect, at right is shown a front view of the same spectra but in stacked form.

As can be clearly seen, the spectra look symmetric, then frequencies are at each side of the zero frequency, having the same magnitude. Each of the positive frequencies corresponding to FWHM of spectra, were represented together (Figure 7) with the SO signal allowing us to observe the temporal overlapping between the instantaneous frequencies and low amplitude intervals of the SO signal where higher frequencies concentrate.

The evolution of the instantaneous frequency in Figure 8 is a continuous one and shows that over the maximum frequency in the Fourier spectrum there exists a continuum of frequencies determining an extra bandwidth more than a discrete frequency. On the other hand, the described effects can be appreciated in a modified Gabor transform. The Gabor transform is represented by the inverse of its elements and in the logarithmic form to highlight the SO part.

Figure 8.

In blue, the evolution of the instantaneous frequency is compared against the SO function (in red). It can be appreciated that the instantaneous frequencies overpass the greatest frequency present in the spectrum: 100 kHz. On the right-side, Gabor transform is modified by applying an adaptive dynamic window to detect regions with superoscillations. Near zero and other intervals in cyan indicate regions with high frequencies.

After the previous analysis and the same procedure, we present, in addition, results corresponding to a fixed window applied to the same SO signal. In cyan over the SO regions it can be seen frequencies of 20 and 100 kHz but with a fixed window, frequencies higher than 100 kHz do not appear. This fact is reflected in the behavior of the instantaneous frequency because in the region of interest (i.e., around t = 0 s.) frequencies are around 34 kHz (Figure 9).

Figure 9.

On the left side, fixed Gabor transforms when a fixed Gaussian window is applied to detect regions with super oscillations. Near t = 0 s and other regularly spaced time intervals, in cyan, are indicated regions with higher frequencies. At the right side, the instantaneous frequency shows a lower span from 26 to 34 kHz.

6. Discrete wavelet transform (DWT)

The variable time-frequency resolution is an important property of the CWT which permits its use for the analysis of the signals consisting of the slowly varying low-frequency components and the rapidly varying high-frequency components.

When the dilation parameter a and the translation parameter b are discrete and take the form $a_j=2^{-j}$ and $b_{\mathit{jk}}=k{2}^{-j}$ (see Eq. (17)) where $k$ and $l$ are integers, the expression for the CWT is a series where the dyadic coefficients $d\left(j,k\right)$correspond to the DWT of $x\left(t\right)$ [10].

$$x\left(t\right)=\sum _j\sum _{k}d\left(j,k\right)2^{-j/2}\psi \left(2^{-j}t-k\right)$$E19

The comparison between CWT and DWT shows that the signal x(t) in both cases is expressed in terms of dilations and translations of a single mother wavelet $\psi$(t). DWT is a powerful tool to analyze time series. Likewise, it is easy to compute and broadly used wavelet transform, also called dyadic wavelet transform. DWT downsamples the signal and does not have the shift-invariant property [11]. In contrast to CWT, it is non-redundant.

In this case the family $\left\{{\psi }_{\mathit{jk}}\left(t\right)\right\}$ where ${\psi }_{\mathit{jk}}\left(t\right)=2^{-j/2}\psi \left(2^{-j}t-k\right)$ for j and k integers, is an orthonormal wavelet basis, thus, the coefficients are the inner product $d\left(j,k\right)=⟨x\left(t\right),{\psi }_{\mathit{jk}}\left(t\right)⟩$ and consequently, the formula is exact. Applying DWT is equivalent to obtaining a decomposition with passband filters in each octave characterizing the behavior in time:

$$x\left(t\right)=\sum _j{w}_j\left(t\right)$$E20

$$w_j\left(t\right)=\sum _{k}d\left(j,k\right){\psi }_{\mathit{jk}}\left(t\right)\text{for each}\mathrm{j}$$E21

In signal processing, several properties to choose wavelets are important to have good time-frequency localization characteristics: number of vanishing moments, support length, regularity, symmetry, and orthogonality.

Many of the wavelets used in signal processing applications have compact support, which means that the wavelet function is zero except for an interval. The support length represents the length of the filter, and the longer the support is, the larger cost of computation is required.

On the other hand, the vanishing moments of any function are a measure of how that function decays toward infinity. This means the more vanishing moments the wavelet has, the smoother it will be. Therefore, the support length and the vanishing moments must be compromised.

The symlet6 orthogonal wavelet family was used in our discrete transforms. It is shown in Figure 10. Symlets are called after symmetrical wavelets, though they are nearly but not exactly symmetrical. For this reason, they are also known as the Daubechies least asymmetric (LA) wavelets. Among other good properties, they are compactly supported wavelets and they have n/2 vanishing moments (n is the size of the filter). It is the highest number of vanishing moments for a given support width.

Figure 10.

Symlets6. In the first panel, the scaling and wavelet functions. In the second and third panels, the decomposition and reconstruction filters respectively. At the bottom the Fourier spectrums of both, the low pass and high pass filters.

7. Multiresolution analysis

The theory that allowed us to obtain wavelet transforms is explained starting from Multiresolution Analysis (MRA). It consists of a sequence of nested spaces of functions ${\mathrm{V}}_{\mathrm{j}}$ that covers the space of finite energy functions (${\cup }_j{V}_j=L^2\left(R\right)$).

Moreover, as the wavelet family $\left\{{\psi }_{\mathit{jk}}\left(t\right)\right\}$ is an orthonormal basis for$L^2\left(R\right)$, each function of finite energy can be expressed as a sum ${\sum }_l{w}_l\left(t\right)$ where $w_l\in W_l$ are mutually orthogonal.

Below is a brief description of these function spaces (see [12, 13] for details).

For each value of the scale parameter j, the wavelets ${\psi }_{j,k}\left(t\right)$ generate a subspace of signals sharing the same localized octave, or the same frequency rank. We can denote these wavelets subspaces with $W_j$. They are orthogonal and contain all the signals $w_j\left(t\right)$ (see Eq. (21)), where the coefficients

$$\sum _{k}d{\left(j,k\right)}^2<\infty$$E22

Taking the union of all these subspaces, the subspaces V_j are obtained:

$$V_j=\bigcup _{l=j}^{\infty }{W}_l$$E23

They are called the scaling subspaces and the following relation is satisfied

$$v_j\left(t\right)=\sum _{l=j}^{\infty }{w}_l\left(t\right)$$E24

where $w_l\left(t\right)$ is in $W_l$. They contain the signals obtained by superposition. Consequently, the following relation is satisfied:

$$v_j\left(t\right)=v_{j+1}\left(t\right)+w_{j+1}\left(t\right)$$E25

for each scale j. This means the information of $v_j\left(t\right)$ is decomposed in the component $w_{j+1}\left(t\right)$ with the details corresponding to high frequencies and the component $v_{j+1}\left(t\right)$ which represents the trend associated with low frequencies. The components $v_{j+1}\left(t\right)$ and $w_{j+1}\left(t\right)$ are just the orthogonal projections of the signal $x\left(t\right)$ onto the subspaces $V_{j+1}$ and $W_{j+1}$, respectively. The decomposition process can be continued:

$$v_j\left(t\right)=v_{j+2}\left(t\right)+w_{j+2}\left(t\right)+w_{j+1}\left(t\right)$$E26

Considering the decomposition until level (or scale J), we have

$$v_j\left(t\right)=v_J\left(t\right)+w_J\left(t\right)+w_{J-1}\left(t\right)+\dots +w_{j+1}\left(t\right)$$E27

In each step, the details and approximation corresponding to scale j are added, and finally, the signal $x\left(t\right)$ can be decomposed as a sum of projections, for each J, we have

$$x\left(t\right)=v_J\left(t\right)+\sum _{j=1}^J{w}_j\left(t\right)$$E28

In this way, the decomposition in successive projections can be expressed as:

$$x\left(t\right)=\sum _{k}a\left(J,k\right){\phi }_{\mathit{Jk}}\left(t\right)+\sum _{j=l}^J\sum _{k}d\left(j,k\right){\psi }_{\mathit{jk}}\left(t\right)$$E29

Where $d\left(j,k\right)$ are the detail coefficients and $a\left(J,k\right)$ are the scaling coefficients corresponding to scale J.

This embedding of subspaces, where $W_{j+1}$ is the orthogonal complement of $V_{j+1}$ in $V_j$ i.e. $V_j=V_{j+1}\oplus W_{j+1}$ constitutes a multiresolution analysis. In the context of a multiresolution analysis (MRA) the expression for the signal decomposition in scale J [11, 14] is:

$$x\left(t\right)=A_j\left(t\right)+\sum _{j=1}^J{D}_j\left(t\right)$$E30

where $A_j\left(t\right)$ corresponds to the approximation (in $V_j$ subspace) and J details $D_j\left(t\right)$ (in $W_j$ subspaces) which contains a time series related to variations in $x\left(t\right)$ at a certain scale.

The decomposition of the signal at different scales is obtained through the pyramidal algorithm which was introduced by Mallat [9]. This algorithm can be described using both linear filtering operations and matrix manipulations. A wavelet filter has some properties: it must sum up to zero, must have unit energy and must be orthogonal to its even shifts. This is the orthonormality property of wavelet filters. The scaling filter is called the father wavelet filter while the wavelet filter is called the mother wavelet filter.

Figure 11 shows multiresolution analysis of level 12 for sym6. The signal and details $D_j\left(t\right)$ for j = 9, 10, 11, 12 and the approximation A₁₂(t) are plotted. Each detail $D_j\left(t\right)$ has a sample mean of zero, while the sample mean of the smooth A₁₂(t) is equal to the sample mean of $x\left(t\right)$.

Figure 11.

DWT decomposition using sym6. In the first panel the signal x(t) is shown. In the second and third panel the details for scales j = 9 D9(t) and j = 10, D10(t) are shown. At right in the first panel and second panel the details D11(t) and D12(t) are shown, corresponding to scales j = 11, 12. In the third row the approximation A12(t) is shown.

It is interesting to observe the Fourier spectrum of the DWT D9. While for $x\left(t\right)$ only three frequencies appear, as can be seen in the inset of Figure 12 for D9 (curve in blue) a new frequency emerges at 130000 Hz. Fourier transforms of the detail appear to enhance the SO content of $x\left(t\right)$.

Figure 12.

DWT decomposition using sym6 wavelet. In red is the SO signal. In blue the detail D9(t), and in the box D9(t) Fourier spectrum, showing a strong component greater than frequencies present in the Fourier spectrum of the original signal.

There are other discrete transforms defined, whose behavior is similar to DWT because it yields decomposition with the property of being shift-invariant for a time series.

The other Discrete Wavelet transform we used gives an orthogonal scaled base additive decomposition. It returns the projections of the signal $x\left(t\right)$ onto the various wavelet subspaces and final scaling space. This means the original signal can be recovered by adding all the projections.

In Figure 13, it is shown how the SO is partially removed from the original signal $x\left(t\right)$ by subtracting mra9(t), corresponding to the 9th scale component.

Figure 13.

MODWTMRA decomposition using sym6. In the first panel, the SO signal (red) and mra9(t) (blue) are shown. In the second panel, the SO signal minus mra9(t) shows that it is capable of “filtrate” the superoscillation.

8. Filtering methods to separate super oscillations

Filtering systems can be also applied to analyze SO. Multiresolution analysis allows separating the SO part of a signal and therefore localizing the oscillations in time, that is, local frequencies. This separation takes place not only in the frequency domain but also in the time domain.

Multiresolution analysis (MRA), recovers the original version when added back together. Each division of the signal represents its variability and can be associated with a physical meaning. MRA is based on the analysis of wavelets. Additionally, MRA can be based on wavelet packets, or non-wavelet techniques. Real-world signals are a mixture of different components, and only a subset of these components is interesting. Multiresolution analysis allows different resolutions for the separated components of the signals.

8.1 Filter banks

When filter banks are applied to a signal it effectively allows it to separate or remove part of the signal. Our result shows that after the application of this filtering, the SO signal is removed.

Filtering can be described as the creation of a continuous wavelet transform (CWT) or filter bank. Signals to be filtered are sampled with 1024 samples whose peak magnitudes are approximately.

Our parameters considering the time axis to design the filter were

T1 = −5e-6; T2 = 5e-6; T3 = 2e-5; T4 = 3e-5; T5 = 4.5e-5; T6 = 5.5e-5; T7 = 7e-5; T8 = 8e-5; F1 = 105e3; F2 = 190e3;

Here [F1, F2] is the region of the spectra where the SO band lies. The filters are normalized so that the peak magnitude is two for all passbands. The highest-frequency passband is designed so that the magnitude falls to half the peak value at the Nyquist frequency (Figure 14).

Figure 14.

In the upper panel at the left is the original SO signal. In the lower panel at the left signal after filtering. Components were removed in the range of 110–150 kKHz using CWT continuous wavelet transform. At the right panel spectrograph after filtering the signal. Removed components can be observed in the range of 110–150 kHz.

9. Conclusions

We proposed a method to determine the extra-bandwidth of a SO function by determining its temporal variation according to the variation of the FWHM. This value proved to be an efficient parameter to estimate the width of the dynamic-adaptive window. Even so, the influence of the width of the window in the determination of errors or uncertainties of the superoscillatory band has not been initiated at this stage of the investigation. The FWHM is, however, a parameter of the spectrum but other levels are possible and in the future its influence will be explored, for example in signals contaminated with noise, using this criterion. Likewise, other wavelet-based decompositions -including wavelet packets and dual-tree complex wavelet transforms- and the synchrosqueezing transform (SST) [15] could be applied to analyze instantaneous frequencies in SO, and might provide complementary insights.