DISCRETE FOURIER TRANSFORM IN A COMPLEX VECTOR SPACE

Abstract

An image-based phase retrieval technique has been developed that can be used on board a space based iterative transformation system. Image-based wavefront sensing is computationally demanding due to the floating-point nature of the process. The discrete Fourier transform (DFT) calculation is presented in “diagonal” form. By diagonal we mean that a transformation of basis is introduced by an application of the similarity transform of linear algebra. The current method exploits the diagonal structure of the DFT in a special way, particularly when parts of the calculation do not have to be repeated at each iteration to converge to an acceptable solution in order to focus an image.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

This application relates to Discrete Fourier Transforms, and in particular, to Discrete Fourier Transforms in phase retrieval for use in waveform analysis.

2. Background

Phase retrieval is an image-based wavefront sensing method that utilizes point-source images (or other known objects) to recover optical phase information. The most famous application of phase retrieval was the diagnosis of the Hubble Space Telescope mirror edge defect discovered soon after the launch of Hubble and subsequent correction using Corrective Optics Space Telescope Axial Replacement (COSTAR.) It may be ironic that a phase retrieval wavefront-sensing method lies at the very heart of the commissioning process for Hubble's successor. The earliest suggestion of using phase retrieval as a wavefront sensing method for JWST was in 1989, nearly a year before the Hubble launch and deployment. Details documenting the Hubble Space Telescope phase retrieval analysis are discussed in the literature.

A number of image-based phase retrieval techniques have been developed that can be classified into two general categories, the (a) iterative-transform and (b) parametric methods. Modifications to the original iterative-transform approach have been based on the introduction of a defocus diversity function or on the input-output method. Various implementations of the focus-diverse iterative-transform method have appeared which deviate slightly by utilizing a single wavelength or by varying the placement and number of defocused image planes. Modifications to the parametric approach include minimizing alternative merit functions as well as implementing a variety of nonlinear optimization methods such as Levenburg-Marquardt, simplex, and quasi-Newton techniques. The concept behind an optical diversity function is to modulate a point source image in a controlled way; in principle, any known aberration can serve as a diversity function, but defocus is often the simplest to implement and exhibits no angular dependence as a function of the pupil coordinates

The discrete Fourier transform (DFT) calculation is presented in “diagonal” form. By diagonal we mean that a transformation of basis is introduced by an application of the similarity transform of linear algebra.

What is needed is to find a more efficient implementation of the DFT for applications using iterative transform methods, particularly phase retrieval.

BRIEF SUMMARY

The image-based wavefront sensing method of and apparatus for focusing by determining phase difference between received images, removing that phase difference by calculating a DFT for each signal and iteratively multiplying the resultant until a convergence is achieved comprises the steps of accepting point-source stellar objects to recover embedded optical phase information, receiving data from a sensor array wherein the data includes embedded phase information, calculating a similarity transform and pre-transforming the data once at the beginning of the calculation.

The image-based wavefront sensing method of and apparatus for further includes iteratively processing the data using arbitrary frequency spacing in a diagonal fashion until a wavefront convergence is achieved to determine a phase difference by implementing the DFT in a 1 dimensional linear array using just N operations, where 0<N<N×log(N)<∞, wherein N×log(N) is the number of operations of an FFT counterpart, back-transforming the data at the end of the calculation in a single step and focusing the received images based on the determined phase.

The general problem of diagonalizing the DFT matrix (DFT eigen-problem) has been around for years. Apparently some aspects of the problem are related to topics in number theory that were even solved by Gauss. Such considerations have not led to more efficient implementations of the DFT over the FFT implementation, but we note that for certain applications, e.g., phase retrieval, that it is possible to exploit the diagonal structure of the DFT in a special way, particularly when parts of the calculation do not have to be repeated at each iteration to converge to an acceptable solution.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

FIG. 1 is a functional block diagram of a system which contains a discrete Fourier transform processor in a complex vector space, according to an example embodiment;

FIG. 2 is a flowchart representing the methodology of a discrete Fourier transform processor in a complex vector space , according to an example embodiment; and

FIG. 3 is another flowchart representing the methodology of a discrete Fourier transform processor in a complex vector space, according to an example embodiment.

DETAILED DESCRIPTION

Detailed illustrative embodiments are disclosed herein. However, specific structural and functional details disclosed herein are merely representative for purposes of describing example embodiments. Example embodiments may, however, be embodied in many alternate forms and should not be construed as limited to only the embodiments set forth herein.

Accordingly, while example embodiments are capable of various modifications and alternative forms, embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that there is no intent to limit example embodiments to the particular forms disclosed, but to the contrary, example embodiments are to cover all modifications, equivalents, and alternatives falling within the scope of example embodiments. Like numbers refer to like elements throughout the description of the figures.

It will be understood that, although the terms first, second, etc. may be used herein to describe various elements, these elements should not be limited by these terms. These terms are only used to distinguish one element from another. For example, a first element could be termed a second element, and, similarly, a second element could be termed a first element, without departing from the scope of example embodiments. As used herein, the term “and/or” includes any and all combinations of one or more of the associated listed items.

The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of example embodiments. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises”, “comprising,”, “includes” and/or “including”, when used herein, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof.

Hereinafter, example embodiments will be described with reference to the attached drawings.

Example embodiments of the present invention include a Discrete Fourier Transform (DFT) in phase retrieval for use in waveform analysis in a space based application.

One such space based application under consideration to use DFT phase retrieval methods includes the James Webb Space Telescope (JWST), one of NASA's great observatories of the Origins program and is scheduled to launch later this decade. The JWST design includes using advanced optical technology and 0-30 grade beryllium mirrors. The primary mirror is a 6.5 meter segmented primary using 18 segments in a hexagonal array. The optical design and manufacturing tolerances are specified for diffraction-limited performance at 2 μm, and the telescope is designed for operation over a temperature range of 30-60 K with an areal density of less than 15 kg/m². As a result, JWST commissioning and periodic optical maintenance incorporate san image-based wavefront sensing and control (WFSC) system to align the mirror segments, minimize the effects of figure error, and position the secondary mirror of the three-mirror anastigmat design. The wavefront sensing method specified for JWST is image-based in the sense that point-source stellar images are collected and processed to recover optical phase information. The primary camera for this function is the JWST Near Infrared Camera (NIRCam), although it has been shown that the Near Infrared Spectrograph (NIRSpec) is suitable for this function as well.

FIG. 1 shows a block diagram of a satellite or platform in which an embodiment of the instant invention is contained. The sensor array 130 is connected both physically and electronically to the platform housing 120 which is in turn likewise connected the processor 110 which contains the primary embodiments of the instant invention. The platform housing includes all the necessary parts (not shown) of a satellite such as fuel, electrical power, communications etc. Image-based wavefront sensing is computationally demanding due to the floating-point nature of the process. With certain data configurations and large array sizes, and particularly for under-sampled data sets, the calculations can take tens of minutes or even hours to reach full convergence for the optical wavefront. In addition, processors are generally susceptible to opto-mechanical instabilities induced by environmental factors such as jitter and turbulence. These factors not only limit fine-phasing performance but also the stability of the phasing results. In such cases, an obvious mitigation strategy is to increase the computational bandwidth of the WFSC system, i.e., decrease the time required for an entire closed-loop WFSC cycle.

To improve convergence times without sacrificing accuracy and to mitigate these environmental factors, we have developed special-purpose floating-point intensive computing hardware. These processors communicate using a series of multi-processor computing architectures specifically tailored to the JWST fine-phasing methodology. The computing system is based on a cluster of 64 DSPs (digital signal processors) operating in parallel under a shared memory architecture. The system is portable and callable from a laptop computer using a typical Ethernet connection, reaching a sustained floating-point performance of approximately 100 Giga-Flops with 64 DSPs for development, testing and hardware redundancy.

An embodiment of the instant invention is illustrated by some practical aspects of calculating in a diagonal DFT basis. Aside from the application to iterative Fourier methods, calculating in the diagonal DFT basis also has an interesting geometrical interpretation as solving for the Fourier transform pair, (F_n, f_n), that are co-linear or “parallel” with one another. Additional discussion on this interpretation is given below.

As an overview we begin by writing down the matrix form of the DFT. Next the hybrid-index notation of Schouten is introduced as an efficient notation for performing the DFT calculations in a complex vector space. The similarity transformation is presented in detail for a low-dimensional DFT example (4×4 array size) and then applied to illustrate how the DFT can be carried out in a diagonal basis using the hybrid-index notation. The main results show that by “pre-transforming” the data once at the beginning of the calculation, and then “back-transforming” the data at the end of the calculation, the DFT can potentially be implemented as a 1-d vector product using just N operations, compared to the N×log(N) operations needed by the FFT counterpart. The “pre-” and “back-” transforms are implemented using the DFT similarity transform and we further note that the calculation can be implemented using arbitrary frequency spacing.

DFT as a Unitary Transform

The DFT calculation can by implemented in matrix form:

$\begin{array}{cc}F\left(v_m\right)\equiv \mathrm{DFT}\left(v_m\right)=\Delta x\sum _{n=0}^{N-1}f\left(x_n\right)V_{\mathrm{mn}},& \left(1\right)\end{array}$

where V_mn is a non-singular M×N Vandermonde array:

V _mn =e ^{−i2πv m (nΔx)} =e ^{−i2π(m·n)/M}=(e^−i2π/M)^m·n=ω_M^m·n,   (2)

with components given by

$\begin{array}{cc}{V}_{\mathrm{mn}}=\left(\begin{array}{ccccc}1& 1& 1& L& 1\\ 1& {\omega }_M& {\omega }_M^2& L& {\omega }_M^{N-1}\\ 1& {\omega }_M^2& {\omega }_M^{2·2}& L& {\omega }_M^{2·\left(N-1\right)}\\ M& M& M& O& M\\ 1& {\omega }_M^{M-1}& {\omega }_M^{\left(M-1\right)·2}& L& {\omega }_M^{\left(M-1\right)\left(N-1\right)}\end{array}\right).& \left(3\right)\end{array}$

Next we define W_mn as

W _mn =V _mn /√{square root over (M)}=(e^−i2π/M)^m·n/√{square root over (M)},   (4)

and then it is straightforward to show that W_mn is a unitary matrix:

W·W ^† =I,   (5)

using boldface matrix notation where I is the identity matrix and the “†” symbol denotes the complex-conjugate transpose or “adjoint.” We comment that the adjoint operation will be presented more generally in the next Section. Next letting

f _n ≡Δxf(x_n),   (6)

