Fast Infinitesimal Fourier Transform for Signal and Image Processing via Multiparametric and Fractional Fourier Transforms

Abstract

The fractional Fourier transforms (FrFTs) is one-parametric family of unitary transformations {F α}2απ=0. FrFTs found a lot of applications in signal and image processing. The identical and classical Fourier transformations are both the special cases of the FrFTs. They correspond to α = 0 (F 0 = I) and α = π/2 (F π/2 = F ), respectively. Up to now, the fractional Fourier spectra F αi = F αi {f } , i = 1, 2, ..., M , has been digitally computed using classical approach based on the fast discrete Fourier transform. This method maps the N samples of the original function f to the N samples of the set of spectra {F αi }iM=1 , which requires M N (2 + log2 N ) multiplications and M N log2 N additions. This paper develops a new numerical algorithm, which requires 2M N multiplications and 3M N additions and which is based on the infinitesimal Fourier transform.

Keywords

• Fast fractional Fourier transform
• infinitesimal Fourier transform
• Schr¨odinger operator
• signal and image analysis

-

The idea of fractional powers of the Fourier operator {F a}4a=0 appeared in the mathematical literature [ 1,2,3,4 ]. The idea is to consider the eigen-value decomposition of the Fourier transform F in terms of the eigen-values λn = ejnπ/2 and eigen-functions in the form of the Hermite functions. The family of FrFT a 4 {F }a=0 is constructed by replacing the n-th eigen-value λn = ejnπ/2 by its a-th power λan = ejnπa/2 for a between 0 and 4. This value is called the transform order. There is the angle parameterization {F α}2απ=0 , where α = πa/2 is a new angle parameter. Since this family depends on a single parameter, the fractional operators {F }a=0 (or {F α}2απ=0) form the Fourier-Hermite onea 4 a⊕b a b = F 4 (or parameter strongly continuous unitary multiplicative group F F α⊕β α β = F 2π ), where a⊕b = (a + b) mod4 (or α⊕β = (α + β) mod2π) and F F

4 2π F 0 = I. The identical and classical Fourier transformations are both the special cases of the FrFTs. They correspond to α = 0 (F 0 = I) and α = π/2 (F π/2 = F ), respectively.

In 1980, Namias reinvented the fractional Fourier transform (FrFT) again in his paper [ 6 ]. He used the FrFT in the context of quantum mechanics as a way to solve certain problems involving quantum harmonic oscillators. He not only stated the standard definition for the FrFT, but, additionally, developed an operational calculus for this new transform. This approach was extended by McBride and Kerr [ 7 ]. Then Mendlovic and Ozaktas introduced the FrFT into the field of optics [ 8 ] in 1993. Afterwards, Lohmann [ 9 ] reinvented the FrFT based on the Wigner-distribution function and opened the FrFT to bulk-optics applications. It has been rediscovered in signal and image processing [ 10 ]. In these cases, the FrFT allows us to extract time-frequency information from the signal. A recent state of the art can be found in [ 11 ]. In the series of papers [ 12,13,14,15,16 ], we developed a wide class of classical and quantum fractional transforms.

In this paper, the infinitesimal Fourier transforms are introduced, and the relationship of the fractional Fourier transform with the Schr¨odinger operator of the quantum harmonic oscillator is discussed. Up to now, the fractional Fourier spectra F αi = F αi {f } , i = 1, 2, ..., M, have been digitally computed using classical approach based on the fast discrete Fourier transform. This method maps the N samples of the original function f to the N M samples of the set of spectra {F αi }iM=1 , which requires M N (2 + log2 N ) multiplications and M N log2 N additions. This paper develops a new numerical algorithm, which requires 2M N multiplications and 3M N additions and which is based on the infinitesimal Fourier transform. 2

Eigen-decomposition and Fractional Discrete

Transforms Let F = [Fk (i)]kN,i−=10 be an arbitrary discrete unitary (N × N )-transform, λn and Ψn (t) n = 0, 1, . . . , N − 1 be its eigen-values and eigen-vectors, respectively.

