DDr.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science Extended Fourier analysis of signals
Abstract.
This summary of the doctoral thesis [8] is created to emphasize the close connection of the proposed spectral analysis method with the Discrete Fourier Transform (DFT), the most extensively studied and frequently used approach in the history of signal processing. It is shown that in a typical application case, where uniform data readings are transformed to the same number of uniformly spaced frequencies, the results of the classical DFT and proposed approach coincide. The difference in performance appears when the length of the DFT is selected greater than the length of the data. The DFT solves the unknown data problem by padding readings with zeros up to the length of the DFT, while the proposed Extended DFT (EDFT) deals with this situation in a different way, it uses the Fourier integral transform as a target and optimizes the transform basis in the extended frequency set without putting such restrictions on the time domain. Consequently, the Inverse DFT (IDFT) applied to the result of EDFT returns not only known readings, but also the extrapolated data, where classical DFT is able to give back just zeros, and higher resolution are achieved at frequencies where the data has been extrapolated successfully. It has been demonstrated that EDFT able to process data with missing readings or gaps inside or even nonuniformly distributed data. Thus, EDFT significantly extends the usability of the DFT based methods, where previously these approaches have been considered as not applicable [10-45] . The EDFT founds the solution in an iterative way and requires repeated calculations to get the adaptive basis, and this makes it numerical complexity much higher compared to DFT. This disadvantage was a serious problem in the 1990s, when the method has been proposed. Fortunately, since then the power of computers has increased so much that nowadays EDFT application could be considered as a real alternative.
A Fourier transform is a powerful tool for signal analysis and representation of a real or complex-valued function of time x ( t ) (hereinafter referred to as the signal) in the frequency domain ๐น(๐) = โซ ๐ฅ(๐ก)๐ โ๐๐๐ก ๐๐ก โโโ , (1.1) ๐ฅ(๐ก) = 12๐ โซ ๐น(๐)๐ ๐๐๐ก ๐๐ โโโ . (1.2) The Fourier transforms orthogonality property provide a basis for the signal selective frequency analysis โซ ๐ โ๐๐ ๐ก ๐ ๐๐๐ก ๐๐ก โโโ = 2๐๐ฟ(๐ โ ๐ ), (2) where ๏ท๏ฌ ๏ท are cyclic frequencies, j is an imaginary number such that j =-1 and ๏ค ( ๏ท - ๏ท ) is the Dirac delta function. Unfortunately, the Fourier transforms calculation according to (1.1) requiring knowledge of the signal x ( t ) as well as performing of integration operation in the infinite time interval. Therefore, for practical evaluation of (1.1) numerically, the observation period and the interval of integration is always limited by some finite value ๏ and the signal is known in the time interval - ๏ / t โค ๏ /
2. The same applies to the Fourier analysis of the signal sampled versions - nonuniformly sampled signal x ( t k ) or uniformly sampled signal x ( kT ) for k = - ๏ฅ ,โฆ, - ๏ฅ . Only a finite length sequence x ( t k ) or x ( kT ), k =0,1,2,โฆ, K -
1, are subject of Fourier analysis, where K is a discrete sequence length, T is sampling period, and the signal observation period is equal to ๏ = t K-1 - t or ๏ = KT.
To avoid aliasing and satisfy the Nyquist limit, uniform sampling of continuous time signals should be performed with the sampling period T โค ๏ฐ / ๏ , where ๏ is the upper cyclic frequency of a signal x ( t ). Although nonuniform sampling has no such a strict limitation on the mean sampling period T s = ๏ / K , in the subsequent analysis we suppose that both sequences, x ( t k ) and x ( kT ), are derived from a band-limited in ๏ signal x ( t ). Let's write the basic expressions of classical and extended Fourier analysis of continuous time signal x ( t ) and its sampled versions x ( t k ) and x ( kT ). r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science โ The formulation of a problem is often more essential than its solution which may be merely a matter of mathematical or experimental skill.โ
Albert Einstein
The classical Fourier analysis dealing with the following finite time Fourier transforms ๐น ๏ (๐) = โซ ๐ฅ(๐ก)๐ โ๐๐๐ก ๐๐ก ๏ /2โ ๏ /2 , (3.1a) ๐น ๏ (๐) = โ ๐ฅ(๐ก ๐ )๐ โ๐๐๐ก ๐ ๐พโ1๐=0 , (3.1b) ๐น ๏ (๐) = โ ๐ฅ(๐๐)๐ โ๐๐๐๐๐พโ1๐=0 , (3.1c) ๐ฅ ๏ (๐ก) = 12๐ โซ ๐น ๏ (๐)๐ ๐๐๐ก ๐๐ ๏ โ ๏ , (3.2) where (3.2) is the inverse Fourier transform obtained from (1.2) for a band-limited in ๏ signal. Transforms (3.1b) and (3.1c) are known as Discrete Time Fourier Transforms (DTFT) of the nonuniformly and uniformly sampled signals. The reconstructed signal x ๏ ( t ) outside the observation period ๏ vanishes quickly reaching values close to zeros. The signal amplitude spectrum is the Fourier transform (3.1) divided by the observation period ๐ ๏ (๐) = 1 ๏ ๐น ๏ (๐). (4) The frequency resolution of the classical Fourier analysis is inversely proportional to the observation period ๏ , thus, the longer interval of signal analysis, the higher resolution is achieved. Obviously, one can get the formula (3.1a) by truncation of infinite integration limits in (1.1) and the DTFT (3.1a) and (3.1b) in a result of replacement of infinite sums by finite ones. This mean, the classical Fourier analysis supposed that the signal outside ๏ is zeros. In other words, the Fourier transform calculation by formulas (3.1) is well justified if applied to time-limited within ๏ signals. On the other hand, a band-limited in ๏ signal cannot be also time-limited and obviously have nonzero values outside ๏๏ฎ Generally, the Fourier analysis results obtained by using the exponential basis tend to the Fourier transform, if ๏โ๏ฅ๏ฌ while in any finite ๏ there may exist another transform basis providing a more accurate estimation of (1.1). The idea of extended Fourier analysis is finding the transform basis, applicable to a band-limited signals registered in the finite time interval ๏ and providing the results as close as possible in terms of the L ๏ฒ -norm (or the Euclidean norm) to the Fourier transform (1.1) defined in the infinite time interval. The formulas for proposed extended Fourier analysis could be written as ๐น ๐ผ (๐) = โซ ๐ฅ(๐ก)๐ผ(๐, ๐ก)๐๐ก ๏ /2โ ๏ /2 , (5.1a) ๐น ๐ผ (๐) = โ ๐ฅ(๐ก ๐ )๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 , (5.1b) ๐น ๐ผ (๐) = โ ๐ฅ(๐๐)๐ผ(๐, ๐๐) ๐พโ1๐=0 , (5.1c) r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science ๐ฅ ๐ผ (๐ก) = 12๐ โซ ๐น ๐ผ (๐)๐ ๐๐๐ก ๏ โ ๏ ๐๐, (5.2) where in general case the transform basis ๏ก ( ๏ท , t ), ๏ก ( ๏ท , t k ) and ๏ก ( ๏ท , kT ) are not equal to the classical ones (3.1). Note that the inverse Fourier transform (5.2) still preserves the exponential basis and the Parseval-Plancherel formula โซ |๐ฅ ๐ผ (๐ก)| ๐๐ก โโโ = โซ |๐น ๐ผ (๐)| ๏ โ ๏ ๐๐ holds for it. To ensure that the results of transforms (5.1) are close to the result of the Fourier transform (1.1) for the signal x ( t ), the following minimum least squares expression will be composed and solved |๐น(๐) โ ๐น ฮฑ (๐)| โ ๐๐๐. (6) Unfortunately, as already stated above, the calculation of F ( ๏ท ) cannot be performed directly for a band-limited signal. So, to compose (6), we should find an adequate substitution. Let's recall that a complex exponent at cyclic frequency ๏ท and with a complex amplitude S( ๏ท ) is defined in the infinite time interval as ๐ฅ(๐ , ๐ก) = ๐(๐ )๐ ๐๐ ๐ก , โโ < ๐ก < โ. (7) The Fourier transform of a signal (7) can be expressed by the Dirac delta function (2) โซ ๐ฅ(๐ , ๐ก)๐ โ๐๐๐ก ๐๐ก โโโ = 2๐๐(๐ )๐ฟ(๐ โ ๐ ). (8) Now, we will use (7) as a signal model with known amplitude spectrum S ( ๏ท ) for frequencies in the range - ๏ โค ๏ท โค ๏ and, in the expression (6), substitute F ( ๏ท ) by the Fourier transform of the signal model (8) and the signals x ( t ), x ( t k ) and x ( kT ) in (5.1) by the signal models (7), correspondingly. Finally, the integral least square error estimators for all the three signal cases get the form โ= โซ |2๐๐(๐ )๐ฟ(๐ โ ๐ ) โ โซ ๐(๐ )๐ ๐๐ ๐ก ๐ผ(๐, ๐ก)๐๐ก ๏ /2โ ๏ /2 | ๐๐ ๏ โ ๏ , (9a) โ= โซ |2๐๐(๐ )๐ฟ(๐ โ ๐ ) โ โ ๐(๐ )๐ ๐๐ ๐ก ๐ ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 | ๐๐ ๏ โ ๏ , (9b) โ= โซ |2๐๐(๐ )๐ฟ(๐ โ ๐ ) โ โ ๐(๐ )๐ ๐๐ ๐๐ ๐ผ(๐, ๐๐) ๐พโ1๐=0 | ๐๐ ๏ โ ๏ . (9c) The solutions of (9) for a definite signal model (7) provide the basis ๏ก ( ๏ท , t ), ๏ก ( ๏ท , t k ) and ๏ก ( ๏ท , kT ) for the extended Fourier transforms (5.1). To control how close the selected signal model amplitudes S ( ๏ท ) are to the signals x ( t ), x ( t k ) and x ( kT ) amplitude spectrum, we will find the formulas for estimate signal amplitude spectrum S ฮฑ ( ๏ท ) in the extended Fourier basis ๏ก ( ๏ท , t ), ๏ก ( ๏ท , t k ) and ๏ก ( ๏ท , kT ). The formula (8) is showing the connection between the signal model Fourier transform and its amplitude spectrum, from where S ( ๏ท ) could be expressed as signal model Fourier transform divided by 2 ๏ฐ๏ค ( ๏ทโ๏ท ). Taking (8) into account, S ฮฑ ( ๏ท ) is calculated as the transforms (5.1) divided by the estimate of 2 ๏ฐ๏ค ( ๏ทโ๏ท ) in the extended Fourier basis, which is determined from (9) in the case ๏=๏ฐ and ๏ท = ๏ท , ๐ ๐ผ (๐) = โซ ๐ฅ(๐ก)๐ผ(๐, ๐ก)๐๐ก ๏ /2โ ๏ /2 โซ ๐ ๐๐๐ก ๐ผ(๐, ๐ก)๐๐ก ๏ /2โ ๏ /2 , (10a) ๐ ๐ผ (๐) = โ ๐ฅ(๐ก ๐ )๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 โ ๐ ๐๐๐ก ๐ ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 , (10b) ๐ ฮฑ (๐) = โ ๐ฅ(๐๐)๐ผ(๐, ๐๐) ๐พโ1๐=0 โ ๐ ๐๐๐๐ ๐ผ(๐, ๐๐)
๐พโ1๐=0 , (10c) r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science and showing that the amplitude spectrum on the frequency ๏ท is estimated as ratio of the signal extended Fourier transform to the transform of exponent with a unit amplitude in the same basis. This is true also for the classical Fourier transform. For example, after substituting exponential basis ๐ผ(๐, ๐ก) = ๐ โ๐๐๐ก in (10a), the denominator becomes equal to ๏ as in formula (4) for the classical Fourier analysis. The values of the denominator in formulas (10) are in inverse ratio to the frequency resolution of the extended Fourier transform. Before finding the extended basis functions for arbitrary S ( ๏ท ), it is reasonable to consider a simple signal model having a rectangular form, S ( ๏ท )=1 for - ๏ โค ๏ท โค ๏ and zeros outside. Then the estimators (9) reduce to โ= โซ |2๐๐ฟ(๐ โ ๐ ) โ โซ ๐ ๐๐ ๐ก ๐ผ(๐, ๐ก)๐๐ก ๏ /2โ ๏ /2 | ๐๐ ๏ โ ๏ , (11a) โ= โซ |2๐๐ฟ(๐ โ ๐ ) โ โ ๐ ๐๐ ๐ก ๐ ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 | ๐๐ ๏ โ ๏ , (11b) โ= โซ |2๐๐ฟ(๐ โ ๐ ) โ โ ๐ ๐๐ ๐๐ ๐ผ(๐, ๐๐) ๐พโ1๐=0 | ๐๐ ๏ โ ๏ . (11c) The solution (11) allows us to establish a relationship between the classical and extended Fourier transforms. In this section the integral least square error estimators (9) and (11) are solved and subsequent analysis of the obtained results is performed to find out only those solutions that can lead to practically realizable algorithms.
The solution of (11a) for continuous time signal x ( t ) is found as a partial derivation ๐โ๐๐ผ(๐,๐) = 0, โ ๏ โค ๐ โค ๏ , and leads to the linear integral equation โซ sin( ๏ (๐ก โ ๐))๐(๐ก โ ๐) ๐ผ(๐, ๐ก)๐๐ก ๏ /2โ ๏ /2 = ๐ โ๐๐๐ . (12) Step by step solution of (12) is given in [4] . Finally, the basis ๏ก ( ๏ท , t ) are obtained by applying a specific function system - a prolate spheroidal wave functions [1] ๏น k ( t ), k =0,1,2,..., and are written as series expansion ฮฑ(๐, ๐ก) = โ ๐ต ๐ (๐) ๏ฌ ๐ ๏น ๐ (๐ก) โ๐=0 . (13) The extended Fourier Transform of continuous time signal x ( t ) are given by ๐น ๐ผ (๐) = โ ๐ต ๐ (๐)๐ ๐โ๐=0 , โ ๏ โค ๐ โค ๏ , (14.1) ๐ฅ ๐ผ (๐ก) = โ ๏น ๐ (๐ก)๐ ๐โ๐=0 , โโ < ๐ก < โ, (14.2) ๐ ๐ผ (๐) = โ ๐ต ๐ (๐)๐ ๐โ๐=0 โ |๐ต ๐ (๐)| , (14.3) r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science where ๐ ๐ = ๏ฌ ๐ โซ ๐ฅ(๐) ๏น ๐ (๐)๐๐ ๏ /2โ ๏ /2 , ๏ฌ ๐ = โซ ๏น ๐2 (๐ก)๐๐ก ๏ /2โ ๏ /2 , ๐ต ๐ (๐) = โ ๐ ๏๏ฌ ๐ ๏ ๏น ๐ (๐ ๏ ๏ )(โ๐) ๐ and the Parseval-Plancherel equality gives โซ |๐ฅ ๐ผ (๐ก)| ๐๐ก โโโ = โซ |๐น ๐ผ (๐)| ๏ โ ๏ ๐๐ = โ |๐ ๐ | . The extended Fourier transform in accordance with (14.1) requesting a calculation of infinite sums, this mean, an infinite quantity of mathematical operations, therefore it's impossible for real world applications. Theoretically, the value of the denominator โ |๐ต ๐ (๐)| in the amplitude spectrum formula (14.3) tends to infinite for K โ๏ฅ , and the extended Fourier transform (14.1) provide a super-resolution - an ability to determine the Fourier transform for the sum of sinusoids or complex exponents, if frequencies of them differ by arbitrary small finite value. In this subsection the minimum least square error estimators (9b,c) and (11b,c) are solved and the extended Fourier transforms for uniformly and nonuniformly sampled complex-valued signals are obtained. The proposed approaches have been developed in articles [5] and [6] , where the derivations for real-valued discrete signals are given.
The following notations are used in the matrix equations: superscripts X -1 , X T , X * and X H denote inverse, transpose, complex conjugate and complex conjugate (Hermitian) transpose of the matrix X ; ./ represents element-by-element division of two matrices with the same size; sum ( X ) means addition of all matrix X elements and the diag ( X ) forms the row vector by extracting the main diagonal elements from quadratic matrix X or it puts the elements of vector X on the main diagonal to form a diagonal matrix. The solutions of (11b,c) can be obtained similarly to (11a), as partial derivatives of ๐โ๐๐ผ(๐,๐ก ๐ ) = 0 and ๐โ๐๐ผ(๐,๐๐) = 0 for l =0,1,2,..., K -1, and leads to the systems of linear equations โ sin( ๏ (๐ก ๐ โ ๐ก ๐ ))๐(๐ก ๐ โ ๐ก ๐ ) ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 = ๐ โ๐๐๐ก ๐ , (15a) โ sin( ๏ (๐ โ ๐)๐)๐(๐ โ ๐)๐ ๐ผ(๐, ๐๐) ๐พโ1๐=0 = ๐ โ๐๐๐๐ . (15b) The solution of (15) in the matrix form is expressed as ๐ ๐ = ๐ โ1 ๐ ๐ , (16) where A ๏ท ( K x1) and E ๏ท ( K x1) are the extended Fourier and the exponential basis. The formulas of Extended Discrete Time Fourier Transform (EDTFT) for signal model S ( ๏ท )=1, - ๏ โค ๏ท โค ๏๏ฌ are derived by substituting of transformation basis (16) into expressions (5) and (10) ๐น ๐ผ (๐) = ๐ฑ๐ โ1 ๐ ๐ , โ ๏ โค ๐ โค ๏ , (17.1) ๐ฅ ๐ผ (๐ก) = ๐ฑ๐ โ1 ๐ ๐ก , โโ < ๐ก < โ, (17.2) ๐ ๐ผ (๐) = ๐ฑ๐ โ1 ๐ ๐ ๐ ๐๐ป ๐ โ1 ๐ ๐ . (17.3) The matrices for nonuniformly sampled signal x ( t k ) are composed as follows x (1x K ): x ( t k ), E ๏ท ( K x1): ๐ โ๐๐๐ก ๐ , R ( K x K ): ๐ ๐,๐ = sin( ๏ (๐ก ๐ โ๐ก ๐ ))๐(๐ก ๐ โ๐ก ๐ ) , E t ( K x1): ๐ ๐ = sin( ๏ (๐กโ๐ก ๐ ))๐(๐กโ๐ก ๐ ) . Uniformly sampled sequence x ( kT ) could be considered as a special case of nonuniform sampling at time moments t k = kT , k =0,1,2,โฆ, K -
1, then the matrices in (16, 17) are formed as x (1x K ): x ( kT ), E ๏ท ( K x1): ๐ โ๐๐๐๐ , R ( K x K ): ๐ ๐,๐ = sin( ๏ (๐โ๐)๐)๐(๐โ๐)๐ , E t ( K x1): ๐ ๐ = sin( ๏ (๐กโ๐๐))๐(๐กโ๐๐) . r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science If sampling of signal x ( kT ) is done with Nyquist rate, T =๏ฐ๏ฏ๏ , the matrix R becomes a unit matrix I and the formula (17.1) coincide with classical DTFT (3.1c), but the formula (17.3) reduces to the well-known relationship between discrete signal Fourier transform and its amplitude spectrum ๐น ๐ผ (๐) = ๐น ๏ (๐) = ๐ฑ๐ ๐ , (18.1) ๐ ๐ผ (๐) = 1๐พ ๐ฑ๐ ๐ . (18.2) Whereas for nonuniformly sampled signal x ( t k ) the matrix R ๏น I , even if mean sampling period T s =๏ฐ๏ฏ๏ and formulas (17) give the results that are close to uniform case and superior to those obtained by the classical nonuniform DTFT (3.1b). The resolution by frequency in both sampling cases equals to 1 / KT , which is a normal frequency resolution. While for oversampled signals, T (or T s ) ๏ผ ๏ฐ๏ฏ๏๏ฌ the EDTFT approach can provide a high frequency resolution and improved spectral estimation quality. Unfortunate an achievement of such results is limited by finite precision in the mathematical calculations and by restrictions on frequency range in the process of signal sampling. The theoretical value of the denominator in (17.3) ๐ ๐๐ป ๐ โ1 ๐ ๐ = ๐พ and the frequency resolution should increase proportionally to the number of samples in the signal observation period ๏ . In the border-case, if the number of samples within ๏ is increasing infinity, K โ๏ฅ , and the discrete time signal tends to the continuous time signal x ( t ), the EDTFT (17.1) gives the same result as (14.1). Now, we will consider the solution of the minimum least square error estimators (9b,c) for arbitrary selected signal model S ( ๏ท ). The derivation formulas for both estimators are like the ones given in the previous section. For example, a partial derivation of (9b) by the basis, ๐โ๐๐ผ(๐,๐ก ๐ ) = 0 for l =0,1,2,..., K -
1, provides the least square solution โซ (2๐๐(๐ )๐ฟ(๐ โ ๐ ) โ โ ๐(๐ )๐ ๐๐ ๐ก ๐ ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 ) ๐ โ (๐ )๐ โ๐๐ ๐ก ๐ ๐๐ = 0 ๏ โ ๏ , (19) Equation (19) can be rewritten as โ (โซ |๐(๐ )| ๐ ๐๐ (๐ก ๐ โ๐ก ๐ ) ๐๐ ๏ โ ๏ ) ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 = 2๐ โซ |๐(๐ )| ๐ โ๐๐ ๐ก ๐ ๐ฟ(๐ โ ๐ )๐๐ ๏ โ ๏ . (20) The filtering feature of the Dirac delta function โซ ๐(๐ฅ)๐ฟ(๐ฅ โ ๐ฅ )๐๐ฅ = ๐(๐ฅ ) โโโ applied to the right part of (20) gives the final form of the system of linear equations โ ( 12๐ โซ |๐(๐ )| ๐ ๐๐ (๐ก ๐ โ๐ก ๐ ) ๐๐ ๏ โ ๏ ) ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 = |๐(๐)| ๐ โ๐๐๐ก ๐ , (21a) โ ( 12๐ โซ |๐(๐ )| ๐ ๐๐ (๐โ๐)๐ ๐๐ ๏ โ ๏ ) ๐ผ(๐, ๐ก ๐ ) ๐พโ1๐=0 = |๐(๐)| ๐ โ๐๐๐๐ , (21b) where |๐(๐)| is the signal model power at ๏ท =๏ท . The system of linear equations (21b) is applicable for uniformly sampled signal x ( kT ) and can be derived from (9c) in a similar way as (21a). The EDTFT basis A ๏ท ( K x1) - ๏ก ( ๏ท , t k ) or ๏ก ( ๏ท , kT ) are found as a solution of (21) ๐ ๐ = |๐(๐)| ๐ โ1 ๐ ๐ . (22) Substituting of basis (22) into expressions (5) and (10), yields the formulas for calculation of the EDTFT in the general case ๐น ๐ผ (๐) = ๐ฑ๐ ๐ = |๐(๐)| ๐ฑ๐ โ1 ๐ ๐ , โ ๏ โค ๐ โค ๏ , (23.1) ๐ฅ ๐ผ (๐ก) = ๐ฑ๐ โ1 ๐ ๐ก , โโ < ๐ก < โ, (23.2) ๐ ๐ผ (๐) = ๐ฑ๐ ๐ ๐ ๐๐ป ๐ ๐ = ๐ฑ|๐(๐)| ๐ โ1 ๐ ๐ ๐ ๐๐ป |๐(๐)| ๐ โ1 ๐ ๐ = ๐ฑ๐ โ1 ๐ ๐ ๐ ๐๐ป ๐ โ1 ๐ ๐ . (23.3) r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science The elements of the matrices R ( K x K ) and E t ( K x1) in the formulas (22, 23) are expressed by integrals ๐ ๐,๐ = 12๐ โซ |๐(๐ )| ๐ ๐๐ (๐ก ๐ โ๐ก ๐ ) ๐๐ ๏ โ ๏ or ๐ ๐,๐ = 12๐ โซ |๐(๐ )| ๐ ๐๐ (๐โ๐)๐ ๐๐ ๏ โ ๏ , (24.1) ๐ ๐ = 12๐ โซ |๐(๐)| ๐ ๐๐(๐กโ๐ก ๐ ) ๐๐ ๏ โ ๏ or ๐ ๐ = 12๐ โซ |๐(๐)| ๐ ๐๐(๐กโ๐๐) ๐๐ ๏ โ ๏ , (24.2) for nonuniformly or uniformly sampled signal cases, respectively. If the signal and its model power spectra are close, |๐ ๐ผ (๐ )| โ |๐(๐ )| , then (24.1) is also an estimate of the autocorrelation function of the sequence x . The inverse transform (23.2) calculated on time moments t=t k or t=kT , k =0,1,2,โฆ, K -
1, returns back the input sequence x undistorted, as the elements of matrices E t become equal to R . Case signal model S ( ๏ท )=1 the formulas (22) and (23) reduce to (16) and (17). The frequency resolution of the EDTFT is in inverse ration to |๐(๐)| ๐ ๐๐ป ๐ โ1 ๐ ๐ and varied in the frequency range - ๏ โค ๏ท โค ๏ . Calculation of the EDTFT by formulas (23) requires knowledge of the signal model spectrum which generally is not known. At the same time, the amplitude spectrum obtained in the previous section by the formula (17.3) can be used as a source of such information. This suggests the following iterative algorithm introduced in [5] , where the signal model spectrum S ( ๏ท ) tends to the signal spectrum S ฮฑ ( ๏ท ): Iteration 1:
Calculate ๐ ๐ผ(1) (๐) (17.3) applying default signal model - S ( ๏ท )=1. Iteration 2:
Calculate ๐ ๐ผ(2) (๐) (23.3) by using the signal model - ๐ ๐ผ(1) (๐ ) . Iteration 3:
Calculate ๐ ๐ผ(3) (๐) (23.3) by using the signal model - ๐ ๐ผ(2) (๐ ) . โฆ Iteration I: Calculate ๐ ๐ผ(๐) (๐) (23.3) by using the signal model - ๐ ๐ผ(๐โ1) (๐ ) . The iterations are repeated until the given maximum iteration number is reached or the power spectrum does not alter from iteration to iteration - |๐ ๐ผ(๐) (๐)| โ |๐ ๐ผ(๐โ1) (๐)| . The EDTFT output F ฮฑ ( ๏ท ) (23.1) is calculated for the last performed iteration I . By default, the signal model S ( ๏ท )=1 is used as input for the EDTFT algorithm. However, additional information about the signal to be analyzed can be applied to create a more realistic signal model for the EDTFT input and to reduce the number of iterations required to reach the stopping iteration criteria. EDTFT considered in the previous section is a function of the continuous frequency ( - ๏ โค ๏ท โค ๏ ), while describing below EDFT algorithm calculate EDTFT on a discrete frequency set - ๏๏ฃ๏ท n ๏ผ๏ for n =0,1,2,โฆ, N -
1. The number of frequency points N ๏ณ K and it should be selected sufficiently great to substitute the integrals (24.1) used for calculation of the matrix R ( K x K ) in the expressions (22, 23) by the finite sums ๐ ๐,๐ = 12๐ โซ |๐(๐ )| ๐ ๐๐ (๐ก ๐ โ๐ก ๐ ) ๐๐ ๏ โ ๏ โ ๏ ๐๐ โ |๐(๐ ๐ )| ๐ ๐๐ ๐ (๐ก ๐ โ๐ก ๐ )๐โ1๐=0 , (25.1) ๐ ๐,๐ = 12๐ โซ |๐(๐ )| ๐ ๐๐ (๐โ๐)๐ ๐๐ ๏ โ ๏ โ ๏ ๐๐ โ |๐(๐ ๐ )| ๐ ๐๐ ๐ (๐โ๐)๐๐โ1๐=0 , (25.2) r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science l , k =0,1,2,โฆ, K -
1. The matrices composed of (25.1) and (25.2),
๐ = [ ๐ (0) ๐ (๐ก โ ๐ก )๐ (๐ก โ ๐ก ) ๐ (0) โฏ ๐ (๐ก ๐พโ1 โ ๐ก )๐ (๐ก ๐พโ1 โ ๐ก )โฎ โฑ โฎ๐ ๐พโ1,0 (๐ก โ ๐ก ๐พโ1 ) ๐
๐พโ1,1 (๐ก โ ๐ก ๐พโ1 ) โฏ ๐
๐พโ1,๐พโ1 (0) ] , (26.1)
๐ = [ ๐ (0) ๐ (๐)๐ (โ๐) ๐ (0) โฏ ๐ ((๐พ โ 1)๐)๐ ((๐พ โ 2)๐)โฎ โฑ โฎ๐
๐พโ1,0 (โ(๐พ โ 1)๐) ๐
๐พโ1,1 (โ(๐พ โ 2)๐) โฏ ๐
๐พโ1,๐พโ1 (0) ] , (26.2) possess Hermitian symmetry, ๐ ๐,๐ = ๐ ๐,๐โ , but (26.2) for a uniformly sampled signal has also a Toeplitz structure. The matrix elements ๐ ๐,๐ represents the autocorrelation function and can be calculated by applying the IDFT to the signal model power spectrum |๐(๐ ๐ )| . The frequency ๏๏ฏ๏ฐ=๏ฒ f u = f N in (25), where f u is the signal upper frequency and f N is the Nyquist rate of a band-limited signal, and it is assumed to be normalized (equal to 1) in DFT calculations. The choice of the frequencies { ๏ท n }={2 ๏ฐ f n } depends on the number of frequencies needed for accurate estimation of (25) as well as for detailed signal spectrum representation, and the limitations on the total amount of computation. Eventually, the uniform set of frequencies in range [ - ๏๏ป๏๏ is preferable in most application cases. The EDFT may be expressed by the iterative algorithm ๐ (๐) = 1๐ ๐๐ (๐) ๐ ๐ป , (27.1) ๐ (๐) = ๐ฑ๐ (๐) = ๐ฑ(๐ (๐) ) โ1 ๐๐ (๐) , (27.2) ๐ (๐) = ๐ฑ(๐ (๐) ) โ1 ๐.๐๐๐๐(๐ ๐ป (๐ (๐) ) โ1 ๐), (27.3) ๐ (๐+1) = ๐๐๐๐ (|๐ (๐) | ), (27.4) for iteration number i =1,2,3,โฆ I , wherein (27.1) is the sum (25) in matrix form. The matrix E ( K x N ) has elements ๐ โ๐2๐๐ ๐ ๐ก ๐ or ๐ โ๐2๐๐ ๐ ๐๐ case sampling of x done uniformly. By default, the diagonal weight matrix W ( i ) ( N x N ) for the first iteration is a unit matrix, W (1) = I . If the other diagonal matrix is used as input for the EDFT algorithm, it should have at least K non-zero elements. For the next iterations W ( i +1) is filled with power spectrum values calculated by (27.4). There could be additional criteria for stopping the iterations before the maximum number of iterations I is reached, for example, the iterations could be interrupted, if the relative change in the power spectrum, | sum ( W ( i +1) ) - sum ( W ( i ) ) |/ sum ( W (2) ) for i >1, is smaller than a given threshold. The IDFT may be applied to output F of each iteration and returns original K samples of uniform or nonuniform sequence ๐ฑ = 1๐ ๐ ๐ ๐ป . (28) Since the length of the frequency set N ๏ณ K , then (28) can be modified to obtain an extrapolated sequence ๐ฑ ๐ผ (1x N ) - x ฮฑ ( t m ) or x ฮฑ ( mT ), m =0,1,2,โฆ, N - ๐ฑ ๐ผ = 1๐ ๐ ๐ ๐๐ป , (29) where exponents matrix E N ( N x N ) has elements ๐ โ๐2๐๐ ๐ ๐ก ๐ or ๐ โ๐2๐๐ ๐ ๐๐ case of uniform ๐ฑ ๐ผ , and ๐ฑ๐ฑ ๐ป โค ๐ฑ ๐ผ ๐ฑ ๐ผ๐ป = ๐ ๐ ๐ป according to Parseval-Plancherel theorem. Reconstructed by the formula (29) sequence is the original sequence plus forward and backward extrapolation of x to length N and/or interpolation if there are gaps inside of x . The maximum frequency resolution is limited by the length N of frequency set, not by the length K of sequence as in the application of classical DFT. It means, the EDFT can increase the frequency resolution N / K times in comparison with the r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science classical DFT. This can be verified by comparing the diagonal elements of the product of IDFT and DFT basis, ๐๐๐๐ ( ๐ ๐ป ๐) , which are equal to K / N at all frequencies, with the relationship, ๐ ๐ป ๐) = ๐ ./๐ โค 1 , corresponding to the IDFT and EDFT basis A (27.2). At the same time there is a restriction on the frequency resolution sum ( F . /S )= NK , which is satisfied by iteration, and to achieve a high resolution at certain frequencies, the EDFT must decrease the resolution on other frequencies. The deviation | sum ( F . /S ) - NK | also could be used as an additional criterion for stopping of iterations, because indicates the possible inaccuracy in the obtained results, mainly caused by the finite precision in calculations. If this happens, the result of the previous EDFT iteration should be considered as a final one. In a border-case N = K , the iterative algorithm output does not depend on weight matrix W and the optimal EDFT basis is found in a non-iterative way (in a result of the first iteration) [7] . In the previous sections, starting with the Fourier integral (1) and using its orthogonality property (2), by establishing and solving the minimum least square error estimators (9), the Extended DFT is obtained analytically. In the next, a comparison with known nonparametric approaches - Capon filter, Generalized (Weighted) Least Squares (GWLS) solution and High-Resolution Discrete Fourier Transform (HRDFT) introduced by Sacchi, Ulrych and Walker in 1998 are done, and the ways and opportunities of derivation of an iterative EDFT algorithm based on these methods are analyzed briefly.
The Capon filter known also as Minimum Variance spectral estimate (see [3 , , , ) can be viewed as the output of a bank of filters with each filter centered at one of the analysis frequencies ๐ฆ ๐ (๐๐) = โ ๐ฅ((๐ โ ๐)๐) ๐พโ1๐=0 โ ๐ (๐๐) = ๐ฑฬ๐ก ๐ , ๐ = 0,1,2, โฆ . (30) In the matrix notation ๐ฑฬ = [๐ฅ(๐๐), ๐ฅ((๐ โ 1)๐), โฆ , ๐ฅ((๐ โ ๐พ + 1)๐)] is the filter input signal and ๐ก ๐ = [โ ๐ (0), โ ๐ (๐), โฆ , โ ๐ ((๐พ โ 1)๐)] ๐ is the filter coefficients. Here the subscript ฯ indicate a dependence on the filterโs center frequency. The Capon filter is designed to minimize the variance on the filter output ๐ ๐ฆ2 = ๐บ{|๐ฆ ๐ (๐๐)| } = ๐บ{๐ฆ ๐๐ป (๐๐)๐ฆ ๐ (๐๐)} = ๐บ{๐ก ๐๐ป ๐ฑฬ ๐ป ๐ฑฬ๐ก ๐ }= ๐ก ๐๐ป ๐บ{๐ฑฬ ๐ป ๐ฑฬ}๐ก ๐ = ๐ก ๐๐ป ๐ ๐ฅ ๐ก ๐ , (31) subject to the constraint that its frequency response at the frequency of interest ฯ has unity gain ๐ป(๐) = โ โ ๐ (๐๐)๐ โ๐๐๐๐๐พโ1๐=0 = ๐ ๐๐ ๐ก ๐ = 1, (32.1) ๐ป(๐) = โ โ ๐โ (๐๐)๐ ๐๐๐๐๐พโ1๐=0 = ๐ก ๐๐ป ๐ ๐โ = 1, (32.2) Where ๐บ{. } denotes the expectation operator and the matrix E ๏ท ( K x1) has elements ๐ โ๐๐๐๐ . The constraints (32.1) and (32.2) must be satisfied by filter (30) and Hermitian transpose filter ๐ฆ ๐๐ป (๐๐) = ๐ก ๐๐ป ๐ฑฬ ๐ป , correspondingly. The matrix ๐ = ๐บ{๐ฑฬ ๐ป ๐ฑฬ} ( K x K ) is the sample autocorrelation matrix and it can be composed of the values of the signal autocorrelation function. For example, so called biased estimate is calculated by ๐ ๐ฅ๐ฅ (๐๐) = 1๐พ โ ๐ฅ((๐ + ๐)๐) ๐พโ๐โ1๐=0 ๐ฅ โ (๐๐), ๐ = 0,1,2, โฆ , ๐พ โ 1, (33) and, considering that ๐ ๐ฅ๐ฅ (โ๐๐) = ๐ ๐ฅ๐ฅโ (๐๐) , the sample autocorrelation matrix is filled as r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science ๐ ๐ฅ = [ ๐ (0) ๐ (โ๐)๐ (๐) ๐ (0) โฏ ๐ (โ(๐พ โ 1)๐)๐ (โ(๐พ โ 2)๐)โฎ โฑ โฎ๐ ๐พโ1,0 ((๐พ โ 1)๐) ๐
๐พโ1,1 ((๐พ โ 2)๐) โฏ ๐
๐พโ1,๐พโ1 (0) ] . (34) Mathematically, the Capon filter coefficients can be obtained by minimizing the variance (31) under the constraints given by (32.1) and (32.2)
๐ฝ = ๐ก ๐๐ป ๐ ๐ฅ ๐ก ๐ โ ๐(๐ ๐๐ ๐ก ๐ โ 1) โ ๏ฌ (๐ก ๐๐ป ๐ ๐โ โ 1) โ ๐๐๐, (35) where ๏ญ , ๏ฌ are Lagrange multipliers. The conditions ๐๐ฝ๐๐ก ๐ = 0 and ๐๐ฝ๐๐ก ๐๐ป = 0 must be fulfilled to determine the minimum of (35). Both requirements lead to the same solution ๐ก ๐ = ๐ ๐ฅโ1 ๐ ๐โ ๐ ๐๐ ๐ ๐ฅโ1 ๐ ๐โ (36) and, traditionally, the Capon power spectrum is computed as ๐ ๐ถ๐๐๐๐ (๐) = ๐ก ๐๐ป ๐ ๐ฅ ๐ก ๐ = 1๐ ๐๐ ๐ ๐ฅโ1 ๐ ๐โ , (37) To obtain an iterative EDFT algorithm from the original Capon filter approach, the sample autocorrelation matrix R x (34) must be substituted by R T = E * WE T . The matrix R T ( K x K ) can also be obtained as a transpose of the EDFT matrix R defined by (26). The elements of quadratic diagonal matrix W ( N x N ) represent an estimate of power at time moment nT =0, determined from one sample at output of each Capon filter |๐ฆ ๐ (0)| = |๐ฑฬ๐ก ๐ | = | ๐ฑฬ(๐ ๐ ) โ1 ๐ ๐โ ๐ ๐๐ (๐ ๐ ) โ1 ๐ ๐โ | (38) where the filter input sequence ๐ฑฬ (30) is related to the EDFT input sequence x as ๐ฅฬ(๐๐) = ๐ฅ((๐พ + ๐ โ 1)๐) or ๐ฅฬ(๐ก ๐ ) = ๐ฅ(๐ก ๐พ+๐โ1 ), k =0, - - - ( K - Eventually, an iterative algorithm can be formed as follows ๐ ๐(๐) = ๐ โ ๐ (๐) ๐ ๐ , (39.1) ๐ ๐ถ๐๐๐๐(๐) = ๐ฑฬ(๐
๐(๐) ) โ1 ๐ โ .๐๐๐๐(๐ ๐ (๐ ๐(๐) ) โ1 ๐ โ ), (39.2) ๐ (๐+1) = ๐๐๐๐ (|๐ ๐ถ๐๐๐๐(๐) | ), (39.3) with the initial condition for W (1) = I and the iteration number i =1,2,3,โฆ I . The estimate of the power spectrum |๐ ๐ถ๐๐๐๐(๐) | coincides with the results of the EDFT, while the phase spectrum, definitely, is different. It should be noted that the calculation of the Capon filter output power by (37) is theoretically well justified, whereas the derivation of (39) requires ad hoc assumptions and substitutions, and actually is a measurement of power obtained from just a one sample at the output of the filter. This leads to conclusion that the approach (39) is simply a filter-bank interpretation of the EDFT, similarly to the DFT which can also be considered as a bank of filters. In addition, an iterative algorithm derived based on Capon filter cannot reveal all the EDFT capacity, such as the ability to estimate DFT (27.2) and restore the signal (28, 29). The Generalized (Weighted) Least Squares approach (see [3, 15 , , ) in the spectral analysis could be based on the following data model ๐ฑ ๐ = ๐ ๐โ ๐ ๐บ๐๐ฟ๐ (๐) + ๐ ๐ , (40) r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science with ๐ ๐ denoting the noise and interference (signals at frequencies other than ฯ ) component, and ๐ ๐โ ๐ ๐บ๐๐ฟ๐ (๐) representing the signal component on the frequency of interest with unknown complex amplitude ๐ ๐บ๐๐ฟ๐ (๐) . The GWLS minimizes [๐ฑ ๐ โ ๐ ๐โ ๐ ๐บ๐๐ฟ๐ (๐)] ๐ป ๐ โ1 [๐ฑ ๐ โ ๐ ๐โ ๐ ๐บ๐๐ฟ๐ (๐)], (41) which is solved by ๐ ๐บ๐๐ฟ๐ (๐) = ๐ ๐๐ ๐ โ1 ๐ฑ ๐ ๐ ๐๐ ๐ โ1 ๐ ๐โ , (42) where Q ( K x K ) is the covariance matrix of the data model component ๐ ๐ . There are two special cases of GWLS called Weighted Least Squares (WLS) and ordinary Least Squares (LS). WLS occur when all the off-diagonal entries of Q are 0, while LS solution is obtained from the GWLS under the assumption that ๐ ๐ at (40) is a white noise, hence Q = I . The problem of GWLS estimator is that, in general, the noise covariance matrix Q is not known, and must be estimated from the data along with the S GWLS ( ฯ ). The initial estimate (the 1 st iteration) could be equal to LS solution, it is (42) with Q = I . Next, to ensure that the GWLS solution works in an iterative way as EDFT do, the covariance matrix should be calculated as Q = R T = E * WE T under the assumption ๐ = ๐๐๐๐(|๐
๐บ๐๐ฟ๐ (๐)| ) . In result, GWLS solution (42) coincides with the EDTFT formula (23.3) ๐ ๐บ๐๐ฟ๐ (๐) = ๐ ๐๐ (๐ ๐ ) โ1 ๐ฑ ๐ ๐ ๐๐ (๐ ๐ ) โ1 ๐ ๐โ = ๐ฑ๐ โ1 ๐ ๐ ๐ ๐๐ป ๐ โ1 ๐ ๐ = ๐ ฮฑ (๐). (43) and, as shown in the Section 3.3.3, can be successfully used to update of the amplitude spectrum iteratively. Although substitution of a noise matrix by R T would be easy done, it is not supported by GWLS data model (40), from where the matrix Q represents the data model component ๐ ๐ only and the signal component ๐ ๐โ ๐ ๐บ๐๐ฟ๐ (๐) must be excluded from it, whereas the matrix R T is calculated for the entire signal ๐ฑ ๐ , including ๐ ๐ and ๐ ๐โ ๐ ๐บ๐๐ฟ๐ (๐) . Furthermore, the derivation of EDFT shows that the signal is restored by applying IDFT to the Extended Fourier transform (28), not as an inverse of the amplitude spectrum (27.3), which is a scaled version of (27.2) by a frequency dependent weight factor, ๐ฑ ๐ = ๐ โ ๐ โ ๐ โ ๐ . Using an estimate S GWLS ( ๏ท )= S ฮฑ ( ๏ท ) in the data model (40) leads to a predetermined split of overall energy at the frequency ฯ in between both components, where the noise part ๐ ๐ may be expressed as a difference of EDFT outputs F ฮฑ ( ๏ท ) and S ฮฑ ( ๏ท ). The conclusion reached is that making the derivation of the Extended DFT algorithm possible, invalidates GWLS minimization expression (41) which require separation of both data model components. The third method considered here is High-Resolution DFT proposed by Sacchi, Ulrych and Walker in [9] . The authors presented an iterative nonparametric approach of spectral estimation, which minimizes the cost function deduced from Bayesโ theorem and, as well as Extended DFT, makes it possible to obtain high-resolution Fourier spectrum. The HRDFT algorithm can be reduced to the following iterative procedure: ๐ (๐) = 1๐ ๐๐ (๐) ๐ ๐ป , (44.1) ๐ ๐ป๐ ๐ท๐น๐(๐) = ๐ฑ(๐ (๐) ) โ1 ๐๐ (๐) (44.2) ๐ (๐+1) = ๐๐๐๐ (| 1๐ ๐ ๐ป๐ ๐ท๐น๐(๐) | ), (44.3) for iteration number i =1,2,3,โฆ I and with the initial condition W (1) = I . r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science The IDFT (28) applied to any iteration output (44.2) returns back the sequence x undistorted. The main difference between approaches is that the HRDFT algorithm lack of formula to estimate of amplitude spectrum (27.3). Instead, as input for the next iteration, it uses the Fourier spectrum estimated in the previous iteration. Thus, the results of HRDFT differ from output of EDFT significantly. HRDFT iterates to the solution where the signal is approximated by K frequencies while the power on other N - K frequencies becomes negligible. Each valuable frequency is resolved with maximum resolution restricted by the length of HRDFT. Also, it still obeys the same limit on the sum of resolutions by frequency ( KN ) as DFT and EDFT. The authors [39] investigated algorithms having different from (44.3) weights for adaptation of the correlation matrix (44.1), although only the amplitude spectrum (10) derived accordingly to the minimum least squares expression (9) and calculated by (27.4) fits perfectly to an iterative update of the matrix R and returns results that are closest to the Fourier transform in the L -norm sense. The EDFT algorithm is validated on the data which are similar to those that have been used in [5 , , . The true spectrum of the first test signal consists of a band-limited noise (flat) in the frequency range [ - - f u =0.5 Hz. Uniform and nonuniform sequences of the length K =64 samples are derived by simulating 10-bit Analog-to-Digital Converter (ADC). Sampling and mean sampling periods of both sequences are equal to 1 second, T=T s =1s . Sampling time points for the nonuniform sequence are generated as, t k = kT + ๏ด k , k =0,1,2,..., K - where { ๏ด k } are uniformly distributed random values in the range [0...0.8s]. Thus, the true spectrum of sequences consists of three non-overlapping in frequency domain components and ADC added floor noise ( ๏ป - The plots in Figures 1 and 2 show the performance of EDFT (black lines) for uniform and nonuniform sequences and allows to compare it with the classical DFT (blue lines). The number of frequencies (the length of DFT) is chosen equal to N =1000, which gives spectral estimates with DFT frequency bin spacing 2 f u / N =0.001 Hz. This means that the range [ - S | ) in a non-iterative way. The input matrix W in this case is composed of values of the true spectrum (red line in the plots), therefore there is no need for further iterations. Non-iterative estimate is very close to the EDFT 15 th iteration depicted in Figures 1b and 2b, where the matrix W = I used in the input and confirms the correctness of the iterative algorithm. The Figure 1c (2c) shows the Power Spectral Density (PSD) calculated by the EDFT as 10log(| F | / N ) and proves the expectations, that PSD estimate on a complex exponent should increase in a value in comparison with the classical DFT if the proposed method achieves a high resolution around this frequency. Figure 1d and 2d plotting the relative frequency resolution for the EDFT 15 th iteration calculated as ๐ข ๐๐พ ๐ ./๐ (1d) or ๐ข ๐ ๐ ๐พ ๐ ./๐ (2d) in respect to the DFT for which, in accordance with (18), it is simply equal to 1 at all frequencies. The value ๏ฒ f u T = ๏ฒ f u T s =1 and this means that the signal is processed in one Nyquist zone. The DFT is showing a normal frequency resolution, whereas the EDFT have ability to increase the resolution (in plot appears values >1) around the powerful signal components and decrease the resolution (in plot appears values <1) at frequencies where the signal has weak power components. The EDFT is called as a high-resolution method and that is true, but with the following remark - it keeps the same 'summary' resolution as the traditional DFT or, in other words, squares under black r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science and blue curves in Figure 1d (2d) are equal. The maximum frequency resolution is limited by value of division N / K . For example, if K =64 and N =1000, then the EDFT can potentially improve the frequency resolution 1000 / ๏ป
16 times. The peak resolution is achieved on a deterministic signal part - at frequency 0.35 Hz. The resolution on 0.3985 Hz exponent does not reach maximum value because of its frequency is not on EDFT grid (0.001 Hz) and the estimate splits between two adjacent frequency bins, which is known artifact of the DFT analysis.
Figure 1. Uniform complex-value sequence - the estimate of: (a) Power spectrum - True (red), DFT (blue) and non-iterative EDFT (black), (b) Power spectrum - True (red), DFT (blue) and EDFT (15 th iteration, black), (c) Power Spectral Density - True (red), DFT (blue) and EDFT (15 th iteration, black), (d) Relative frequency resolution - DFT (blue) and EDFT (15 th iteration, black). r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science Figure 2. Nonuniform complex-value sequence - the estimate of: (a) Power spectrum - True (red), DFT (blue) and non-iterative EDFT (black), (b) Power spectrum - True (red), DFT (blue) and EDFT (15 th iteration, black), (c) Power Spectral Density - True (red), DFT (blue) and EDFT (15 th iteration, black), (d) Relative frequency resolution - DFT (blue) and EDFT (15 th iteration, black). Pulse signal ([0...0.25] Hz) is processed by EDFT with about the same resolution as DFT ( ๏พ - - r.sc.comp. Vilnis Liepiลลก Email : [email protected]
The extended summary of the doctoral thesis in Computer Science can handle uniform and nonuniform sequences with the same quality, while the efficiency of classical DFT gets worse in case of nonuniform data. Figure 3 explains the difference in performance between uniform and nonuniform inputs, where the spectra of uniform and nonuniform sequences are analyzed in the extended frequency range, [ - N =2000 and f u =1 Hz. This means that the step by frequency remains the same as in the previous plots. The true spectrum of sequences at frequencies above 0.5 Hz consists only of floor noise ( ๏ป - / ( ๏ฒ f u T s )=0.5 and it is half the normal resolution because of analysis is performed in two Nyquist zones. Nevertheless, squares under blue and black plots in Figure 3c are equal to one's depicted in Figure 2d. The maximum increase in the frequency resolution 2000 / ๏ป๏ณ๏ฑ times is achieved on a complex exponent at frequency 0.35 Hz by the EDFT. The EDFT also increases resolution in half to process transient and random signal components with the normal frequency resolution equal to 1, as it is indicated by the red dotted lines in Figure 3c. Hence the conclusion that EDFT can handle nonuniformly sampled signals in Figure 3. The estimates obtained in the extended frequency range: (a) Power spectrum of uniform sequence - True (red), DFT (blue) and EDFT (black), (b) Power spectrum of nonuniform sequence - True (red), DFT (blue) and EDFT (black), (c) Relative frequency resolution of nonuniform sequence - DFT (blue) and EDFT (black). r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science multiple Nyquist zones only if the overall spectrum of band-limited signal components does not exceed one Nyquist zone. Since the spectrum of uniform sequence (red color line in Figure 1) does not cover the entire Nyquist zone the EDFT should be able to handle it with mean sampling period T s greater than T but less than 2 T . The increase of T s could be achieved by skipping samples from the uniform sequence randomly. The resulting sequence is considered as nonuniformly sampled because the distance between adjacent readings become unequal. The power spectra in Figure 4 show an example of the impact of sample skipping on the performance of DFT and EDFT. Input sequences are modeled by removing 16 and 24 samples randomly from the uniform 64-point data and leads to increase T s =64 / T= T s =64 / T= T s can be achieved for pure deterministic signals [7] . Figure 4. The power spectrum - True (red), DFT (blue) and EDFT (black), of the 64-point sequence without losses (a) and with randomly skipped (b) 16 and (c) 24 samples. r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science Let's validate above expectations for signal consisting only of sinusoids in a white noise and generate a real-value sequence of the length K =64 samples as a sum of four sine waves with predefined amplitudes 0.5, 1, 2 and 3, having arbitrary initial phases and randomly selected frequencies on the EDFT grid (0.001 Hz). Moreover, the signal sampling is also performed on the grid with T =1 second lag by choosing 64 time points within NT =1000 second interval on a random basis and results in about a 16-fold increase of the average sampling period, T s = NT / K= Figure 5. Real-value sequence (a) - 64 nonuniform samples, SNR 20 dB,
Amplitude spectrum (b) - True values (red cycles) and estimates by DFT (blue), EDFT (black),
Relative frequency resolution (c) by DFT (blue line) and EDFT (black),
Original (blue cycles) and interpolated sequence (d) by Inverse EDFT - 1000 uniform samples. r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science The true frequencies and amplitudes (red cycles) as well as the amplitude estimates of DFT (blue line) and EDFT (back line) are depicted in the Figure 5b and showed that DFT cannot recognize weaker power sinusoids while the Extended DFT picks up all of them and estimates their amplitudes and phases precisely. The performance difference is explained in the Figure 5c by comparison of the resolution of both DFTs with respect to the normal frequency resolution (equal to one). The relative resolution of the DFT (blue line) is calculated as 1 / ( ๏ฒ f u T s )= K / N= N - K samples at the input of the DFT could be considered as zeroed by the rectangular windows. The relative resolution of the EDFT (black line) is calculated as ๐ข ๐ ๐ ๐พ ๐ ./๐ = ๐ ./๐ and it increases N / K times reaching the value close to one at frequencies of sinusoids. Thus, signal processing with just a normal frequency resolution allows EDFT not only estimate the parameters of the signal components correctly, but also IDFT applied to its output returns a sequence of length N consisting of the original K and N - K interpolated samples (see Figure 5d). It should be noted that only a deterministic part of the signal is interpolated by EDFT, whereas a white Gaussian noise stays localized in time around the sampling points (blue cycles). The next sequence used in the computer simulations is well-known Marple&Kay data set taken from [3] . It is 64-point real sample sequence of a process consisting of two unit power sine waves with frequencies of 0.2 and 0.21 Hz, a third one with a power of 0.1 (20 dB down) at 0.1 Hz and a Figure 6. The power spectrum obtained for Marple&Kay data set by (a) DFT, (b) EDFT, (c) HRDFT. r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science colored noise in the frequency range [0.2โฆ0.5] Hz (see red color lines in Figure 6). The signal upper frequency is f u =0.5 Hz and the length of the DFT is selected N =1000. Only 500 positive frequencies are displayed, because Marple&Kay sequence is a real-valued and negative frequencies, if they are depicted, gives a symmetrical pattern to zero frequency. The Figure 5 shows the power spectra of the DFT, EDFT and HRDFT approaches in a single picture, while separately, these plots have been presented in [5] and [9] . The performance of other well-known spectral analysis methods for Marple&Kay data set could be found in [3] , including Minimum Variance approach, named in the Section 5.1 as a traditional Capon filter (37). The simulation results in the Figure 6a,b demonstrate that the classical DFT and EDFT can evaluate not only the spectrum of sinusoids, but also the shape of continuous spectrum of other signal components, whereas HRDFT on Figure 6c is suitable mostly for the estimation of a line spectrum. The plot in Figure 6a shows that due to limited frequency resolution the classical DFT cannot resolve sine waves at the frequencies 0.2 and 0.21. Although the first EDFT iteration coincides with DFT, in the further iterations EDFT is able to increase the frequency resolution around the powerful signal components and all three sine waves are clearly distinguished after the 15 th iteration in the Figure 6b. All the three DFTs have one common feature - the ability to get back 64 samples of Marple&Kay data set by applying IDFT to the output of each of these methods. Since the length of the DFT is chosen equal to 1000, the inverse transform (29) returns 1000 -
64 additional samples, which are
Figure 7. Marple&Kay sequence (blue) and extrapolated data (black) by inverse (a) DFT, (b) EDFT, (c) HRDFT. r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science plotted in Figure 7 (black). The samples 65, 66, 67, ... are considered as a forward extrapolation, but samples 1000, 999, 998, ... as a backward extrapolation of known 64-sample sequence (blue). Of course, the Marple&Kay sequence outside of giving data set is unknown, and plots on Figure 6 are just three possible versions of its extrapolation. The classical DFT (Fig.7a) suggests that Marple&Kay sequence outside of given 64 samples will be zeros, HRDFT (Fig.7c) shows that the extrapolated data even will increase in power, while EDFT (Fig.7b) expects that the sequence beyond will have approximately the same power, which only gradually decreases in time. Figure 8. White Gaussian noise (blue) and extrapolated data (black) by inverse (a) DFT, (b) EDFT, (c) HRDFT.
At the end of the computer modeling, we check extrapolated sequences obtained at the output of IDFT if Marple&Kay input data are replaced by the same size white Gaussian noise (Figure 8). According to the theory the PSD of white Gaussian noise should be constant (flat) across the entire frequency range and the readings in a such sequence are uncorrelated random variables, therefore they cannot be extrapolated. In practice, because of finite length sequences and pseudo-random generators used in the simulations, the above expectations are satisfied only approximately. The classical DFT, as the case Marple&Kay data illustrated in Figure 7a, yields zeros outside of given 64-point sequence also in Figure 8a, that this time is perfectly consistent with the theory. Extrapolate by the EDFT (Fig.8b) vanish quickly, and this still agrees with the theory if practical considerations are taken into account. The HRDFT (Fig.8c) in contrary to DFT and EDFT extends the white Gaussian noise up to a length of 1000 samples showing a strong correlation in the input sequence and this is very unlikely to be true. r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science Any approach that claims that it is a high frequency resolution method in accordance with the Uncertainty Principle must make certain assumptions about the data outside of the observation period even if by itself it is not able to recover the signal. The advantage of the proposed method over similar ones is that EDFT based on a solution that satisfies the minimum least squares criteria (6), making it an accurate, reliable and stable. Run MATLAB program EDFT_FIG.m available on mathworks.com to recreate the computer simulations presented in this section.
The EDFT package consisting of programs written in a simple MATLAB code and created to demonstrate the Extended DFT capabilities described in the previous sections. Each function contains commented (%) help text section where its syntax, algorithm, usage and features are described. The programs NEDFT.m and the inverse transform INEDFT.m can be applied for uniform or nonuniform input/output data and frequency sets. function [F,S,Stopit]=nedft(X,tk,fn,I,W) % NEDFT Nonuniform Extended Discrete Fourier Transform. % % SYNTAX % a. Mandatory inputs/outputs % F=nedft(X,tk,fn) % Function NEDFT returns discrete Fourier transform F(fn) of input sequence X sampled at % arbitrary selected time moments tk, where the frequencies fn also may be selected arbitrary. % If size of fn is less than sixe of X, inputs X and tk will be truncated. % b. Mandatory and optional inputs/outputs % [F,S,Stopit]=nedft(X,tk,fn,I,W) % I Optional input parameter I can be used for limiting maximum number of iterations. % If I is not specified in input arguments, default value for I is set by parameter faster, % 'Miteration', that is, nedft(X,tk,fn)=nedft(X,tk,fn,Miteration). To complete iterations % the value for 'Miteration' should be decreased. % W Input weight vector W, if specified, override the default values W=ones(size(fn)). % W must have at least one nonzero element. % S The second output argument S represents the Amplitude spectrum. Peak values of abs(S) can be % used for estimate amplitudes of sinusoids in the input sequence X. % Stopit is an informative output parameter. The first row of Stopit showing the number of performed % iteration, the second row indicate breaking of iteration reason and may have the following values: % 0 - Maximum number of iteration performed. % 1 - Sum of outputs division sum(F./S) is not equal to llength(fn)*length(X) within Relative % deviation 'Rdeviat'. The calculations were interrupted because of results could be inaccurate. % 2 - Relative threshold 'Rthresh' reached. To complete iteration process faster, % the value for 'Rthresh' should be increased. % ALGORITHM % Input: % X - input sequence % E - complex exponents matrix (Fourier transform basis) - E=exp(-i*2*pi*tk.'*fn); % I - (optional) number of maximum iteration. % W - (optional) weight vector W. The default values W = ones(size(fn)). % Output F and S for each NEDFT iteration are calculated by following formulas: % 1. R=E*diag(W/length(fn))*E'; % 2. F=W.*(X*inv(R)*E); % S=(X*inv(R)*E)./diag(E'*inv(R)*E).'; % 3. W=S.*conj(S); - the weight vector W for the next iteration. % A special case: if length(X) is equal to length(fn), the NEDFT output do not depend on selected % weight vector W and is calculated in non-iterative way. % TIPS for selection of mandatory NEDFT inputs X(tk) and fn: % 1. Input sequence X(tk) for NEDFT can be sampled uniformly or nonuniformly. Uniform sampling can be % considered as a special case of nonuniform sampling, where tk=[0,1,...,K-1]*T and T is sampling period. % Nonuniform sampling can be realized in many different ways, like as: % - uniform sampling with randomly missed samples (known as sparse data); % - uniform sampling with missed data segments (known as gapped data); % - uniform sampling with jitter: tk=([0,1,...,K-1] + jitter*rand(1,K))*Ts, where value for jitter is % selected in range [0...1[ and Ts is the mean sampling period; % - additive nonuniform sampling: tk=tk-1 + (1+jitter*(rand-0.5))*Ts, k=1,...K-1, t0=0; % - signal dependent sampling, e.g, level-crossing sampling, etc... . % 2. Frequencies for fn can be selected arbitrary. This mean, that user can choose not only the length % of NEDFT (number of frequencies in fn), but also the way how to distribute frequencies along the r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science % frequency axis. On other hand, to get adequate sequence X representation, frequencies fn should be % selected to cover overall range, where the input sequence X spectrum is supposed to be found, % otherwise, in result of NEDFT, all components having spectra outside fn will be incorporated. % Note that fn should contain negative frequencies too, and for a real value X(tk) analysis each % positive frequency in fn should have corresponding negative one, for example % fn=[-ceil((N-1)/2):floor((N-1)/2)]/(N*T), where N=length(fn). % 3. Frequencies for vector fn can be added in any order. Therefore it is possible to combine different % frequency sets in one or just add individual frequencies of interest to fn, e.g, fn=[fn1 fn2 f1 f2], % where fn1 and fn2 are different frequency sets, f1,f2 - specific frequencies. NEDFT outputs will be % calculated accordingly - F(fn)=[F(fn1) F(fn2) F(f1) F(f2)], S(Fn)=[S(Fn1) S(fn2) S(f1) S(f2)]. % FEATURES % 1. NEDFT output F(fn) is the discrete Fourier transform of sequence X(tk). The Power Spectral Density % function of nonuniform sequence X(tk) can be estimated by the following formula: abs(F).^2/(N*Ts). % 2. In general, the function Y=inedft(F,fn,tn) (see attached program) is used to calculate the reconstructed % sequence Y(tn). If frequencies fn are selected on the same grid as used by FFT algorithm, then ifft(F) % can be applied to get uniformly re-sampled and extrapolated to length(fn) version of input sequence X(tk). % 3. NEDFT output S(fn) estimate amplitudes and phases of sinusoidal components in sequence X(tk). % 4. NEDFT can increase frequency resolution length(fn)/length(X) times. Division of outputs 1/(Ts*(F./S)) % demonstrate the frequency resolution of NEDFT. The following is true for any NEDFT iteration: % 0
The extended summary of the doctoral thesis in Computer Science if abs(ERE*W(l,:).'/N/K-1)>Rdeviat,Stopit(:,l)=[it;1];S(l,W(l,:)==0)=zeros(1,sum(W(l,:)==0));break,end; % Calculate weight (W) for the next iteration. W(l,:)=S(l,:).*conj(S(l,:)); % Stopit 2: Break iterations if relative threshold reached. SW(it)=sum(W(l,:)); if it>1,thit=abs(SW(it-1)-SW(it))/SW(1); if thit<=Rthresh,Stopit(:,l)=[it;2];break,end,end end % ... end iterations. end %========================== Adjust size of NEDFT output =========================================== if trf==1,F=F.';S=S.';end % Adjust size of NEDFT outputs. function Y=inedft(F,fn,tn) %INEDFT Inverse Nonuniform Extended Discrete Fourier Transform. % % Y=inedft(F,fn,tn) is the inverse discrete Fourier transform of the vector F estimated % by NEDFT function at arbitrary frequency set fn: % F(fn) -> Y(tn), % where time moments tn for reconstructed sequence Y can be uniformly or nonuniformly % spaced in time. In the special case of uniform vectors fn and tn, the INEDFT function % can be replaced by well-known MATLAB function IFFT. % % If input arguments are matrixes, the INEDFT operation is applied to each column. % % See also IFFT, EDFT, NEDFT. % % Vilnis Liepins ([email protected]) %========================== Check INEDFT input arguments ===================================== if nargin<3,error('Not enough input arguments. See help inedft.'),end if size(F,1)==1, trf=1; F=F.'; tn=tn.'; else trf=0; fn=fn.'; end [N L]=size(F); if size(fn,2)~=N, error('Sizes of input arguments F and fn must be equal. See help inedft.'),end if size(tn,2)~=L, error('Incorrect size of input argument tn. See help inedft.'),end %========================== Calculate INEDFT for each X column l ================================ for l=1:L E=exp(i*2*pi*tn(:,l)*fn(l,:)); Y(:,l)=E*F(:,l)/N; end %========================== Adjust size of INEDFT output ======================================= if trf==1,Y=Y.';end From the viewpoint of calculation complexity is reasonable to use the same frequency grid as Fast Fourier Transform (FFT.m in MATLAB library). This allows to apply the FFT algorithm in EDFT calculations, which considerably reduce computation time, because each FFT requiring number of operations proportional to N log( N ) rather than N [2] , however the efficiency of FFT algorithm could be in between of these two values since it depends also on the value of N . EDFT.m program is designed as a faster realization of Extended DFT, where the algorithm described in [7] is implemented. The code is applicable for uniformly sampled signals or data with gaps in it. The inverse transform to EDFT.m is MATLAB library program IFFT.m. function [F,S,Stopit]=edft(X,N,I,W) % EDFT Extended Discrete Fourier Transform. % % Function EDFT produce discrete N-point Fourier transform F and amplitude spectrum S of the data % vector X. Data X may contain NaN (Not-a-Number). Inverse transform is Matlab built-in function ifft. % % SYNTAX % [F,S,Stopit]=edft(X,N) for N>length(X) calculate F and S iteratively (see ALGORITHM below). % If data X do not contain NaN and N<=length(X)or N is not specified, EDFT return % the same results as fast Fourier transform: F=fft(X,N) or F=fft(X) and S=F/N. % [F,S,Stopit]=edft(X,N,I) performs edft(X,N) with limit I for maximum number of iterations. % Default value for I is set by parameter 'Miteration', that is, edft(X,N)=edft(X,N,Miteration). % To complete iteration process faster, the value for 'Miteration' should be decreased. % [F,S,Stopit]=edft(X,N,I,W) execute edft(X,N,I) with initial conditions defined by weight % vector W. Default values for W are ones(size(F)). W must have at least length(X) nonzero elements. % Stopit is an informative (optional) output parameter. The first row of Stopit showing the number r.sc.comp. Vilnis Liepiลลก
Email : [email protected]
The extended summary of the doctoral thesis in Computer Science % of performed iteration, the second row indicate breaking of iteration reason and may have % the following values: % 0 - Maximum number of iteration performed. If length(X)<=N, only one EDFT iteration % is performed (I=1). % 1 - Sum of outputs division sum(F./S) is not equal to length(X)*N within Relative deviation % 'Rdeviat'. The calculations were interrupted because of results could be inaccurate. % If this occur in the first iteration then outputs F and S are zeros. % 2 - Relative threshold 'Rthresh' reached. To complete iteration process faster, the value % for 'Rthresh' should be increased. % ALGORITHM % Input: % X - input data. % N - length of discrete Fourier transform. % I - (optional) number of maximum iteration. If not specified, I=Miteration. % W - (optional) weight vector W. If not specified, W = ones(1,N) is used for the first iteration. % E - Fourier transform basis matrix: E=exp(-i*2*pi*(0:length(X)-1)'*(0:N-1)/N); % If part of unknown data in X are replaced by NaN then the time vector (0:length(X)-1) % is changed to exclude time moments where NaN inserted. % Output F and S for each EDFT iteration are calculated by following formulas: % 1. R=E*diag(W/N)*E'; % EDFT using function ifft to calculate R faster. % 2. F=W.*(X*inv(R)*E); % S=(X*inv(R)*E)./diag(E'*inv(R)*E).'; % Levinson-Durbin recursion used for inverse of toeplitz R. % Function fft applied to speed up matrix multiplications. % 3. W=S.*conj(S); W used as input to the next EDFT iteration. % A special case: If length(X) is equal to N, the EDFT output do not depend on selected weight % vector W and is calculated in non-iterative way. % FEATURES % 1. EDFT output F is the N-point Fourier transform of data X. % The Power Spectral Density (PSD) function can be calculated by the following formula: % abs(F).^2/(N*T), T - sampling period. % 2. EDFT can extrapolate input data X to length N. That is, if apply EDFT for % N>length(X), get the results: F=edft(X,N)=edft(Y)=fft(Y); Y=ifft(F), where Y is input X % plus non-zero forward and backward extrapolation of X to length N. % 3. EDFT output S estimate amplitudes and phases of sinusoidal components in input data X. % 4. EDFT can increase frequency resolution N/length(X) times. Division of outputs % 1/(T*F./S) demonstrate the frequency resolution of EDFT. The following is true or any EDFT % iteration: 0 The extended summary of the doctoral thesis in Computer Science for l=1:L, %========================== Check for a special cases ========================================= if KK(l)==N|KK(l)==0, % If length(X)=N or X(:,l) has all NaNs then F(:,l)=fft(X(:,l),N); % EDFT output (F,S) equals to FFT. S(:,l)=F(:,l)/N; Stopit(:,l)=[1; 0]; elseif K==1&N~=1, % Special case, the length(X)=1, F(:,l)=fft(X(:,l),N).'; % EDFT output (F,S) equals to FFT. S(:,l)=F(:,l)/N; Stopit(:,l)=[1; 0]; elseif isempty(find(X(:,l)))&KK(l)>0, % If input X(:,l) has all zeros or zeros&NaN Stopit(:,l)=[1; 0]; % then EDFT output (F,S) is zeros. %========================== Basic EDFT algorithm started ====================================== elseif KK(l)==K, % Input X(:,l) does not contain NaN %========================== Apply FASTER algorithm ========================================= for it=1:I, % Start iterations... % Calculate correlation vector (r). r=ifft(W(:,l)); % Perform inverse of correlation matrix R. % a=[1; toeplitz(conj(r(1:K-1)))\(-r(2:K))];V=a.'*conj(r(1:K)); % use Matlab \ operator a=-r(2)/r(1); V=r(1)-r(2)*conj(r(2))/r(1); % Levinson-Durbin recursion. for n=1:K-2, alfa=[1 a.']*r(n+2:-1:2);rho=-alfa/V;V=V+rho*conj(alfa);a=[a+rho*conj(flipud(a));rho]; end a=[1;a]; % Calculate ERE=diag(E'*inv(R)*E) and XR=X*inv(R). XR=zeros(K,1);RE=zeros(K,1);rc=a; for k=1:K/2, k0=K-k+1;k1=2:K-2*k+1; k2=k+1:K-k; k3=k:K-k+1; RE(1)=RE(1)+2*rc(k);RE(k0-k+1)=RE(k0-k+1)+2*rc(k0);RE(k1)=RE(k1)+4*rc(k2); XR(k)=XR(k)+rc(k3)'*X(k3,l);XR(k0)=XR(k0)+(flipud(rc(k3))).'*X(k3,l); XR(k2)=XR(k2)+rc(k2)*X(k,l)+flipud(conj(rc(k2)))*X(k0,l); rc(k2)=rc(k2-1)+conj(a(k+1))*a(k2)-a(k0)*flipud(conj(a(k2+1))); end if round(K/2)>K/2, RE(1)=RE(1)+rc(k+1);XR(k+1)=XR(k+1)+X(k+1,l)*rc(k+1); end ERE=real(fft(RE,N));W(:,l)=W(:,l)/real(V); % Stopit 1: Break iterations if sum(F./S) is not equal to N*K or NaN. stit=abs(ERE.'*W(:,l)/N/K-1);if (stit>Rdeviat)|isnan(stit), Stopit(:,l)=[it-1; 1]; break,end % Calculate outputs for iteration (it): N-point EDFT (F) and Amplitude Spectrum (S). F(:,l)=fft(XR,N); S(:,l)=F(:,l)./ERE; F(:,l)=F(:,l).*W(:,l); % Calculate weight (W) for the next iteration. W(:,l)=S(:,l).*conj(S(:,l)); % Stopit 2: Break iterations if relative threshold reached. SW(it)=sum(W(:,l)); if it>1, thit=abs(SW(it-1)-SW(it))/SW(1);if thit<=Rthresh, Stopit(:,l)=[it; 2]; break,end end end % ... end iterations. %========================== End of FASTER algorithm ========================================= else % Input X(:,l) contains NaN %========================== Apply SLOWER algorithm ========================================= INVR=zeros(K);ER=zeros(K,1); X(find(~Xnan(:,l)),l)=zeros(K-KK(l),1); % Replace NaN by 0 in X t=find(Xnan(:,l)); % Sample numbers vector (t) for it=1:I, % Start iterations... % Calculate correlation matrix R by applying ifft and inverse R by Matlab backslash operator. RT=ifft(W(:,l)); R=toeplitz(RT(1:K)); INVR(t,t)=R(t,t)\eye(KK(l)); % Calculate ERE=diag(E'*inv(R)*E).' by applying fft. ER(1)=trace(INVR); for k=1:K-1 ER(k+1,1)=sum(diag(INVR,k)+conj(diag(INVR,-k))); end ERE=real(fft(ER,N)); % Stopit 1: Break iterations if sum(F./S) is not equal to N*KK or NaN. stit=abs(ERE.'*W(:,l)/N/KK(l)-1);if (stit>Rdeviat)|isnan(stit), Stopit(:,l)=[it-1; 1];break,end % Calculate outputs for iteration (it): N-point EDFT (F) and Amplitude Spectrum (S). F(:,l)=fft(conj(INVR)*X(:,l),N); S(:,l)=F(:,l)./ERE; F(:,l)=F(:,l).*W(:,l); % Calculate weight (W) for the next iteration. r.sc.comp. Vilnis Liepiลลก Email : [email protected] The extended summary of the doctoral thesis in Computer Science W(:,l)=S(:,l).*conj(S(:,l)); % Stopit 2: Break iterations if relative threshold reached. SW(it)=sum(W(:,l)); if it>1, thit=abs(SW(it-1)-SW(it))/SW(1);if thit<=Rthresh,Stopit(:,l)=[it; 2];break,end end end % ... end iterations. %========================== End of SLOWER algorithm ======================================= end end %========================== Adjust size of EDFT output ======================================= if trf==1,F=F.';S=S.';end The next program demonstrates the applicability of the Extended DFT in 2-dimensional signal processing. The EDFT2.m program is based on the MATLAB library program FFT2.m where FFT.m calls are simply replaced by EDFT.m. The inverse transform to EDFT2.m is the MATLAB library program IFFT2.m. function f = EDFT2(x, mrows, ncols) % EDFT2 Two-dimensional Extended Discrete Fourier Transform. % % EDFT2(X) returns the two-dimensional Fourier transform of matrix X. % Before running EDFT2 unknown data (if any) inside of X should be replaced by NaN (Not-a-Number). % If X is a vector, the result will have the same orientation. % EDFT2(X,MROWS,NCOLS) performing size MROWS-by-NCOLS Fourier transform without padding of % matrix X with zeros. Inverse transform is Matlab built-in function ifft2. % % See also EDFT, IFFT2. % % EDFT2 is created on basis of FFT2 (J.N. Little 12/18/1985) by Vilnis Liepins. % No input. if nargin==0, error('Not enough input arguments. See help edft2.'), end [m, n] = size(x); % Basic algorithm. if (nargin == 1) & (m > 1) & (n > 1) % f = fft(fft(x).').'; f = edft(edft(x).').'; return; end % Padding for vector input. if nargin < 3, ncols = n; end if nargin < 2, mrows = m; end mpad = mrows; npad = ncols; if m == 1 & mpad > m, x(2, 1) = 0; m = 2; end if n == 1 & npad > n, x(1, 2) = 0; n = 2; end if m == 1, mpad = npad; npad = 1; end % For row vector. % Transform. %f = fft(x, mpad); %if m > 1 & n > 1, f = fft(f.', npad).'; end f = edft(x, mpad); if m > 1 & n > 1, f = edft(f.', npad).'; end The first version of EDFT (file GDFT.m) was submitted to file-exchange server on 10/7/1997 as MATLAB 4.1 code. The renewed code version uploaded on 8/5/2006 and available online mathworks.com and researchgate.net. Please note that programs have not been tested on the latest MATLAB versions and therefore have opportunities to performance improvements (see for example [25 , ). [1] D. Slepian, H.O. Pollak. Prolate Spheroidal Wave Functions, Fourier Analysis and Uncertainty - I. Bell System Technical Journal, Vol.40 No.1, pp.43-63, 1961. [2] James W. Cooley, John W. Tukey. An algorithm for the machine calculation of complex Fourier series. Math. Comput., Vol.19, pp.297โ301, 1965. [3] S.M. Kay, S.L. Marple. Spectrum analysis - a modern perspective. Proc. IEEE, Vol.69, No.11, 1981. r.sc.comp. Vilnis Liepiลลก Email : [email protected] The extended summary of the doctoral thesis in Computer Science [4] Vilnis Liepins. A method for spectral analysis of band-limited signals, Automatic Control and Computer Sciences. Vol.27, No.5, pp.56-63, 1993. Available on researchgate.net. [5] Vilnis Liepins. A method of spectrum evaluation applicable to analysis of periodically and non-regularly digitized signals. Automatic Control and Computer Sciences, Vol.27, No.6, pp.57-64, 1993. Available on researchgate.net. [6] Vilnis Liepins. A spectral estimation method of nonuniformly sampled band-limited signals. Automatic Control and Computer Sciences, Vol.28, No.2, pp.66-73, 1994. Available on researchgate.net. [7] Vilnis Liepins. An algorithm for evaluation a discrete Fourier transform for incomplete data. Automatic Control and Computer Sciences, Vol.30, No.3, pp.27-40, 1996. Available on researchgate.net. [8] Vilnis Liepins. High-resolution spectral analysis by using basis function adaptation approach. Doctoral Thesis for Scientific Degree of Dr.Sc.Comp. / in Latvian /, University of Latvia, 1997. Available on researchgate.net. [9] M.D. Sacchi, T.J. Ulrych, C. Walker. Interpolation and extrapolation using a high-resolution discrete Fourier transform. IEEE Trans. on Signal Processing, Vol.46, No.1, pp.31-38, 1998. [10] Modris Greitans. Multiband signal processing by using nonuniform sampling and iterative updating of autocorrelation matrix. Proceedings of the 2001 International Conference on Sampling Theory and Application, May 13-17, 2001, Orlando, Florida, USA, pp.85-89. [11] Modris Greitans. Spectral analysis based on signal dependent transformation. The 2005 International Workshop on Spectral Methods and Multirate Signal Processing, (SMMSP 2005), June 20-22, 2005, Riga, Latvia. [12] Jayme Garcia Arnal Barbedo, Amauri Lopes, Patrick J. Wolfe. High Time-Resolution Estimation of Multiple Fundamental Frequencies. Proceedings of the 8th International Conference on Music Information Retrieval, ISMIR 2007, Sept.23-27, Vienna, Austria, pp.399-402. [13] David Dolenc, Barbara Romanowicz, Paul McGill, William Wilcock. Observations of infragravity waves at the ocean-bottom broadband seismic stations Endeavour (KEBB) and Explorer (KXBB). Geochemistry, Geophysics, Geosystems, Vol.9, Issue 5, May 2008. [14] Quanming Zhang, Huijin Liu, Hongkun Chen, Qionglin Li, Zhenhuan Zhang. A Precise and Adaptive Algorithm for Interharmonics Measurement Based on Iterative DFT. IEEE Trans on Power Delivery, Vol.23, Issue 4, pp.1728-1735, Oct.2008. [15] Petre Stoica, Jian Li, Hao He. Spectral Analysis of nonuniformly Sampled Data: A New Approach Versus the Periodogram. IEEE Trans on Signal Processing, Vol.57, Issue 3, pp.843-858, March 2009. [16] Eric Greenwood, Fredric H. Schmitz. Separation of Main and Tail Rotor Noise Sources from Ground-Based Acoustic Measurements Using Time-Domain De-Dopplerization. 35th European Rotorcraft Forum 2009, Sept.22-25, Hamburg, Germany. [17] Jayme Garcia Arnal Barbedo, Amauri Lopes, Patrick J. Wolfe. Empirical Methods to Determine the Number of Sources in Single-Channel Musical Signals. IEEE Transactions on Audio, Speech and Language Processing, Vol.17, Issue 7, pp.1435-1444, Sept.2009. [18] Tarik Yardibi, Jean Li, Petre Stoica, Ming Xue, Arthur B. Baggeroer. Source localization and sensing: A nonparametric iterative adaptive approach based on weighted least squares. IEEE Transactions on Aerospace and Electronic Systems, Vol.46, pp.425-443, Jan.2010. [19] M. Caciotta, S. Giarnetti, F. Leccese, Z. Leonowicz. Comparison between DFT, adaptive window DFT and EDFT for power quality frequency spectrum analysis. Modern Electric Power Systems (MEPS), 2010 Proceedings of the International Symposium, Sept.20-22, 2010, Wroclaw, pp.1-5. r.sc.comp. Vilnis Liepiลลก Email : [email protected] The extended summary of the doctoral thesis in Computer Science [20] Li Li, Yan Zheng, Wang Xing-zhi. Inter-harmonic Analysis Using IGG and Extended Fourier. Proceedings of the Chinese Society of Universities for Electric Power System and its Automation, 22(3), 2010. [21] Erik Gudmundson, Andreas Jakobsson, Jรถrgen Jensen, Peter Stoica. An Iterative Adaptive Approach for Blood Velocity Estimation Using Ultrasound. EUSIPCO 2010, Aug 23-27, Aalborg, Denmark, pp.348-352. [22] Modris Greitans, Rolands Shavelis. Reconstruction of sequences of arbitrary shaped pulses from its low pass or band pass approximations using spectrum extrapolation. EUSIPCO 2010, Aug. 23-27, Aalborg, Denmark, pp.1607-1611. [23] Juggrapong Treetrong. Fault Detection of Electric Motors Based on Frequency and Time-Frequency Analysis using Extended DFT. International Journal of Control and Automation, Vol.4, No.1, March 2011. [24] Jesper Rindom Jensen, Mads Grรฆsbรธll Christensen, Sรธren Holdt Jensen. A Single Snapshot Optimal Filtering Method for Fundamental Frequency Estimation. 36th International Conference on Acoustics, Speech and Signal Processing (ICASSP), Prague, Czech Republic, May 22-27, 2011, pp. 4272-4275. [25] Ming Xue, Luzhou Xu, Jian Li. IAA spectral estimation: fast implementation using the Gohberg-Semencul factorization. ICASSP 2011, May 22-27, Prague, Czech Republic, pp.3251-3261. [26] Bonifatius Wilhelmus Tilma. Supervisor: M.K. Smit; Co-promotor: E.A.J.M. Bente. Integrated tunable quantum-dot laser for optical coherence tomography in the 1.7ยตm wavelength region. Eindhoven, Technische Universiteit Eindhoven, Diss., Jun.2011. [27] Elmar Mair, Michael Fleps, Michael Suppa, and Darius Burschka. Spatio-temporal initialization for IMU to camera registration. In Proceedings of the IEEE International Conference on Robotics and Biomimetics (ROBIO), Dec.2011, pp.557-564. [28] George-Othan Glentis, Andreas Jakobsson. Superfast Approximative Implementation of the IAA Spectral Estimate. IEEE Transactions on Signal Processing, Vol.60, Issue 1, Jan.2012, pp.472-478. [29] B. W. Tilma, Yuqing Jiao, J. Kotani, B. Smalbrugge, H. P. M. M. Ambrosius, P. J. Thijs, X. J. M. Leijtens, R. Notzel, M. K. Smit, E. A. J. M. Bente. Integrated Tunable Quantum-Dot Laser for Optical Coherence Tomography in the 1.7 ๏ญ m Wavelength Region. IEEE Journal of Quantum Electronics, Vol.48, No.2, Feb.2012, pp. 87-98. [30] T. Odstrcil, M. Odstrcil, O. Grover, V. Svoboda, I. ฤuran, and J. Mlynรกล. Low cost alternative of high speed visible light camera for tokamak experiments. Review of Scientific Instruments. Oct.2012, Vol.83, Issue 10. [31] Elmar Mair. Co-promotor: Gregory Donald Hager. Efficient and Robust Pose Estimation Based on Inertial and Visual Sensing. Muฬnchen, Technische Universitaฬt Muฬnchen, Diss., 2012. [32] Akash K Singh. Quantum-Dot Laser OCT. International Journal of Engineering Research and Applications (IJERA), Vol.2, Issue 6, Nov.- Dec.2012, pp.340-371. [33] Elliot Briggs, OFDM Physical Layer Architecture and Real-Time Multi-Path Fading Channel Emulation for the 3GPP Long Term Evolution Downlink. PhD Thesis. Texas Tech University, Dec. 2012. [34] Kwadwo S. Agyepong. Fang-Han Hsu, Edward R. Dougherty and Erchin Serpedin. Spectral Analysis on Time-Course Expression Data: Detecting Periodic Genes Using a Real-Valued Iterative Adaptive Approach. Advances in Bioinformatics, Vol.2013. [35] Matthew James Kelley. Terahertz time domain spectroscopy of amino acids and sugars. Dissertation (Ph.D.), California Institute of Technology, March 2013. [36] Vsevolod Kharyton, Ronnie Bladh. Using Tiptiming and Strain Gauge Data for the Estimation of Consumed Life in a Compressor Blisk Subjected to Stall-Induced Loading. ASME r.sc.comp. Vilnis Liepiลลก Email : [email protected] The extended summary of the doctoral thesis in Computer Science Turbo Expo 2014: Turbine Technical Conference and Exposition, Vol.7B, June 16โ20, 2014, Duฬsseldorf, Germany. [37] Hui Wang, Cancan Liu, Bing Zhu, Juanjuan Cai & Yutian Wang. Multi-pitch estimation based on correlation and spectrum analysis. ICCT, Future Communication Technology, 2014. [38] D. Krause, W. B. Hussein, M. A. Hussein, T. Becker. Ultrasonic sensor for predicting sugar concentration using multivariate calibration. Ultrasonics, 54(6), 1703-1712, Aug.2014. [39] Petre Stoica, Dave Zachariah, Jian Li. Weighted SPICE: A Unifying Approach for Hyperparameter-Free Sparse Estimation, Digital Signal Processing, Vol.33, pp.1-12, Oct.2014. [ D.T. Michel, A.K. Davis, W. Armstrong, R. Bahr, R. Epstein, V.N. Goncharov, M. Hohenberger, I.V. Igumenshchev, R.Jungquist, D.D. Meyerhofer, P.B. Radha, T.C. Sangster, C. Sorce and D.H. Froula. Measurements of the ablationfront trajectory and low-mode nonuniformity in direct-drive implosions using x-ray self-emission shadowgraphy. High Power Laser Science and Engineering, Vol.3, E19, July 2015. [41] Jan Kober, Zdenek Prevorovsky, Milan Chlada. In situ calibration of acoustic emission transducers by time reversal method. Sensors and Actuators A: Physical, Vol.240, pp.50โ56, April 2016. [42] Stephanie Tsuei, Mark B. Milam. Trajectory generation for constrained differentially flat systems with time and frequency domain objectives. IEEE 55th Conference on Decision and Control (CDC), pp.4172-4177, December 2016. [43] L.Luzi, A.Stevens, H.Yang, M.D.Browning. Resolution Versus Error for Computational Electron Microscopy. Microscopy and Microanalysis, 23(S1), pp.88-89, August 2017. [44] Stavroula Karatasou and Mat Santamouris, Multifractal Analysis of High-Frequency Temperature Time Series in the Urban Environment, Climate 2018, 6(2), 50. [45][45]