Equation (1) takes the form of a unitary transformation:

$\begin{array}{cc}F\left(v_m\right)\equiv F_m=\sum _{n=0}^{N-1}W_{\mathrm{mn}}f_n,& \left(7\right)\end{array}$

with the inverse transformation given by (generally but not always we take M=N):

$\begin{array}{cc}{f}_n=\sum _{m=0}^{M-1}{\left(W_{\mathrm{mn}}\right)}^{†}F_m.& \left(8\right)\end{array}$

Given (5) and (7), the DFT can be analyzed using the transformation theory of complex vector spaces. The hybrid-index notation of Schouten supplies a natural framework for this approach and is presented in the next Section.

Hybrid Index Notation

As illustrated above, the DFT is a unitary transformation in a complex vector space. Therefore, the calculations can be performed elegantly and efficiently using the tensor-transformation formalism of Schouten. In particular, Schouten's hybrid-index notation provides a natural framework for the DFT calculations. The goal in using this approach is to give a more general representation of the calculations by identifying components with quantities of a geometric manifold. Often, but not always, such identification can lead to alternative or more efficient calculation strategies as discussed later. Therefore, the goal in the next Section is to sort out the details of the hybrid-index notation for the DFT, and thus identify the DFT with quantities in a geometric manifold. In addition, by adopting a notation that makes the transformation rules of a quantity explicit, the meaning of the calculations can often be clarified and the notation serves as an efficient bookkeeping device for checking consistency throughout the calculations.

DFT in Hybrid-Index Notation

The DFT is defined here as a unitary transform between conjugate Fourier domains. The 1-d DFT transform is then expressed using Schouten's hybrid-index notation as

F ^m =W ^m _·n f ⁿ,   (9)

using the Einstein sum convention defined over repeated indices in an upper/lower combination:

$\begin{array}{cc}\sum _n^{}W_{\mathrm{mn}}f_n\equiv W_{·n}^mf^n.& \left(10\right)\end{array}$

The lower/upper indices denote the covariant and contravariant transformation rules—the upper indices are not exponents in this notation. Exponents of ( a quantity will be denoted using parentheses, or brackets, respectively, for either a transformation of coordinates, or a “point transformation.” Allowable transformations of coordinates for this application will be discussed later. The distinction between the transformations of coordinates and point transformations are fundamental to Schouten's kernel-index method. To elaborate briefly, coordinate transformations are denoted by a transformation of basis, say, m→m′, where the new basis set is denoted by m′. The components in the new basis are thus

u ^m ′=A _m ^m ′u _m,   (11)

where A_m^M′ is the transformation matrix and the u^m′ are the vector components in the new basis set (u is arbitrary).

A point transformation, on the other hand, is denoted by

U ^m =P _·n ^M u ⁿ,   (12)

and is the form adopted for the DFT in (9). Another important distinction between a transformations of basis and a point transformation is that the indices for the point transformation are ordered and the positions are marked using over- or under-dots as in P_·n^m (12). The main take-away is that for coordinate transformations, the kernel symbol, u, remains the same as in (11). But for point transformations, the kernel symbol changes (as in (12)) but the indices remain the same. So initially at least in this context, the DFT transformation in (9) is defined as a point transformation, but later we find that when calculating the DFT in a diagonal basis, it is natural to reconsider (9) using a common kernel symbol for the Fourier transform pair, specifically, (F^m, fⁿ)→(f_m, fⁿ). We note also that it is possible to mix the two transformations (basis and point) as is required when considering the allowable transformations of basis in digital signal processing. These details will be discussed in more detail later.

Continuing with an explanation of the formalism, the complex conjugate of (9) is given by

F ^m = W ^m _n f ⁿ ≡ f ^m,  (13)

using the over-bar to denote the complex conjugation of components and kernel symbols in accordance with Schouten's hybrid-index notation and noting for consistency that

F ^m =F^mF^m==F^m.   (14)

The conjugation of indices and kernel objects is fundamental to the study of transformations in a complex vector space, and as illustrated below, is perfectly suited for calculations of the DFT transform.

We choose an initial convention that the data values are labeled by contravariant components, f^m, and then by complex conjugation of the kernel symbol and indices in various combinations, several other types of components are represented in Schouten's notation by the “raising” and “lowering” of indices using the “fundamental” or metric tensor, α_kl:

f^m,f_m, f^m, f_m,   (15)

where the latter two follow by complex conjugation of the first two. Specifically, we mean that

f ^m f _m =a _mm f ^m ; f ^m f _m =a _mm f ^m ,   (16)

As mentioned above, co- and contravariant transformation rules can be switched back and forth using the metric or fundamental tensor, a_kl, to raise or lower indices as appropriate. The metric tensor also defines how the differential element of distance is calculated in the manifold as well as the inner product between two vectors, say, u^k and v^l:

u|v =a_klū^kv^l=ū_lv^t=ū^kv_k,   (17)

pointing out that in a linear complex vector space, the “bra-ket” notation of Dirac is easily incorporated in this approach as illustrated on the LHS of (17). Illustrating the use of the metric for raising and lowering indices we have:

W _·n ^m =a _ln W ^{m l} and fⁿ=ā^{n k}f_k,   (18)

and then substituting these expressions back into (9) we see that an alternative form of the DFT is obtained:

F ^m =W _·n ^m f ⁿ=(a_lnW^{m l})(a^{n k}f_k)=a_lna^{n k}W^{m l}f_k=δ_l^{· k}W^{m l}f_kF^m=W^{m n}f_n.   (19)

Adjoint in a Geometric Manifold

In general, the adjoint operation applied to the W_·n^m, denoted by {tilde over (W)}_·n^m=(W_·n^m)^†, is defined through a two-step operation using both the metric tensor and complex conjugation:

• • (a) each index is raised or lowered as appropriate using the metric:

a _nk a ^{l m} W _·l ^k =W _n ^{· m},   (20)

• • (b) the result in (a) is then complex conjugated to give

(W_·n^m)^†===W_n^·m, tm (21)

and finally we write

{tilde over (W)} _·n ^m=≡(W_·n^m)^†= W_n^·m,   (22)

For comparison, in matrix boldface notation:

(W)^\= W^T=a Wa⁻¹=W^{31 1},   (23)

where “T” defines the usual transpose operation. We note especially that the example in (23) illustrates a degree of non-generality in the way the adjoint is usually calculated in matrix notation. Specifically, defining the adjoint using components in a geometric manifold (using steps (a) and (b) above) can give different results than simply complex conjugating and then performing the transpose operation using the “T” as in (23). In special coordinates, a_nk≐δ_nk, these two are identical, but in general, that is not always the case for arbitrary a_nk.

Since the W_·n^m are unitary, the inverse,

${\stackrel{-1}{W}}_{·n}^m,$

is defined using the adjoint and is expressed (from (22)):

$\begin{array}{cc}{\tilde{W}}_{·n}^m={\stackrel{-1}{W}}_{·n}^m={\left(W_{·n}^m\right)}^{†}={\stackrel{\_}{W}}_n^{·m},& \left(24\right)\end{array}$

such that

W _·k ⁿ =δ_·k^l.  (25)

We further comment that Schouten did not initially use an over-bar on the kernel symbol as in we have in (25). Considering the W_·n^m using two upper indices, W^{m n} (see (18) and (19)) the corresponding adjoint is similarly calculated:

(W^{m n})^\=a_mka_lnW^{k l}==W_{m n}=.  (26)

The inverse DFT transformation to (9) is thus defined:

{tilde over (W)} _·m ⁿ F ^m=W_·k^mf^kδ_kⁿf^k=fⁿ,   (27)

as expected, or when using the form of (19) the inverse DFT is expressed

{tilde over (W)} _{m n} F ^m=(W^{m k}f_k)=δ_n^kf_k=f_n.   (28)

A few additional comments from this Section: in tensor transformation theory, there is geometric meaning attributed to these differing but complementary transformation rules. Specific transformations are defined by an allowable Group of transformations, G_a, as discussed further in Schouten. The corresponding allowable transformations of basis for the DFT application (digital signal processing) can include fractional and integer shifts of the data origin in addition to scaling and phase transformations. The details for these various transformations in this notation will be discussed in a later set of notes.

Matrix Multiplication Using Index Notation

When calculating with components, index order is important for proper matrix multiplication and thus we adopt the following convention for a rank-2 quantity:

$\begin{array}{cc}{A}_{·\underset{\underset{\mathrm{col}}{\downarrow }}{n}}^{m\to \mathrm{row}},& \left(29\right)\end{array}$

and the inverse of A_·n^m is notated as,

${\stackrel{-1}{A}}_{·n}^m,$

such that

$A_{·n}^m{\stackrel{-1}{A}}_{·k}^n={\delta }_{·k}^m$

where δ_·k^m is the Kronecker delta. To illustrate, the matrix multiplication of two matrices A and B is notated as:

C _·k ^m =A _·n ^m B _·k ⁿ, (30)

and to be specific we calculate (30) using Mathematica with 2×2 arrays:

Table [Sum [(A[[um,n]]B[[n,lk]]), {n,1,NN}], {um,1,2}, {lk,1, 2}, ]]  (31)

where um denotes “upper-m” and lk denotes “lower-k;” the upper/lower dummy variable, n, is left un-notated in Mathematica with regard to its position. The result is of course (30):

$\begin{array}{cc}{C}_{·k}^m=A_{·n}^mB_{·k}^n=\mathrm{AB}=\left(\begin{array}{cc}{a}_{·1}^1b_{·1}^1+a_{·2}^1b_{·1}^2& a_{·1}^1b_{·2}^1+a_{·2}^1b_{·2}^2\\ a_{·1}^2b_{·1}^1+a_{·2}^2b_{·1}^2& a_{·1}^2b_{·2}^1+a_{·2}^2b_{·2}^2\end{array}\right),& \left(32\right)\end{array}$

where AB is the conventional matrix product using matrix boldface notation. For comparison, transposing both the m-n and n-k index positions in (30) reverses the order and gives the BA matrix product:

$\begin{array}{cc}{C}_{·k}^m=A_n^{·m}B_k^{·n}=\mathrm{BA}=\left(\begin{array}{cc}{b}_{·1}^1a_{·1}^1+b_{·2}^1a_{·1}^2& b_{·1}^1a_{·2}^1+b_{·2}^1a_{·2}^2\\ b_{·1}^2a_{·1}^1+b_{·2}^2a_{·1}^2& b_{·1}^2a_{·2}^1+b_{·2}^2a_{·2}^2\end{array}\right),& \left(33\right)\end{array}$