.

Let U = Ψ0(i)|Ψ1(i)|..|ΨN−1(i) be the matrix of the F -transform eigen-vectors. Then U−1·F ·U = Diag {λn}. Hence, we have the following eigen-decomposition: F = [Fk(i)] = U · Λ · U−1 = U · Diag {λn} · U−1.

Definition 1. [ 12 ]. For an arbitrary real numbers a0, . . . , aN−1, we introduce the multi-parametric F -transform

F (a0,...,aN−1) := U diag λ0a0 , . . . , λaNN−−11 U−1.

(1) If a0 = . . . = aN−1 ≡ a then this transform is called fractional F -transform [ 12,13,14,15,16 ]. For this transform we have

F a := U diag λ0a, . . . , λaN−1

U−1 = UΛaU−1.

(2) The zero-th-order fractional F -transform is equal to the identity transform: F 0 = UΛ0U−1 = UU−1 = I , and the first-order fractional Fourier transform operator F 1 T=hFe fiasmeqiluieasl ntoFt(hαe0,.i.n.,iαtNia−l1)Fo-(tαra0,n..s.,fαoNr−m1)F∈R1 N=aUndΛU{F−a1}.a∈R form multi- and one-parameter continuous unitary groups, respectively, with multiplication rules F (a0,...,aN−1)F (b0,...,bN−1) = F (a0+b0,...,aN−1+bN−1) and F aF b = F a+b.

Indeed, F aF b = UΛaU−1 · UΛbU−1 = UΛa+bU−1 = F a+b and F

(a0,...,aN−1)F (b0,...,bN−1) = = U diag λ0a0 , . . . , λaNN−−11

U−1 · U ndiag λb00 , . . . , λbNN−−11 o U−1 = = U ndiag λ0a0+b0 , . . . , λaNN−−11+bN−1 o U−1 = F (a0+b0,...,aN−1+bN−1).

Let F = [Fk (i)]kN,i−=10 be a discrete Fourier (N × N )-transform (DFT), then λn = ejπn/2 ∈ {±1, ±j} , where j = √−1 and {Ψn (t)}nN=−01 are the Kravchuk polynomials.

Definition 2. The multi-parametric and fractional DFT are

F (a0,...,aN−1) := U ndiag ejπ0a0/2, ejπ1a1/2, . . . , ejπ(N−1)aN−1/2 o U−1,

F a := U ndiag ejπna/2 o U−1 and

F (α0,...,αN−1) := U ndiag ej0α0 , ej1α1 , . . . , ej(N−1)αN−1 o U−1,

F α := U diag ejnα U−1 in a- and α-parameterizations, respectively, where α = πa/2.