and the corresponding Mathematica code for (33) is

Table[Sum[(A[[n,um]]B[[lk,n]]), {n,1,NN}], {um,1,2}, {lk,1,2}]  (34)

Metric Tensor Deriving from Unitary Invariance

In summary of Schouten's notation, various geometric quantities are defined that are distinguished by their transformation rules using index position. A group of allowable transformations, G_a, is introduced and then several invariant quantities may be constructed based on invariance of a given symmetry group. One such quantity is the metric tensor, a_kl, defining “distance” in the manifold and also defining the inner product between any two vectors. Therefore, one goal in these notes is to solve for the DFT metric tensor, a_kl, by postulating nothing more than unitary invariance under the action of the point transformation, W_·n^m, as illustrated below. Specifically, given that the W_·n^m are defined as a unitary transformation, we introduce a metric tensor defined by invariance under the action of the W_·n^m transformation such that

W _·l ^m W _k ^{· n} a _nm =a _kl,   (35)

or in matrix form:

W a W ⁻¹ =a.  (36)

Quantities that are invariant under the action of the unitary group are hermitian and thus the a_kl satisfy

a _kl=(a kl)^\=ā_lk,   (37)

The introduction of a hermitian metric a_kl into the E_n geometric manifold defines the manifold as a U_n which differs from the earlier notation used by Schouten where this manifold was labeled as the {tilde over (R)}_n.

General form of Parseval's Theorem

We can show that Parseval's theorem is really implied by unitary invariance, and not necessarily “energy conservation” as it is commonly referred to in the engineering literature. To begin, we consider the magnitude squared of F^m→|F|², giving

F|F =|F|²=a_mnF^mFⁿ.  (38)

Substituting the DFT transformation rule from (9) into (38) we have

|F|²=a_nmW_·k^mW_l^{· n}f^kf^l,   (39)

but since the a_kl are invariant under this unitary transformation (see (35)) then the RHS of (39) simplifies and we get

|F|²=a_klf^kf^l.   (40)

Therefore, by comparison with (38), Parseval's theorem is seen to hold for any solution of the a_klas implied by unitary invariance (35), and thus in general:

a _mn F ⁿ F ⁿ =a _kl f ^k f ^l.   (41)

The common form of Parseval's theorem is expressed in vector notation as

|F|²=|f|².  (42)

To give some perspective on these results we look forward a bit and comment that the standard DFT approach adopts a special solution (coordinate system) for the a_kl given by the Kronecker delta (as shown later):

a _nk≐δ_nk  (43)

and so in index notation (42) takes the familiar “Euclidean” form:

δ_mnF^mFⁿ=δ_klf^kf^l,   (44)

or as it is commonly written:

|F ₁|²+|F₂|²+L+|F_N|²=|f₁|²+|f₂|²+L+|f_N|².  (45)

The point is that for a general a_kl, Parseval's theorem (41) will not necessarily take the Euclidean form as demonstrated later below (see (67)), although unitary invariance implies that (41) still holds.

A couple of additional comments: as discussed later the δ_kl is a solution for the a_kl and corresponds to a positive definite signature solution. The approach taken here is motivated by that fact that in general, the eigen-problem requires the use of a general metric, and by no means does this necessarily correspond to the Kronecker delta in a unitary complex vector space, U_n. But even so, we later adopt the Kronecker delta solution to make progress on the main point of calculating the DFT in a diagonal basis, and for easy comparison with previously published results.

Metric Solutions

By postulating unitary invariance of the metric we can solve for the a_kl using (35) and (37). The W_·k^m are known through identification with (4):

$\begin{array}{cc}{W}_{·k}^m\equiv \frac{1}{\sqrt{M}}\left({}^{-\mathrm{2\pi }\left(\mathrm{mk}\right)/M}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right).& \left(46\right)\end{array}$

As discussed above, the inverse transformation of (46) is defined through the adjoint operation in (20). But at this point in the discussion we have not yet solved for the metric a_kl, which is needed to form the adjoint. Therefore, we first solve for the

${\stackrel{-1}{W}}_{·k}^n$

such that

$\begin{array}{cc}{W}_{·n}^m{\stackrel{-1}{W}}_{·k}^n={\delta }_{·k}^m,& \left(47\right)\end{array}$

and the

${\stackrel{-1}{W}}_{·k}^n$

are given by:

$\begin{array}{cc}{\stackrel{-1}{W}}_{·k}^m={\stackrel{\_}{W}}_k^{·m}=\frac{1}{\sqrt{M}}\left({}^{\mathrm{2\pi }\left(\mathrm{mk}\right)/M}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right).& \left(48\right)\end{array}$

Continuing to the general ×2 metric, a_kl defined as

$\begin{array}{cc}{a}_{\stackrel{\_}{k}l}=\left(\begin{array}{cc}{z}_{\stackrel{\_}{1}1}& z_{\stackrel{\_}{1}2}\\ z_{\stackrel{\_}{2}1}& z_{\stackrel{\_}{2}2}\end{array}\right)=\left(\begin{array}{cc}{x}_{11}+y_{11}& x_{12}+y_{12}\\ x_{21}+y_{21}& x_{22}+y_{22}\end{array}\right),& \left(49\right)\end{array}$

where the complex notation is dropped on indices of real components. The condition for unitary invariance (35) then gives the following simultaneous system of 4 equations and 8 unknowns:

$\begin{array}{cc}\frac{1}{2}\left(\begin{array}{cc}{z}_{\stackrel{\_}{1}1}+z_{\stackrel{\_}{1}2}+z_{\stackrel{\_}{2}1}+z_{\stackrel{\_}{2}2}& z_{\stackrel{\_}{1}1}-z_{\stackrel{\_}{1}2}+z_{\stackrel{\_}{2}1}-z_{\stackrel{\_}{2}2}\\ z_{\stackrel{\_}{1}1}+z_{\stackrel{\_}{1}2}-z_{\stackrel{\_}{2}1}-z_{\stackrel{\_}{2}2}& z_{\stackrel{\_}{1}1}-z_{\stackrel{\_}{1}2}-z_{\stackrel{\_}{2}1}+z_{\stackrel{\_}{2}2}\end{array}\right)=\left(\begin{array}{cc}{z}_{\stackrel{\_}{1}1}& z_{\stackrel{\_}{1}2}\\ z_{\stackrel{\_}{2}1}& z_{\stackrel{\_}{2}2}\end{array}\right).& \left(50\right)\end{array}$

The system (50) is under-determined and can therefore have many possible solutions but will be further constrained by the hermitian condition (37). Solving (50) eliminates 4 of the unknowns (z₁₁ and z₁₂) in terms of the other 4 unknowns (z₂₁ and z₂₂) showing that a_kl has the general form:

$\begin{array}{cc}{a}_{\stackrel{\_}{k}l}=\left(\begin{array}{cc}2z_{\stackrel{\_}{2}1}+z_{\stackrel{\_}{2}2}& z_{\stackrel{\_}{2}1}\\ z_{\stackrel{\_}{2}1}& z_{\stackrel{\_}{2}2}\end{array}\right).& \left(51\right)\end{array}$

Next applying the hermitian (symmetry) condition (35) to (51) gives the condition that the 2×2 a_kl must be real valued:

$\begin{array}{cc}\left(\begin{array}{cc}2y_{21}+y_{22}& y_{21}\\ y_{21}& y_{22}\end{array}\right)=\left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right)\Rightarrow \mathrm{Im}\left(a_{\stackrel{\_}{k}l}\right)=0,& \left(52\right)\end{array}$

but as expected there are multiple solutions for the a_kl in terms of the 2 final parameters x₂₁ and x₂₂. The general solution for the 2×2 case is thus reduced from 8 to 2 unspecified real components and is given by

$\begin{array}{cc}{a}_{\stackrel{\_}{k}l}=\left(\begin{array}{cc}2x_{21}+x_{22}& x_{21}\\ x_{21}& x_{22}\end{array}\right).& \left(53\right)\end{array}$

As a check, substituting (53) back into (35) verifies that this a_kl is invariant under the action of the W_·k^m for arbitrary (but real) numerical values of x₂₁ and x₂₂ in (53).

We have shown that an infinite number of solutions exist for the 2×2 metric, a_kl, as implied by unitary invariance in terms of two real parameters. We note that one of the simplest and most important solutions is a special case of (53) when x₂₁=0 and x₂₂=1 giving the Kronecker delta as mentioned earlier:

$\begin{array}{cc}{\underset{1}{a}}_{\stackrel{\_}{k}l}{}_{\underset{x_{22}\to 1}{x_{21}\to 0}}\stackrel{*}{=}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\stackrel{*}{=}{\delta }_{\stackrel{\_}{k}l},& \left(54\right)\end{array}$

where the complex conjugation is kept for consistency, and the numeral “1” beneath the kernel symbol in Schouten's notation labels this case as a distinct solution, different from other solutions, say

${\underset{2}{a}}_{\stackrel{\_}{k}l}.$

The labeling of kernel symbols in this way is used extensively in the eigen-problem discussed later.

We comment that since (54) is positive definite, as discussed at greater length below in the Section titled Generalized DFT, there will be no distinction between the other various quantities that may be derived from the a_kl in the U₂, and thus, no numerical distinction exists for components labeled by upper/lower indices, or even the complex conjugation of indices. This is seen in the solution (54) above.

As illustrated above there are an infinite number of solutions for the a_kl. Several specific solutions are shown in (55) based on various choices for the x₂₁ and x₂₂ components:

$\begin{array}{cc}{\underset{1}{a}}_{\stackrel{\_}{k}l}{}_{\underset{x_{22}\to 1}{x_{21}\to 0}}\stackrel{*}{=}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right);{\underset{2}{a}}_{\stackrel{\_}{k}l}{}_{\underset{x_{22}\to -1}{x_{21}\to 1}}\stackrel{*}{=}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right){\underset{3}{a}}_{\stackrel{\_}{k}l}{}_{\underset{x_{22}\to 0}{x_{21}\to 1}}\stackrel{*}{=}\left(\begin{array}{cc}2& 1\\ 1& 0\end{array}\right);{\underset{4}{a}}_{\stackrel{\_}{k}l}{}_{\underset{x_{22}\to 1}{x_{21}\to -\frac{1}{2}}}\stackrel{*}{=}\left(\begin{array}{cc}0& -\frac{1}{2}\\ -\frac{1}{2}& 1\end{array}\right),& \left(55\right)\end{array}$

Given some specific solutions for the a_kl in (55), we examine various other quantities that are derived in Schouten's notation. For instance, consider the second solution in (55) that is a multiple of the W_·n^m:

$\begin{array}{cc}{a}_{\stackrel{\_}{k}l}\equiv {\underset{2}{a}}_{\stackrel{\_}{k}l}\stackrel{*}{=}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)\Rightarrow {\stackrel{\_}{a}}^{k\stackrel{\_}{l}}\stackrel{*}{=}\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right),& \left(56\right)\end{array}$

and then for generality define a Complex data vector:

f ⁿ=(z¹,z²)=(x¹+iy¹,x²+iy²)  (57)

The DFT transformation (9) is thus given by

$\begin{array}{cc}{F}^m=W_{·n}^mf^n=\frac{1}{\sqrt{2}}\left(z^1+z^2,z^1-z^2\right).& \left(58\right)\end{array}$

As a consistency check on both the solution for the a_kl of (56), and the legitimacy of constructing the inverse transformation by use of the adjoint operation in (20), i.e.,

${\stackrel{-1}{W}}_{·n}^m={\left(W_{·n}^m\right)}^{†}\equiv \text{?},\text{?}\text{indicates text missing or illegible when filed}$

we can check the inverse transformation calculated by way of the adjoint:

$\begin{array}{cc}\begin{array}{c}{\left(W_{·n}^m\right)}^{†}=a_{\stackrel{\_}{n}k}a^{l\stackrel{\_}{m}}W_{·l}^k\\ =\text{?}\frac{1}{\sqrt{2}}\text{?}\frac{1}{2}\text{?}\\ =\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)\\ ={\stackrel{-1}{W}}_{·n}^m,\end{array}\text{?}\text{indicates text missing or illegible when filed}& \left(59\right)\end{array}$

in agreement with (48). The inverse transformation is also for consistency giving

$\begin{array}{cc}{\tilde{W}}_{·n}^mF^m=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)\left(\begin{array}{c}{z}^1+z^2\\ z^1-z^2\end{array}\right)=\left(z^1,z^2\right)=f^n.& \left(60\right)\end{array}$

Several other quantity are derived from the fⁿ using the a_kl:

f ⁿ=(z¹,z²); f_n=a_nkf^k=(z¹+z²,z¹−z²) = f^k=(z¹, z²); f_k=a_nkfⁿ=(z¹+ z², z¹− z²),  (61)

and comparing the RHS's in (61), as expected we have f_n=, showing consistency using the hybrid-index notation. Correspondingly, various other F^m quantities are given by

$\begin{array}{cc}{F}^m=\frac{1}{\sqrt{2}}\left(z^1+z^2,z^1-z^2\right);F_{\stackrel{\_}{m}}=a_{\stackrel{\_}{m}k}F^k=\sqrt{2}\left(z^1,z^2\right){\stackrel{\_}{F}}^{\stackrel{\_}{m}}=\frac{1}{\sqrt{2}}\left({\stackrel{\_}{z}}^{\stackrel{\_}{1}}+{\stackrel{\_}{z}}^{\stackrel{\_}{2}},{\stackrel{\_}{z}}^{\stackrel{\_}{1}}-{\stackrel{\_}{z}}^{\stackrel{\_}{2}}\right);{\stackrel{\_}{F}}_m=a_{\stackrel{\_}{n}m}{\stackrel{\_}{F}}^{\stackrel{\_}{n}}=\sqrt{2}\left({\stackrel{\_}{z}}^{\stackrel{\_}{1}},{\stackrel{\_}{z}}^{\stackrel{\_}{2}}\right)& \left(62\right)\end{array}$

and components for the W^{m n} and the W_mn are also calculated:

$\begin{array}{cc}{W}^{m\stackrel{\_}{n}}=a^{k\stackrel{\_}{n}}W_{·k}^m=\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right);W_{\stackrel{\_}{m}n}=a_{\stackrel{\_}{k}n}W_{\stackrel{\_}{m}}^{·\stackrel{\_}{k}}=\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right).& \left(63\right)\end{array}$

Using the f_n from (61) and then comparing with the F^m in (62) we can check the forward DFT for consistency using the W^{m n} of (63):

$\begin{array}{cc}\begin{array}{c}{W}^{m\stackrel{\_}{n}}f_{\stackrel{\_}{n}}=\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right)\left(\begin{array}{c}{z}^1+z^2\\ z^1-z^2\end{array}\right)\\ =\frac{1}{\sqrt{2}}\left(z^1+z^2,z^1-z^2\right)\\ =F^m.\end{array}& \left(64\right)\end{array}$

The adjoint operation on the W^{m n} and W_mn of (63) is similarly verified:

$\begin{array}{cc}\begin{array}{c}{\left(W_{·n}^m\right)}^{†}=\text{?}\\ =a_{m\stackrel{\_}{k}}a_{\stackrel{\_}{n}l}W^{m\stackrel{\_}{n}}\\ =\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)\left(\begin{array}{cc}\frac{1}{\sqrt{2}}& 0\\ 0& \frac{1}{\sqrt{2}}\end{array}\right)\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)\\ =\left(\begin{array}{cc}\sqrt{2}& 0\\ 0& \sqrt{2}\end{array}\right),\end{array}\text{?}\text{indicates text missing or illegible when filed}& \left(65\right)\end{array}$

such that

W ^{m n} =δ_kⁿ.   (66)

Finally, Parseval's theorem is demonstrated:

a _mn F ^m F ⁿ =a _kl f ^k f ^l =|z ¹|²−|z²|^s,   (67)

which is not in Euclidean form but is nevertheless invariant.

Comments on a General DFT

Further restrictions may be placed on the solution values for the x₂₁ and x₂₂ in (53), in accordance with the classification of the a_kl as being positive definite, negative definite, or indefinite, corresponding to the number of values (or signature)>0, <0, or mixed, respectively, along the diagonal of the a_kl). Note that when a_kl is definite that there is no distinction between the various quantities that may be derived from a_kl in the U₂ and thus, no numerical distinction exists for components labeled by upper/lower indices, or even the complex conjugation of indices; these notations are merely formal in this instance. We comment further that since the a_kl are required when solving the eigen-problem in general (more details on this below) there can in principle at least, be different solutions for the eigen-problem corresponding to different solutions for the a_kl. Note that all earlier work on the DFT eigen-problem uses the a_kl≐δ_kl solution without considering the possibility that the a_kl might be defined in a more general way based on unitary invariance. Thus, different solutions for the a_kl may lead to possibly different DFT eigen-quantities, but all the while preserving unitary invariance and the hermitian character of the a_kl. Thus, from the perspective of the transformation theory of complex vector spaces, we can regard some properties of the DFT as deriving from the characteristics of the a_kl, i.e., that (a) a_kl is invariant under unitary transformations and (b) is hermitian. Condition (b) really derives from (a) but it is convenient when solving for the a_kl components (as illustrated above) to regard (b) as an independent condition.

Similarly as in the 2×2 example, we can solve for the a_kl corresponding to higher dimensional data sets. Of special interest are solutions in which the a_kl are complex valued since these solutions provide further consistency checks on the proper use of Schouten's' notation. We proceed by deriving the 4×4 W_·k^m and its inverse,

${\stackrel{-1}{W}}_{·k}^m,$

from (46):

$\begin{array}{cc}{W}_{·k}^m=\frac{1}{2}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& -i& -1& i\\ 1& -1& 1& -1\\ 1& i& -1& -i\end{array}\right);\begin{array}{c}{\stackrel{-1}{W}}_{·k}^m=\text{?}\\ ={\stackrel{\_}{W}}_k^{·m}\\ =\frac{1}{2}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& i& -1& -i\\ 1& -1& 1& -1\\ 1& -i& -1& i\end{array}\right),\end{array}\text{?}\text{indicates text missing or illegible when filed}& \left(68\right)\end{array}$

and a general 4×4 metric is defined by

$\begin{array}{cc}{a}_{\stackrel{\_}{k}l}\left(\begin{array}{cccc}{z}_{\stackrel{\_}{1}1}& z_{\stackrel{\_}{1}2}& z_{\stackrel{\_}{1}3}& z_{\stackrel{\_}{1}4}\\ z_{\stackrel{\_}{2}1}& z_{\stackrel{\_}{2}2}& z_{\stackrel{\_}{2}3}& z_{\stackrel{\_}{2}4}\\ z_{\stackrel{\_}{3}1}& z_{\stackrel{\_}{3}2}& z_{\stackrel{\_}{3}3}& z_{\stackrel{\_}{3}4}\\ z_{\stackrel{\_}{4}1}& z_{\stackrel{\_}{4}2}& z_{\stackrel{\_}{4}3}& z_{\stackrel{\_}{4}4}\end{array}\right)& \left(69\right)\end{array}$

The condition for unitary invariance (35) then gives a simultaneous system of 16 equations and 32 unknowns which can be further reduced to 6 free parameters by imposing the hermitian condition (37) on the a_kl. After considerable algebra the general solution to (35) for the a_kl, constrained by the hermitian condition (37) is given by

$\begin{array}{cc}{a}_{\stackrel{\_}{k}l}=\left(\begin{array}{cccc}{x}_{33}+2\left(x_{41}+x_{43}\right)& x_{41}-{\mathrm{iy}}_{43}& \left(x_{33}+x_{41}-x_{42}+x_{43}-x_{44}+2{\mathrm{iy}}_{43}\right)& x_{41}-{\mathrm{iy}}_{43}\\ x_{41}+{\mathrm{iy}}_{43}& x_{44}& x_{43}+{\mathrm{iy}}_{43}& x_{42}\\ \left(x_{33}+x_{41}-x_{42}+x_{43}-x_{44}+2{\mathrm{iy}}_{43}\right)& x_{43}-{\mathrm{iy}}_{43}& x_{33}& x_{43}-{\mathrm{iy}}_{43}\\ x_{43}+{\mathrm{iy}}_{43}& x_{42}& x_{43}+{\mathrm{iy}}_{43}& x_{44}\end{array}\right).& \left(70\right)\end{array}$

As before, there are an infinite number of both real and complex valued 4×4 solutions for the a_kl in terms of the 6 free parameters in (70). Several example solutions corresponding to basic numerical values for the x and y components are listed in (71):