The parameters (a0, . . . , aaNa−r1e)paenridodaiccianneabcehapnayrarmeaeltevralwuietsh. pHeorwioedve4rs,inthcee operators F (a0,...,aN−1) and F F 4 = I. Hence, F (a0,...,aN−1)F (b0,...,bN−1 = F and F aF b = a⊕b F 4 , where ai⊕bi = (ai + bi) mod4, ∀i = 0, 1, ..., N − 1. Therefore, the ranges 4 of (a0, . . . , aN−1) and a are (Z/4Z)N = [ 0, 4 ]N = [ −2, 2 ]N and Z/4Z = [ 0, 4 ] = [ −2, 2 ], respectively.

In the case of α-parameterization, we have αi ⊕βi = (αi + βi) mod2π, ∀i = 2π 0, 1, ..., N −1. So, the ranges of (α0, . . . , αN−1) and α are (Z/2πZ)N = [0, 2π]N = [−π, π]N and Z/2πZ = [0, 2π] = [−π, π], respectively. (a0⊕b0,...,aN−1⊕4 bN−1)

4

Canonical FrFT

The continuous Fourier transform is a unitary operator F that maps squareintegrable functions on square-integrable ones and is represented on these functions f (x) by the well-known integral

1 Z F (y) = (F f ) (y) = √2π x∈R f (x)e−jyxdx. (3) Relevant properties are that the square F 2f (x) = f (−x) is the inversion operator, and that its fourth power F 4f (x) = f (x) is the identity. Hence, F 3 = F −1. Thus, the operator F generates a cyclic group of the order 4. In 1961, Bargmann extended the Fourier transform in his paper [ 5 ] where he gave definition of the FrFT that was based on the Hermite polynomials as an integral transformation. If Hn (x) is a Hermite polynomial of order n, where Hn (x) = (−1)n ex2 ddxnn ex2 , then for n ∈ N0, functions Ψn(x) = √2n1n!√π Hn(x)e−x2/2 are the eigen-functions of the Fourier transform

F [Ψn (x)] = 1 Z +∞ 2π −∞

Ψn (x) e2πjyxdx = λnΨn (y) = e−j π2 nΨn (y) with λn = jn = e−j π2 n being the eigen-value corresponding to the n-th eigenfunction. According to Bargmann, the fractional Fourier transform F α = [Kα (x, y)] is defined through its eigen-functions as

Kα (x, y) := U diag e−jαn

∞ U−1 = X e−jαnΨn (x) Ψn (y) .

(4) n=0 n=0

Hence,

n=0 Kα (x, y) := X e−jαnΨn (x) Ψn (y) = e−(x2+y2) X∞ e−jαn2Hnnn(!√x)πHn(y) = ∞ 1 = √ √ π 1 − e−2jα · exp ( 2xye−jα − e−2jα x2 + y2 ) 1 − e−2jα

( exp − x2 + y2 ) 2 (5) where Kα (x, y) is the kernel of the FrFT. In the last step we used the Mehler formula [ 19 ] n=0 X∞ e−jαnHn(x)Hn(y) 2nn!√π

1 = √ √ π 1 − e−2jα exp ( 2xye−jα − e−2jα x2 + y2 ) 1 − e−2jα

Expression (5) can be rewritten as

Kα(x, y) = r 1 − j cot α 2π exp (x2 + y2) cos α − 2xy where α 6= πZ (or a 6= 2Z). Obviously, functions Ψn(x) are eigen-functions of the fractional Fourier transform F α [Ψn(x)] = ejnαΨn(x) corresponding to the n-th eigen-values ejnα, n = 0, 1, 2, . . .. The FrFT F α is a unitary operator that maps square-integrable functions f (x) on square-integrable ones F α(y) = (F αf ) (y) =

f (x)Kα(x, y)dx = Z

x∈R e− 2j ( π2 αˆ−α) Z = p2π |sin α| R f (x) exp e− 2j ( π2 αˆ−α)ejy2 2csoisnαα Z hf (x)ej x22 cot αi e−jxydx = = Aα(y) · F {f (x) · Bα(x)} (y), where Aα(y) = e− 2j ( π2 αˆ−α)ejy2 2csoisnαα , Bα(x) = ej x22 cot α.

√2π|sin α|

Let us introduce the uniform discretization of the angle parameter α on M discrete values {α0, α1, ..., αi, αi+1, ..., αM−1} , where αi+1 = αi + Δα, αi = iΔα and Δα = 2π/M.

The set of M spectra {F α0 (y) , F α1 (y) , ..., F αM−1 (y)} can be computed by applying the following sequence of steps for all {α0, α1, ..., αM−1}: 1. Compute products f (x)Bαk (x), which require N multiplications. 2. Compute the Fast Fourier Transform (N log2 N multiplications and additions).

3. Multiply the result by Aα(y) (N multiplications).

This numerical algorithm requires M N log2N additions and M N (2 + log2N ) multiplications. 4 Infinitesimal Fourier Transform In order to construct fast multi-parametric F -transform and fractional Fourier transform algorithms we turn our attention to notion of a semigroup and its generator (infinitesimal operator). Let L2(R, C) be a space of complex-valued functions (signals), and let Op(L2) be the Banach algebra of all bounded linear operators on L2(R, C) endowed with the operator norm. A family {U(α)}α∈R ⊂ Op(L2) is called the Hermite group on L2(R, C) if it satisfies the Abel functional equations U(α + β) = U(α)U(β), α, β ∈ R and U(0) = I, and the orbit maps α → F α = U(α) {f } are continuous from R into L2(R, C) for every f ∈ L2(R, C). Definition 3. The infinitesimal generator A(0) of the group {U(α)}α∈R and the infinitesimal transform U(dα) are defined as follows [ 18,19 ]: A(0) =

, U(dα) = I + dU(0) = I + A(0)dα.

Obviously, But

U(α0 + dα) = U(α0) + dU(α0) = U(α0) + ∂U(α) ∂α α0

dα = = U(α0) + A(α0)dα.

U(α0 + dα) = U(dα0)U(α0) = [I + dU(0)] U(α0) = = U(α0) +

U(α0)dα = ∂U(α) ∂α

α=0 = U(α0) + A(0)U(α0)dα = [I + A(0)] U(α0)dα.

Hence, A(α0) = A(0)U(α0) and F α0+dα(y) = hI + A(0)iF α0 (y)dα. 2 Define now the linear operator H = 12 ddx2 − x2 + 1 . It is known that 1 d2

HΨn(x) = 2 dx2 − x2 + 1 Ψn(x) = nΨn(x).

From (4) and (6) we have

j ∂ ∂α

F α(y) α=0

∂ = j ∂α {F αF } (y) α=0

∞ = X nΨn(y) n=0

Z

R

Ψn(x)f (x)dx, ∞ HF α(x) = X nΨn(y) n=0

Z

R Ψn(x)f (x)dx.

Therefore, j ∂F∂αα(x) = HF α(y), ∂FFαα((xx)) = −jH∂α. The solution of this equation is given by F α(x) = e−jαHF Obviously,

and F α = e−jαH = e−jαh 21 ddx22 −x2+1 i.

F α+dα = F dαF α ' (I + dF α) exp [−jαH] = = I + ∂F α dα exp (−jαH) = (I − jHdα) exp (−jαH) ,

∂α where the operator 1 d2 F dα = (I − jHdα) = I − j 2 dx2 − x2 + 1 dα (6) (7) is called the infinitesimal Fourier transform or the generator of the fractional Fourier transforms [ 17,18 ].

Let us introduce operators (Mxf ) (x) := xf (x) and (MyF ) (y) := yF (y). Using the Fourier transform (3), the first of ones may be written as Mx = F −1 j ddy F . Obviously, x2 = Mx2 = −F 1 ddy22 F . Then

Discretization of x-domain with the interval discretization Δx is equal to the periodization of y-domain = f (n)e−j 2Nπ n − 2f (n) + f (n)ej 2Nπ n = 2f (n) cos 2Nπ n − 1 .

In the N -diagonal basis we have

. .

. . 1/2  . . .  1/2 . .  cos(3Ω) − 3/2 1/2 .  1/2 . . . 1/2  . 1/2 −1/2 f (0) f (1) f (2) f (3) . .

. f (N − 1) where Ω = 2π/N .

Let us introduce the uniform discretization of the angle parameter α on M discrete values {α0, α1, ..., αi, αi+1, ..., αM−1} , where αi+1 = αi + Δα, αi = iΔα and Δα = 2π/M. Then

F αi+1 (y) = F αi+Δα(y) ≈ F αi (k) + jΔα× × cos 2π N k − 3/2 F αi (k) + (f (k − 1) − 2f (k) + f (k + 1)) . More fine approximations O(h2k) also can be used [ 19 ].      ,     (8) 5

Conclusions

In this work, we introduce a new algorithm of computing for Fractional Fourier transforms based on the infinitesimal Fourier transform. It requires 2M N multiplications and 3M N additions vs. M N (2 + log2 N ) multiplications and M N log2 N additions in the classical algorithm. Presented algorithm can be utilized for fast computation in most applications of signal and image processing. We have presented a definition of the infinitesimal Fourier transform that exactly satisfies the properties of the Schrodinger Equation for quantum harmonic oscillator. Acknowledgment. This work was supported by grants the RFBR Nos. 13-0712168, and 13-07-00785.

References

1. Mitra , S.K. : Nonlinear Image Processing . Academic Press Series in Communications, Networking, and Multimedia, San Diego, New York, 248 ( 2001 )

2. Wiener , N.: Hermitian polynomials and Fourier analysis . J. Math. Phys. 8 , 70 - 73 ( 1929 ).

3. Condon , E.U.: Immersion of the Fourier transform in a continuous group of functional transforms . Proc. Nat. Acad. Sci . 23 , 158 - 164 ( 1937 )

4. Kober , H.: Wurzeln aus der Hankel- und Fourier und anderen stetigen Transformationen . Quart. J. Math. Oxford Ser . 10 , 45 - 49 ( 1939 )

5. Bargmann , V.: On a Hilbert space of analytic functions and an associated integral transform . Part 1. Commun. Pure Appl . Math. 14 , 187 - 214 ( 1961 )

6. Namias , V. : The fractional order Fourier transform and its application to quantum mechanics . J. Inst. Math. Appl . 25 , 241 - 265 ( 1980 )

7. McBride , A.C. , Kerr , F.H. : On Namias' fractional Fourier transforms . IMA J. Appl. Math . 39 , 159 - 265 ( 1987 ).

8. Ozaktas , H.M. , Mendlovic , D. : Fourier transform of fractional order and their optical interpretation . Opt. Commun . 110 , 163 - 169 ( 1993 )

9. Lohmann , A.W. : Image rotation, Wigner rotation, and the fractional order Fourier transform . J. Opt. Soc. Am. A . 10 , 2181 - 2186 ( 1993 )

10. Almeida , L.B. : The fractional Fourier transform and time-frequency representation . IEEE Trans. Sig. Proc. 42 , 3084 - 3091 ( 1994 )

11. Ozaktas , H.M. , Zalevsky , Z. , Kutay M.A. : The fractional Fourier transform . Wiley, Chichester, ( 2001 )

12. Rundblad-Labunets , E. , Labunets , V. , Astola , J. , Egiazarian , K. : Fast fractional Fourier-Clifford transforms . Second International Workshop on Transforms and Filter Banks . Tampere, Finland, TICSP Series, 5 , 376 - 405 ( 1999 )

13. Rundblad , E. , Labunets , V. , Astola , J. , Egiazarian , K. , Polovnev , S. : Fast fractional unitary transforms . Proc. of Conf. Computer Science and Information Technologies , Yerevan, Armenia, 223 - 226 ( 1999 )

14. Creutzburg , R. , Rundblad , E. , Labunets , V. : Fast algorithms for fractional Fourier transforms . Proc. of IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing . Antalya, Turkey, June 20-23, 383 - 387 ( 1999 )

15. Rundblad , E. , Labunets , V. , Astola , J. , Egiazarian , K. , Smaga , A. : Fast fractional Fourier and Hartley transforms . Proc. of the 1999 Finnish Signal Processing Symposium . Oulu, Finland, 291 - 297 ( 1999 )

16. Labunets , E. , Labunets , V. : Fast fractional Fourier transforms . Proc. of Eusipco-98 , Rhodes, Greece, 8 - 11 Sept. 1757 - 1760 ( 1998 )

17. Goldstein , J.: Semigroups of Linear Operators and Applications , Oxford U. Press, New York, ( 1985 )

18. Griffiths , R.B.: Consistent Quantum Theory . http://www.cambridge.org

19. Candan C. : The Discrete Fractional Fourier Transform , M.S. thesis, Bilkent Univ., Ankara, Turkey, ( 1998 )