$\begin{array}{cc}{\underset{1}{a}}_{k\stackrel{\_}{l}}{}_{\underset{x_{43}\to 0,y_{43}\to 0}{\underset{x_{41}\to 0,x_{42}\to 0}{x_{33}\to 1,x_{44}\to 1}}}\stackrel{*}{=}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right);{\underset{2}{a}}_{k\stackrel{\_}{l}}{}_{\underset{x_{43}\to 0,y_{43}\to 1}{\underset{x_{41}\to 0,x_{42}\to 0}{x_{33}\to 1,x_{44}\to 1}}}\stackrel{*}{=}\left(\begin{array}{cccc}1& -i& 2i& -i\\ i& 1& i& 0\\ -2i& -i& 1& -i\\ i& 0& i& 1\end{array}\right){\underset{3}{a}}_{k\stackrel{\_}{l}}{}_{\underset{x_{43}\to 1,y_{43}\to 0}{\underset{x_{41,}\to -1x_{42}\to 0}{x_{33}\to 1,x_{44}\to 1}}}\stackrel{*}{=}\left(\begin{array}{cccc}1& -1& 0& -1\\ -1& 1& 1& 0\\ 0& 1& 1& 1\\ -1& 0& 1& 1\end{array}\right);{\underset{4}{a}}_{k\stackrel{\_}{l}}{}_{\underset{x_{43}\to 0,y_{43}\to 0}{\underset{x_{41}\to -1,x_{42}\to 0}{x_{33}\to 1,x_{44}\to 1}}}\stackrel{*}{=}\left(\begin{array}{cccc}-1& -1& -1& -1\\ -1& 1& 0& 0\\ -1& 0& 1& 0\\ -1& 0& 0& 1\end{array}\right).& \left(71\right)\end{array}$

Once again we pick a solution and then illustrate a few of the basic quantities that are derived for the DFT in Schouten's notation. For generality, consider the complex-valued solution of (71) for the a_kl:

$\begin{array}{cc}{a}_{\stackrel{\_}{k}l}\equiv {\underset{2}{a}}_{\stackrel{\_}{k}l}\stackrel{*}{=}\left(\begin{array}{cccc}1& -i& 2i& -i\\ i& 1& i& 0\\ -2i& -i& 1& -i\\ i& 0& i& 1\end{array}\right),& \left(72\right)\end{array}$

and then the inverse is given by

$\begin{array}{cc}{\stackrel{-1}{a}}_{\stackrel{\_}{k}l}=a^{k\stackrel{\_}{l}}\stackrel{*}{=}\frac{1}{7}\left(\begin{array}{cccc}1& -2-i& -2+2i& -2-i\\ -2+2i& 5& 2+i& -2\\ -2-2i& 2-i& 1& 2-i\\ -2+i& -2& 2+i& 5\end{array}\right).& \left(73\right)\end{array}$

To illustrate some specific quantities we define a general complex “data” vector:

f ⁿ=(z¹,z²,z³,z⁴),  (74)

and then the DFT transformation (9) is given by

$\begin{array}{cc}{F}^m=W_{·n}^mf^n=\frac{1}{2}\left(\begin{array}{c}{z}^1+z^2+z^3+z^4\\ z^1-{\mathrm{iz}}^2-z^3+{\mathrm{iz}}^4\\ z^1-z^2+z^3-z^4\\ z^1+{\mathrm{iz}}^2-z^3-{\mathrm{iz}}^4\end{array}\right).& \left(75\right)\end{array}$

Applying the inverse DFT transformation shows that as expected:

$\begin{array}{cc}\begin{array}{c}{\tilde{W}}_{·m}^nF^m=\frac{1}{4}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& i& -1& -i\\ 1& -1& 1& -1\\ 1& -i& -1& i\end{array}\right)\left(\begin{array}{c}{z}^1+z^2+z^3+z^4\\ z^1-{\mathrm{iz}}^2-z^3+{\mathrm{iz}}^4\\ z^1-z^2+z^3-z^4\\ z^1+{\mathrm{iz}}^2-z^3-{\mathrm{iz}}^4\end{array}\right)\\ =\left(z^1,z^2,z^3,z^4\right)\\ =f^n.\end{array}& \left(76\right)\end{array}$

Other quantities are derived from the fⁿ using the a_kl of (72), e.g., consider

$\begin{array}{cc}{f}_{\stackrel{\_}{k}}=a_{\stackrel{\_}{k}n}f^n=\frac{1}{2}\left(\begin{array}{c}{z}_1-{\mathrm{iz}}_2+2{\mathrm{iz}}_3-{\mathrm{iz}}_4\\ {\mathrm{iz}}_1+z_2+{\mathrm{iz}}_3\\ -2{\mathrm{iz}}_1-{\mathrm{iz}}_2+z_3-{\mathrm{iz}}_4\\ {\mathrm{iz}}_1+{\mathrm{iz}}_3+z_4\end{array}\right),& \left(77\right)\end{array}$

and we further check that f_n=, by calculating:

$\begin{array}{cc}\begin{array}{c}{a}_{\stackrel{\_}{k}n}\stackrel{\_}{\left(f^k\right)}=a_{\stackrel{\_}{k}n}{\stackrel{\_}{f}}^{\stackrel{\_}{k}}\\ ={\stackrel{\_}{f}}_n\\ =\frac{1}{2}\left(\begin{array}{c}{\stackrel{\_}{z}}_1+i{\stackrel{\_}{z}}_2-2i{\stackrel{\_}{z}}_3+i{\stackrel{\_}{z}}_4\\ -i{\stackrel{\_}{z}}_1+{\stackrel{\_}{z}}_2-i{\stackrel{\_}{z}}_3\\ -2i{\stackrel{\_}{z}}_1+i{\stackrel{\_}{z}}_2+{\stackrel{\_}{z}}_3+i{\stackrel{\_}{z}}_4\\ -i{\stackrel{\_}{z}}_1-i{\stackrel{\_}{z}}_3+{\stackrel{\_}{z}}_4\end{array}\right).\end{array}& \left(78\right)\end{array}$

Comparing f_n in (77) to the complex conjugation of f_n in (78) we see that as expected:

f _n =,   (79)

demonstrating consistency using the hybrid-index notation. Correspondingly, we can lower indices on the F^m to get

$\begin{array}{cc}\begin{array}{c}{F}_{\stackrel{\_}{k}}=a_{\stackrel{\_}{k}m}F^m\\ =\frac{1}{2}\left(\begin{array}{c}{z}^1+\left(1-2i\right)z^2+\left(1+4i\right)z^3+\left(1-2i\right)z^4\\ \left(1+2i\right)z^1-{\mathrm{iz}}^2-\left(1-2i\right){\mathrm{iz}}^3+{\mathrm{iz}}^4\\ \left(1-4i\right)z^1-\left(1+2i\right)z^2+z^3-\left(1+2i\right)z^4\\ \left(1+2i\right)z^1+{\mathrm{iz}}^2+\left(2i-1\right)z^3-{\mathrm{iz}}^4\end{array}\right).\end{array}& \left(80\right)\end{array}$

Parseval's theorem is invariant as expected:

$\begin{array}{cc}{\underset{2}{a}}_{\stackrel{\_}{m}n}{\stackrel{\_}{F}}^{\stackrel{\_}{m}}F^n={\underset{2}{a}}_{\stackrel{\_}{k}l}{\stackrel{\_}{f}}^{\stackrel{\_}{k}}f^l={z^1}^2+{z^2}^2+{z^3}^2+{z^4}^2,& \left(81\right)\end{array}$

but the Euclidean form in (81) deriving from this particular

${\underset{2}{a}}_{\stackrel{\_}{m}n}$

is not necessarily expected. As pointed out earlier, Parseval's theorem can take different forms for the different a_kl solutions. E.g., considering solution 4 of the a_kl in (71) we have

$\begin{array}{cc}\begin{array}{c}{\underset{4}{a}}_{\stackrel{\_}{m}n}{\stackrel{\_}{F}}^{\stackrel{\_}{m}}F^n={\underset{4}{a}}_{\stackrel{\_}{k}l}{\stackrel{\_}{f}}^{\stackrel{\_}{k}}f^l\\ =-{z^1}^2+{z^2}^2+{z^3}^2+{z^4}^2-2z^1\left(z^2+z^3+z^4\right),\end{array}& \left(82\right)\end{array}$

which is invariant but not real valued.

Eigenvalue-Eigenvector Problem

Given the formalism above, we can examine the eigenvalue-eigenvector (eigen-problem) from a more general viewpoint than what is normally discussed. The additional generality stems from the possibility that the metric is not necessarily the Kronecker delta, constrained only by its invariance under unitary transformations and hermitian symmetry. We have shown from a general viewpoint that the a_kl≐δ_kl solution is an important “base” solution without additional worries that perhaps something was overlooked. We have also shown that other solutions for the a_kl exist based simply on unitary invariance, and thus, more general calculations of the DFT and its various quantities are easily extended to arbitrary basis sets. Some advantages in calculation may be possible in these other general basis sets and we will reserve this topic for a future set of notes. For now, we examine a specific basis set that follows by solution to the eigen-problem as discussed below.

Before examining the specifics of the eigen-problem for the DFT, some generalities are discussed. The eigenvalue-eigenvector problem seeks solutions to the following set equations:

$\begin{array}{cc}{T}_{\stackrel{\_}{m}n}{\underset{a}{v}}^n=\underset{a}{\sigma }{\underset{a}{v}}_{\stackrel{\_}{m}},& \left(83\right)\end{array}$

for the vectors

${\underset{a}{v}}^n$

of a general rank-2 quantity, T_mn. The scalar values,

$\underset{a}{\sigma },$

are the eigenvalues and the

${\underset{a}{v}}^n$

are the eigenvectors. The “a” index below the kernel object, in this case

$,\underset{a}{v}$

and

$\underset{a}{\sigma }$

does not refer to a component, but rather, the a-th enigenvector,

$\underset{a}{v}$

corresponding to the a-th eigenvalue,

$\underset{a}{\sigma }$

as separate and distinct solutions. We therefore choose indices from the set, {a, b, c, d}, to label the kernel objects corresponding to distinct eigen-solutions, and then reserve the index range, {k, l, m, n}, to label components of a given basis.

Continuing to the DFT application we substitute W_mn of (26) into (83):

$\begin{array}{cc}{W}_{\stackrel{\_}{m}n}{\underset{a}{v}}^n=\underset{a}{\sigma }{\underset{a}{v}}_{\stackrel{\_}{m}},& \left(84\right)\end{array}$

and make special note that we have switched to a common kernel symbol for the Fourier transform pair, (F^m, fⁿ)→(v_m, vⁿ), with the substitutions:

v ⁿ ≡f ⁿ and F_m≡v_m.   (85)

Looking ahead, we mention now that the solution to (84) results in a transformation of the W_mn into a diagonal basis by use of the similarity transform. In the context of the present application, and noting expression (9) for the Fourier transform, (84) then has an interesting geometrical interpretation as solving for the Fourier transform pair, (v_m, vⁿ), that are “parallel” to one another, or equivalently, left invariant under the action of the DFT (aside from the scaling factor σ). To give greater emphasis to this last comment, consider the overall structure of the DFT transformation in a diagonal basis for the 4×4 case (for now using the notation of (9) to emphasize this last point):

$\begin{array}{cc}\begin{array}{c}{W}_{\stackrel{\_}{m}n}{\underset{a}{v}}^n=\underset{a}{\sigma }{\underset{a}{v}}_{\stackrel{\_}{m}}\Rightarrow \left(\begin{array}{cccc}\underset{1}{\sigma }& 0& 0& 0\\ 0& \underset{2}{\sigma }& 0& 0\\ 0& 0& \underset{3}{\sigma }& 0\\ 0& 0& 0& \underset{4}{\sigma }\end{array}\right)\left(\begin{array}{c}{f}^{1^{\prime }}\\ f^{2^{\prime }}\\ f^{3^{\prime }}\\ f^{4^{\prime }}\end{array}\right)\\ =\left(\begin{array}{c}\underset{1}{\sigma }f^{1^{\prime }}\\ \underset{2}{\sigma }f^{2^{\prime }}\\ \underset{3}{\sigma }f^{2^{\prime }}\\ \underset{4}{\sigma }f^{2^{\prime }}\end{array}\right)\\ =\underset{a}{\sigma }\left(\begin{array}{c}{F}_{1^{\prime }}\\ F_{2^{\prime }}\\ F_{3^{\prime }}\\ F_{4^{\prime }}\end{array}\right),\end{array}& \left(86\right)\end{array}$

where the primes on indices denote the data and its transform in the diagonal basis as discussed further below. Thus as mentioned above, the Fourier transform pair are co-linear when calculated in the diagonal basis:

$\begin{array}{cc}\left(\begin{array}{c}{F}_{1^{\prime }}\\ F_{2^{\prime }}\\ F_{3^{\prime }}\\ F_{4^{\prime }}\end{array}\right)=\left(\begin{array}{c}\left(\underset{1}{\sigma }/\underset{a}{\sigma }\right)f^{1^{\prime }}\\ \left(\underset{2}{\sigma }/\underset{a}{\sigma }\right)f^{2^{\prime }}\\ \left(\underset{3}{\sigma }/\underset{a}{\sigma }\right)f^{3^{\prime }}\\ \left(\underset{4}{\sigma }/\underset{a}{\sigma }\right)f^{4^{\prime }}\end{array}\right),& \left(87\right)\end{array}$

differing only by a scalar factor. This example illustrates the geometrical significance of the DFT transform in a diagonal basis and motivates using the same kernel symbol for the Fourier transform pair in (85), since the quantities are co-linear. This difference in notation is to be contrasted with the approach taken so far where the kernel symbols have been kept distinct, and the W_·n^m transformation was regarded as a point transformation.

Continuing with the eigen-problem calculation we use the metric tensor to show that (84) is equivalently expressed in the form:

$\begin{array}{cc}\left(W_{\stackrel{\_}{m}n}-\underset{a}{\sigma }a_{\stackrel{\_}{m}n}\right){\underset{a}{v}}^n=0,& \left(88\right)\end{array}$

or also as

$\begin{array}{cc}\left(W_{·n}^m-\underset{a}{\sigma }{\delta }_{·n}^m\right){\underset{a}{v}}^n=0.& \left(89\right)\end{array}$

And thus the eigen-problem, when cast in the form of (89), does not depend on a specific choice for the metric. Thus we solve the DFT eigen-problem in the form of (89). The eigenvalues for the W_·n^m are then solutions to

$\begin{array}{cc}\det \left(W_{·n}^m-\underset{a}{\sigma }{\delta }_{·n}^m\right)=0.& \left(90\right)\end{array}$

Eigen-Problem for the DFT Matrix

The eigenvector-eigenvalue problem for the DFT matrix in the form of (89) has been examined extensively and so we review some of these basic results before examining the diagonal DFT calculation in detail. Since

W⁴=I,  (91)

the eigenvalues for the DFT matrix for N>3 are non-unique and are multiples of the so-called “roots of unity:”

$\begin{array}{cc}\underset{a}{\sigma }={}^{\mathrm{\pi }a/2};a=1,\dots ,4;\Rightarrow \underset{a}{\sigma }\in \left\{\pm 1,\pm i\right\}.& \left(92\right)\end{array}$

Specifically, (92) then becomes

$\begin{array}{cc}\underset{1}{\sigma }=1,\underset{2}{\sigma }=-1,\underset{3}{\sigma }=i,\underset{4}{\sigma }=-i.& \left(93\right)\end{array}$

The eigenvalues of the DFT transform are not all distinct for N≧4, so there are many possible sets of eigenvectors satisfying (89). As a consequence of having these non-unique eigenvalues, the corresponding eigenvectors for these degenerate eigenvalues are not necessarily orthogonal. The importance of having an orthogonal set of eigenvectors for the present application is not known, but should this turn out to be an issue, the Gramm-Schmidt procedure can always be applied or an orthogonal set can be obtained. Further consequences and examples of this analysis for higher dimensional DFTs will be discussed in a future set of notes. For now we examine a low-dimensional example to map out some of the basic features.

Eigenvalues

Considering M=N=4, the DFT transform and its inverse is given in (68) which is reproduced below:

$W_{·k}^m=\frac{1}{2}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& -i& -1& i\\ 1& -1& 1& -1\\ 1& i& -1& -i\end{array}\right);$ ${\stackrel{-1}{W}}_{·k}^m=W{\%}_{·k}^m={\stackrel{\_}{W}}_k^{·m}=\frac{1}{2}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& i& -1& -i\\ 1& -1& 1& -1\\ 1& -i& -1& i\end{array}\right).$

Considering solutions of (90) we have

$\det \left(W_{·n}^m-\underset{a}{\sigma }{\delta }_{·n}^m\right)=0,$

and for the 4×4 example has the solution:

$\begin{array}{cc}{\left(\underset{a}{\sigma }-1\right)}^2\left(\underset{a}{\sigma }+1\right)\left(\underset{a}{\sigma }+i\right)=0\Rightarrow \left\{\underset{1}{\sigma }=1,\underset{2}{\sigma }=-1,\underset{3}{\sigma }=1,\underset{4}{\sigma }=-i\right\},& \left(94\right)\end{array}$

with one degenerate eigenvalue,

$\underset{1}{\sigma }=\underset{3}{\sigma }=1.$

Multiplicity of Eigenvalues for the DFT Matrix

As discussed above, the significance of (91) for the DFT eigenvalue problem is that degenerate eigenvalues appear for N≧4. The “multiplicities” of these eigenvalues are defined as the number of times each eigenvalue appears. For the 4×4 example in (94) we therefore have

m ₁=2,m_2′=1,m₃=1,m₄=0,   (95)

where the m_k ordering corresponds to the ordering of (93). These results can be generalized to higher dimensions by the following Table of multiplicities for each eigenvalue:

TABLE 1
M = N	m1$\mathrm{for}\underset{1}{\sigma }=1\text{:}$	m2$\mathrm{for}\underset{2}{\sigma }=-1\text{:}$	m3$\mathrm{for}\underset{3}{\sigma }=i\text{:}$	m4$\mathrm{for}\underset{4}{\sigma }=-i\text{:}$
4m	m + 1	m	m	m − 1
4m + 1	m + 1	m	m	m
4m + 2	m + 1	m + 1	m	m
4m + 3	m + 1	m + 1	m + 1	m

By inspection of Table 1, the 4×4 case corresponds to the m=1 value in the top row giving the multiplicities of eigenvalues listed in (95). Similarly, it can be shown that for the 5×5 case, the eigenvalues are solutions of

$\begin{array}{cc}{\left(\underset{a}{\sigma }-1\right)}^2\left(\underset{a}{\sigma }+1\right)\left(\underset{a}{{\sigma }^2}+1\right)=0\Rightarrow \left\{\underset{1}{\sigma }=1,\underset{2}{\sigma }=-1,\underset{3}{\sigma }=i,\underset{4}{\sigma }=-i,\underset{5}{\sigma }=1\right\},& \left(96\right)\end{array}$

corresponding to the 2^nd row of the Table with m=1. From (96) (or equivalently Table 1) the multiplicities for the 5×5 case are thus

m ₁=2,m₂=1,m₃=1,m₄=1.  (97)

To illustrate a higher dimensional example, consider say M=N=131=4.32+3, corresponding to the third row of the Table giving the following multiplicities:

m ₁=33,m₂=22,m₃=22,m₄=32.   (98)

Eigenvectors

Given the eigenvalues in (94) the eigenvectors corresponding to each eigenvalue are constructed using (89). It is instructive to consider this process for the eigenvectors corresponding to the degenerate first and third eigenvalues,

$\underset{1}{\sigma }=\underset{3}{\sigma }=1.$

The first eigenvector is obtained as a solution to the set of linear equations:

$\begin{array}{cc}\left(W_{·n}^m-\underset{a}{\sigma }{\delta }_{·n}^m\right){\underset{a}{v}}^n{}_{a\to 1}=0\Rightarrow \left(\begin{array}{cccc}1-\underset{1}{\sigma }& 1& 1& 1\\ 1& -i-\underset{1}{\sigma }& -1& i\\ 1& -1& 1-\underset{1}{\sigma }& -1\\ 1& i& -1& -i-\underset{1}{\sigma }\end{array}\right)\left(\begin{array}{c}{\underset{1}{v}}^1\\ {\underset{1}{v}}^2\\ {\underset{1}{v}}^3\\ {\underset{1}{v}}^4\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right),& \left(99\right)\end{array}$

which simplifies to

$\begin{array}{cc}\left(\begin{array}{c}-{\underset{1}{v}}^1+{\underset{1}{v}}^2+{\underset{1}{v}}^3+{\underset{1}{v}}^4\\ {\underset{1}{v}}^1-\left(2+i\right){\underset{1}{v}}^2-{\underset{1}{v}}^3i{\underset{1}{v}}^4\\ {\underset{1}{v}}^1-{\underset{1}{v}}^2-{\underset{1}{v}}^3-{\underset{1}{v}}^4\\ {\underset{1}{v}}^1+i{\underset{1}{v}}^2-{\underset{1}{v}}^3-\left(2+i\right){\underset{1}{v}}^4\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right).& \left(100\right)\end{array}$

The solution to (100) for

${\underset{1}{v}}^n$

is therefore given in terms of two arbitrary components, say

${\underset{1}{v}}^3\mathrm{and}{\underset{1}{v}}^4\text{:}$

$\begin{array}{cc}\left({\underset{1}{v}}^1,{\underset{1}{v}}^2,{\underset{1}{v}}^3,{\underset{1}{v}}^4\right)=\left({\underset{1}{v}}^3+2{\underset{1}{v}}^4,{\underset{1}{v}}^4,{\underset{1}{v}}^3,{\underset{1}{v}}^4\right).& \left(101\right)\end{array}$

We point out that two components are required because of the two-fold degeneracy of the

$\underset{1}{\sigma }=\underset{3}{\sigma }=1$

eigenvalue. In general, the eigenvectors corresponding to unique eigenvalues are specified in terms of a single component. Continuing we arbitrarily choose

${\underset{1}{v}}^4=1$

and

${\underset{1}{v}}^3=0,$

to obtain a realization for the first eigenvector:

$\begin{array}{cc}\begin{array}{c}\left({\underset{1}{v}}^1,{\underset{1}{v}}^2,{\underset{1}{v}}^3,{\underset{1}{v}}^4\right)=\left({\underset{1}{v}}^3+2{\underset{1}{v}}^4,{\underset{1}{v}}^4,{\underset{1}{v}}^3,{\underset{1}{v}}^4\right){}_{\underset{{\underset{1}{v}}^3\to 0}{{\underset{1}{v}}^4\to 1}}\Rightarrow \underset{1}{v^n}\\ =\left(2,1,0,1\right).\end{array}& \left(102\right)\end{array}$

Note also that since,

$\underset{1}{\sigma }=\underset{3}{\sigma }=1,$

are degenerate, the solution to (99) is trivially the same as (100). Therefore, a solution for

${\underset{3}{v}}^n$

follows immediately by choosing complementary but arbitrary numerical values for the

${\underset{1}{v}}^3$

and

${\underset{1}{v}}^4$

in (102), e.g.,

${\underset{3}{v}}^4=0$

and

$\begin{array}{cc}{\underset{3}{v}}^3=1:\left({\underset{3}{v}}^1,{\underset{3}{v}}^2,{\underset{3}{v}}^3,{\underset{3}{v}}^4\right)=\left({\underset{3}{v}}^3+2{\underset{3}{v}}^4,{\underset{3}{v}}^4,{\underset{3}{v}}^3,{\underset{3}{v}}^4\right){}_{{\underset{1}{v}}^3\to 1}^{\text{?}}\Rightarrow {\underset{3}{v}}^n=\left(1,0,1,0\right).\text{?}\text{indicates text missing or illegible when filed}& \left(103\right)\end{array}$

Similarly, we can solve for the remaining eigenvectors corresponding to their eigenvalues,

$\underset{2}{\sigma }=-1$

and

$\underset{4}{\sigma }=-i.$

In summary, a set of engenvectors (but not the only set) for the 4×4 DFT transformation is summarized below:

$\begin{array}{cc}{\underset{1}{v}}^n=\left(\begin{array}{c}2\\ 1\\ 0\\ 1\end{array}\right),{\underset{2}{v}}^n=\left(\begin{array}{c}-1\\ 1\\ 1\\ 1\end{array}\right),{\underset{3}{v}}^n=\left(\begin{array}{c}1\\ 0\\ 1\\ 0\end{array}\right),{\underset{4}{v}}^n=\left(\begin{array}{c}0\\ -1\\ 0\\ 1\end{array}\right).& \left(104\right)\end{array}$

Since two of the eigenvalues are degenerate, (102) and (104) collectively do not form an orthogonal system as can be verified by forming the inner product between any two vectors using the “base” solution,

${\underset{1}{a}}_{\stackrel{\_}{k}l}\stackrel{*}{=}{\delta }_{\stackrel{\_}{k}l}$

of (71):

$\begin{array}{cc}〈\underset{m}{v}\underset{n}{v}〉={\underset{1}{a}}_{\stackrel{\_}{k}l}{\stackrel{\_}{\underset{m}{v}}}^{\stackrel{\_}{k}}{\underset{n}{v}}^l={\delta }_{\stackrel{\_}{k}l}{\stackrel{\_}{\underset{m}{v}}}^{\stackrel{\_}{k}}{\underset{n}{v}}^l=\left(\begin{array}{cccc}3& 0& 1& 0\\ 0& 2& 0& 0\\ 1& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right),& \left(105\right)\end{array}$

where the off-diagonal entries correspond to the

$\underset{1}{\sigma }=\underset{3}{\sigma }=1$

value. For reference, we can still construct an orthonormal set from (104):

$\begin{array}{cc}{\underset{1}{u}}^n=\frac{1}{\sqrt{6}}\left(\begin{array}{c}2\\ 1\\ 0\\ 1\end{array}\right),{\underset{2}{u}}^n=-\frac{1}{2}\left(\begin{array}{c}-1\\ 1\\ 1\\ 1\end{array}\right),{\underset{3}{u}}^n=\frac{1}{2\sqrt{3}}\left(\begin{array}{c}1\\ -1\\ 3\\ -1\end{array}\right),{\underset{4}{u}}^n=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\ -1\\ 0\\ 1\end{array}\right).& \left(106\right)\end{array}$

and verify as expected that

$\begin{array}{cc}{\delta }_{\stackrel{\_}{k}l}{\stackrel{\_}{\underset{m}{u}}}^{\stackrel{\_}{k}}{\underset{n}{u}}^l={\delta }_{\mathrm{mn}}.& \left(107\right)\end{array}$

Similarity Transform

Given the eigenvectors of (104) we form the “similarity” transform using the eigenvectors as columns. We denote this as a transformation of basis and we find by construction that the similarity transform, S_m^m′, and its inverse,

${\stackrel{-1}{S}}_m^{m^{\prime }},$

are given by

$\begin{array}{cc}{S}_m^{m^{\prime }}=\left(\begin{array}{cccc}2& -1& 1& 0\\ 1& 1& 0& -1\\ 0& 1& 1& 0\\ 1& 1& 0& 1\end{array}\right);{\stackrel{-1}{S}}_m^{m^{\prime }}=\frac{1}{4}\left(\begin{array}{cccc}1& 1& -1& 1\\ -1& 1& 1& 1\\ 1& -1& 3& -1\\ 0& -2& 0& 2\end{array}\right).& \left(108\right)\end{array}$

When using the kernel-index method for transformations of basis, say x^m′=A_m^m′x^m. that the A_m^m′ quantities do not use dots above or below an index to mark index placement. This differs from the notation for point transformations, which does use the dots e.g., F^m=W_·n^mfⁿ. The reason for this is to discourage the “raising” or “lowering” of indices on the A_m^m′ to avoid an inconsistency when using the metric tensor. This inconsistency can occur because the metric a_kl also transforms under the A_m^m′, and thus has different components in the new coordinate system. For point transformations however, no coordinate transformation occurs and thus indices can be raised or lowered using a single (unambiguous) metric.

We then verify that applying the similarity transformation to the W_·n^m gives as expected:

$\begin{array}{cc}{\mathrm{SWS}}^{-1}\equiv S_m^{m^{\prime }}W_{·n}^m{\stackrel{-1}{S}}_{n^{\prime }}^n=W_{·n^{\prime }}^{m^{\prime }}=\left(\begin{array}{cccc}\underset{1}{\sigma }& 0& 0& 0\\ 0& \underset{2}{\sigma }& 0& 0\\ 0& 0& \underset{3}{\sigma }& 0\\ 0& 0& 0& \underset{4}{\sigma }\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& -i\end{array}\right).& \left(109\right)\end{array}$

As a check on index ordering we run the following line of Mathematica code with Si defined as the inverse of the Sin (108):

Table [Sum[S[[m,ump]]W[[m,n]]Si[[lnp,n]], {n,l,M}, {m,l,M}], {ump,l,M}, {lnp,l,M}]  (110)

where ump denotes “upper-m-prime” and lnp is “lower-m-prime.”

It is of further interest to point out that if the similarity transform is constructed using the orthonormal eigenvectors of (106) as columns (and correspondingly it inverse), then the application of the “orthonormal” similarity transform to the W_·n^m also gives (109). Therefore, the similarity transform that leaves the W_·n^m in diagonal form is certainly not unique, so many possible choices can be made for the S_n^m′. So as indicated by calculations thus far, it does not seem to matter whether or not the eigenvectors forming the S_n^m′ are orthogonal, and this can greatly simplify the initial calculations needed for calculating in the diagonal basis.

DFT Calculation in Diagonal Form

The diagonal form of the DFT in (109) shows that the DFT calculation can be performed both efficiently and elegantly as simply a multiplication of the entries on the main diagonal (the eigenvalues), to the corresponding data components that have been pre-multiplied using the similarity matrix. This means that after pre-multiplying the data, the DFT transform is reduced to a vector product consisting of N operations rather than a matrix multiply with N² operations (or N log N when using the FFT). In reality though, this reduction in computation cost is only partly true in terms of the performance gains, i.e., the diagonal DFT representation does obviously lead to a speedy calculation for one part of the DFT operation, but in general, this gain is offset by the added computational cost of having to pre-transform the data and then back-transform to bring the Fourier transform into the original basis as illustrated below. But the key is that these latter operations may not need to be performed at every step in an iterative transform application, such as with phase retrieval. Additional details on this calculation approach are outlined below.

To begin, we write down the DFT calculation in the diagonal basis after switching the notation back to distinct kernel symbols for the data and its transform:

f ^m ′=W _·n ^m ′f ⁿ′,   (111)

where the m′ represents the transformation of basis under the action of the similarity transform as illustrated above in (109). Since the only non-zero values for the W_·n′^m′ are along the diagonal, we form a vector,

$\underset{m^{\prime }}{\sigma },$

form the eigenvalues (i.e., the diagonal entries of the W_·n′^m′) as in (94). Continuing with the calculation of (111, the starting data, f_k, are pre-multiplied using the inverse similarity transform,

${\stackrel{-1}{S}}_n^{n^{\prime }},$

to give

$\begin{array}{cc}{f}^{n^{\prime }}={\stackrel{-1}{S}}_n^{n^{\prime }}f^n.& \left(112\right)\end{array}$

The diagonal DFT is then the product of the N components of the

$\underset{m^{\prime }}{\sigma }$

and the fⁿ′ defined by (112):

$\begin{array}{cc}{F}^{m^{\prime }}=\underset{m^{\prime }}{\sigma }·f_{m^{\prime }},& \left(113\right)\end{array}$

where no sum is implied over the (m′,m′) pair. To complete the calculation, the “diagonal” DFT components in (113) need to be transformed back to the original DFT basis, F^m. To obtain this result we express the LHS of (111) in terms of components in the original basis:

$\begin{array}{cc}{F}^{m^{\prime }}={\stackrel{-1}{S}}_m^{m^{\prime }}F^m,& \left(114\right)\end{array}$

and therefore the original Fourier transform components are readily obtained by back-transforming the F^m′ of (114) using the similarity transform:

$\begin{array}{cc}{F}^n=S_{n^{\prime }}^nF^{n^{\prime }}=S_{n^{\prime }}^n{\stackrel{-1}{S}}_m^{n^{\prime }}F^m={\delta }_m^nF^m.& \left(115\right)\end{array}$

The inverse DFT is also calculated in a diagonal basis since the inverse of a diagonal matrix is found by populating its diagonal entries with the inverse of the individual components. Therefore, given the W_·n^m′,the inverse is simply (no sum is implied on the RHS over the (m′, m′) pair):

$\begin{array}{cc}{\stackrel{-1}{W}}_{•n^{\prime }}^{m^{\prime }}=\left(1/\underset{m^{\prime }}{\sigma }\right){\delta }_{n^{\prime }}^{m^{\prime }},& \left(116\right)\end{array}$

and then inverse DFT transformation in the diagonal basis gives as expected:

$\begin{array}{cc}\left(1/\underset{m^{\prime }}{\sigma }\right)·F^{m^{\prime }}=\left(1/\underset{m^{\prime }}{\sigma }\right)·\underset{m^{\prime }}{\sigma }·f^{m^{\prime }}=f^{m^{\prime }},& \left(117\right)\end{array}$

The overall point is that for applications where it is feasible to perform the DFT in a diagonal basis, the similarity transform and its inverse can be pre-calculated. The multiplication of the similarity transform with the data can also be pre-calculated. Then at an appropriate point, the DFT results can be transformed back to the original basis using the inverse similarity transform, preferably at the end of an iterative loop. The approach is summarized in FIG. 2.

Example: 4×4 Diagonal DFT Calculation

The steps above are illustrated for the 4×4 example. Consider the general complex data vector of (74):

f ⁿ=(z¹,z²,z³,z⁴).  (118)

The inverse similarity transform of (108) times the data vector (118) is given by

$\begin{array}{cc}{f}^{m^{\prime }}={\stackrel{-1}{S}}_m^{m^{\prime }}f^m=\frac{1}{4}\left(\begin{array}{c}{z}_1+z_2-z_3+z_4\\ -z_1+z_2+z_3+z_4\\ z_1-z_2+3z_3-z_4\\ 2\left(-z_2+z_4\right)\end{array}\right),& \left(119\right)\end{array}$

and then from (94):

$\begin{array}{cc}\underset{m^{\prime }}{\sigma }=\left(1,-1,1,-i\right).& \left(120\right)\end{array}$

The DFT in the transformed basis follows by multiplication of the eigenvalues in (120) with the f^m′ components in (119) to give the Fourier transform in the diagonal basis:

$\begin{array}{cc}{F}^{m^{\prime }}=\underset{m^{\prime }}{\sigma }·f^m=\frac{1}{4}\left(\begin{array}{c}{z}_1+z_2-z_3+z_4\\ z_1-z_2-z_3-z_4\\ z_1-z_2+3z_3-z_4\\ -2i\left(-z_2+z_4\right)\end{array}\right).& \left(121\right)\end{array}$

Finally, we back-transform the F^m′ components in (121) by multiplying with the similarity transformation of (108) to give the DFT result in the original non-diagonal basis:

$\begin{array}{cc}\begin{array}{c}{F}^m=S_{m^{\prime }}^mF^{m^{\prime }}\\ =\left(\begin{array}{cccc}2& -1& 1& 0\\ 1& 1& 0& -1\\ 0& 1& 1& 0\\ 1& 1& 0& 1\end{array}\right)\frac{1}{4}\left(\begin{array}{c}{z}_1+z_2-z_3+z_4\\ z_1-z_2-z_3-z_4\\ z_1-z_2+3z_3-z_4\\ -2i\left(-z_2+z_4\right)\end{array}\right)\\ =\frac{1}{2}\left(\begin{array}{c}{z}^1+z^2+z^3+z^4\\ z^1-{\mathrm{iz}}^{2}-z^3+{\mathrm{iz}}^4\\ z^1-z^2+z^3-z^4\\ z^1+{\mathrm{iz}}^2-z^3-{\mathrm{iz}}^4\end{array}\right),\end{array}& \left(122\right)\end{array}$

in agreement with the original DFT transform, F^m=W_·n^m fⁿ of (75). Similarly, the inverse transformation can be calculated in the diagonal basis as discussed above using (117).

The diagonal form of the DFT in (113) shows that the calculation can be performed, efficiently and elegantly as simply a vector product. This means that the DFT can be reduced to N multiplications of the data with the eigenvalues. But as emphasized above, this observation is only partly true, i.e., although the performance gains from the diagonal representation are indeed realized, these are offset by the computational cost of setting up the data in a transformed basis as illustrated above. But even with this in mind, certain applications spend most of the time in the forward and inverse transform part of the calculations (see FIG. 2), so the main point is to explore applications where the initial and final data transformations are needed only once, and then subsequent DFT calculations can proceed as simply N multiplications of the entries along the main diagonal. The performance gain for a single DFT calculation can then be approximately log(N) times faster. Assuming for simplicity a power of 2 for the data size (N=2ⁿ), we can estimate a performance improvement for large n to be n log(2)˜10¹, or about an order of magnitude faster than the equivalent FFT calculation, and for arbitrary frequency spacing.

FIG. 3 shows an example embodiment of the present invention and is thus directed to an image-based wavefront sensing method of and means for focusing by determining phase difference between received images, removing that phase difference by calculating a DFT for each signal and iteratively multiplying the resultant until a convergence is achieved by accepting point-source stellar objects to recover embedded optical phase information, receiving data from a sensor array wherein the data includes embedded phase information, calculating a similarity transform and pre-transforming the data once at the beginning of the calculation.

FIG. 3 further shows, in an example embodiment of the present invention, an image-based wavefront sensing method of and means for iteratively processing the data using arbitrary frequency spacing in a diagonal fashion until a wavefront convergence is achieved to determine a phase difference by implementing the DFT in a 1 dimensional linear array using just N operations, where 0<N<N×log(N)<∞, wherein N×log(N) is the number of operations of an FFT counterpart, back-transforming the data at the end of the calculation in a single step and focusing the received images based on the determined phase.

While the invention is described with reference to an exemplary embodiment, it will be understood by those skilled in the art that various changes may be made and equivalence may be substituted for elements thereof without departing from the scope of the invention. In addition, many modifications may be made to the teachings of the invention to adapt to a particular situation without departing from the scope thereof. Therefore, it is intended that the invention not be limited the embodiments disclosed for carrying out this invention, but that the invention includes all embodiments falling with the scope of the appended claims. Moreover, the use of the terms first, second, etc. does not denote any order of importance, but rather the terms first, second, etc. are used to distinguish one element from another.

Claims

1. An image-based wavefront sensing method of focusing by determining phase difference between received images, removing that phase difference by calculating a discrete Fourier transform (DFT) for each signal and iteratively multiplying the resultant until a convergence is achieved comprising the steps of: accepting point-source images to recover embedded phase information:receiving data from an external source wherein the data includes embedded phase information;calculating a similarity transform;iteratively processing the data in a diagonal fashion until a wavefront convergence is achieved to determine a phase difference; andfocusing the received images based on the determined phase.

2. The wavefront sensing method of claim 1 wherein calculating a similarity transform further includes pre-transforming the data once at the beginning of the calculation.

3. The wavefront sensing method of claim 2 wherein iteratively processing further includes a single step of back-transforming the data at the end of the calculation,

4. The wavefront sensing method of claim 3 wherein said point source images include stellar objects.

5. The wavefront sensing method of claim 4 wherein said embedded phase information is of an optical nature.

6. The wavefront sensing method of claim 5 wherein said external source includes a sensor array.

7. The wavefront sensing method of claim 6 wherein the step of iteratively processing further includes the step of determining the phase difference by utilizing a discrete Fourier transform (DFT).

8. The wavefront sensing method of claim 7 including implementing said DFT in a 1 dimensional linear array using just N operations, where 0<N<N×log(N)<∞, wherein N×log(N) is the number of operations of an FFT counterpart and N is an integer.

9. The wavefront sensing method of claim 8 including iteratively processing using arbitrary frequency spacing.

10. An image-based wavefront apparatus for focusing by determining phase difference between received images, removing that phase difference by calculating a discrete Fourier transform (DFT) for each signal and iteratively multiplying the resultant until a convergence is achieved comprising: means for accepting point-source images to recover embedded phase information:means for receiving data from an external source wherein the data includes embedded phase information;means for calculating a similarity transform;means for iteratively processing the data in a diagonal fashion until a wavefront convergence is achieved to determine a phase difference; andmeans for focusing the received images based on the determined phase.

11. The wavefront sensing apparatus of claim 10 wherein said means for calculating a similarity transform further includes means for pre-transforming the data once at the beginning of the calculation.

12. The wavefront sensing apparatus of claim 11 wherein said means for iteratively processing further includes means for of back-transforming the data at the end of the calculation in a single step,

13. The wavefront sensing apparatus of claim 12 wherein said point source images include stellar objects.

14. The wavefront sensing apparatus of claim 13 wherein said embedded phase information is of an optical nature.

15. The wavefront sensing apparatus of claim 14 wherein said external source includes a sensor array.

16. The wavefront sensing apparatus of claim 15 wherein the means for iteratively processing further includes means for determining the phase difference by utilizing a discrete Fourier transform (DFT).

17. The wavefront sensing apparatus of claim 16 further including means for implementing said DFT in a 1 dimensional linear array using just N operations, where 0<N<N×log(N)<∞, wherein N×log(N) is the number of operations of an FFT counterpart.

18. The wavefront sensing apparatus of claim 17 further including means for iteratively processing using arbitrary frequency spacing.