Fast Infinitesimal Fourier Transform for Signal and Image Processing via Multiparametric and Fractional Fourier Transforms Ekaterina Ostheimer1 , Valeriy Labunets2 , and Stepan Martyugin3 1

Capricat LLC, 1340 S., Ocean Blvd., Suite 209, Pompano Beach, 33062 Florida, USA [email protected], 2 Ural Federal University, pr. Mira, 19,Yekaterinburg, 620002, Russian Federation [email protected], 3 SPA Automatics, named after Academician N.A. Semikhatov, Mamina Sibiryaka, 145, Yekaterinburg, 620002, Russian Federation [email protected]

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 }M i=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

1

Introduction 4

The idea of fractional powers of the Fourier operator {F a }a=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 4 {F a }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 trans2π form order. There is the angle parameterization {F α }α=0 , where α = πa/2 is a new angle parameter. Since this family depends on a single parameter,

19

4

2π

the fractional operators {F a }a=0 (or {F α }α=0 ) form the Fourier-Hermite onea ⊕b parameter strongly continuous unitary multiplicative group F a F b = F 4 (or α⊕β F α F β = F 2π ), where a⊕b = (a + b) mod4 (or α⊕ β = (α + β) mod2π) and 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 M set of spectra {F αi }i=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 N −1

Let F = [Fk (i)]k,i=0 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 = Ψ (i)|Ψ (i)|.|Ψ (i) be the matrix of the F-transform eigen-vectors. 0

1

N −1

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  a −1
F (a0 ,...,aN −1 ) := U diag λa0 0 , . . . , λNN−1 U−1 . (1)

20

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 λa0 , . . . , λaN −1 U−1 = UΛa U−1 . (2)

The zero-th-order fractional F-transform is equal to the identity transform: F 0 = UΛ0 U−1 = UU−1 = I , and the first-order fractional Fourier transform operator F 1 = F is equal n to the initial F-transform F 1 = UΛU−1 . o α ,...,α N−1 ) The families F ( 0 and {F a }a∈R form multi- and (α0 ,...,αN−1 )∈RN 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 a F b = F a+b . Indeed, F a F b = UΛa U−1 · UΛb U−1 = UΛa+b U−1 = F a+b and

F (a0 ,...,aN −1 ) F (b0 ,...,bN −1 ) = n  o  a −1
b −1 = U diag λa0 0 , . . . , λNN−1 U−1 · U diag λb00 , . . . , λNN−1 U−1 = n  o a −1 +bN −1 = U diag λ0a0 +b0 , . . . , λNN−1 U−1 = F (a0 +b0 ,...,aN −1 +bN −1 ) . N −1

Let F = [Fk (i)]k,i=0 be a discrete Fourier (N × N )-transform (DFT), then √ N −1 λn = ejπn/2 ∈ {±1, ±j} , where j = −1 and {Ψn (t)}n=0 are the Kravchuk polynomials. Definition 2. The multi-parametric and fractional DFT are n  o F (a0 ,...,aN −1 ) := U diag ejπ0a0 /2 , ejπ1a1 /2 , . . . , ejπ(N −1)aN −1 /2 U−1 , n  o F a := U diag ejπna/2 U−1 and

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

in a- and α-parameterizations, respectively, where α = πa/2. The parameters (a0 , . . . , aN −1 ) and a can be any real values. However, the operators F (a0 ,...,aN −1 ) and F a are periodic in each parameter with period 4 since (a0 ⊕b0 ,...,aN −1 ⊕bN −1 ) 4 F 4 = I. Hence, F (a0 ,...,aN −1 ) F (b0 ,...,bN −1 = F 4 and F a F b = a⊕b F 4 , where ai ⊕bi = (ai + bi ) mod4, ∀i = 0, 1, ..., N − 1. Therefore, the ranges 4

N

N

N

of (a0 , . . . , aN −1 ) and a are (Z/4Z) = [0, 4] = [−2, 2] 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 [−π, π] and Z/2πZ = [0, 2π] = [−π, π], respectively.

21

N

N

= [0, 2π]

=

3

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 Z 1 F (y) = (Ff ) (y) = √ f (x)e−jyx dx. (3) 2π x∈R
Relevant properties are that the square
F 2 f (x) = f (−x) is the inversion operator, and that its fourth power F 4 f (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) = 2 2 n dn x2 , then for n ∈ N0 , functions Ψn (x) = √ n1 √ Hn (x)e−x /2 are (−1) ex dx ne 2 n! π

the eigen-functions of the Fourier transform Z +∞ π 1 F [Ψn (x)] = Ψn (x) e2πjyx dx = λn Ψn (y) = e−j 2 n Ψn (y) 2π −∞ π

with λn = j n = 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 ∞ X 
K α (x, y) := U diag e−jαn U−1 = e−jαn Ψn (x) Ψn (y) .

(4)

n=0

Hence, K α (x, y) :=

∞ X

e−jαn Ψn (x) Ψn (y) = e−(x

n=0

1 =√ √ · exp π 1 − e−2jα

(

2

+y 2 )

∞ X e−jαn Hn (x)Hn (y) √ = 2n n! π n=0

2xye−jα − e−2jα x2 + y 2 1 − e−2jα


)

(

x2 + y 2 exp − 2


)

,

(5) where K α (x, y) is the kernel of the FrFT. In the last step we used the Mehler formula [19] (
) ∞ X 2xye−jα − e−2jα x2 + y 2 1 e−jαn Hn (x)Hn (y) √ =√ √ exp . 1 − e−2jα 2n n! π π 1 − e−2jα n=0 Expression (5) can be rewritten as r    1 − j cot α j  2 α 2 K (x, y) = exp (x + y ) cos α − 2xy , 2π 2 sin α

22

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 Z α α F (y) = (F f ) (y) = f (x)K α (x, y)dx = x∈R

j

ˆ ) e− 2 ( 2 α−α =p 2π |sin α| π

Z

f (x) exp

R




 j  2 x + y 2 cos α − 2xy dx. 2 sin α

There exist several algorithms for fast calculation of spectrum of the fractional Fourier transform F α (y). But all of them are based on the following transform of the FrFT: j

F α (y) = (F α f ) (y) =

ˆ ) ejy 2 sin α e− 2 ( 2 α−α p 2π |sin α| π

2 cos α

Z h

f (x)ej

R

x2 2

cot α

i

e−jxy dx =

= Aα (y) · F {f (x) · Bα (x)} (y), where Aα (y) =

j − e 2

(

) ejy 2 sin α 2π|sin α|

π α−α ˆ 2

√

2 cos α

, Bα (x) = ej

x2 2

cot α

.

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 log2 N additions and M N (2 + log2 N ) 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).

23

Definition 3. The infinitesimal generator A(0) of the group {U(α)}α∈R and the infinitesimal transform U(dα) are defined as follows [18,19]: ∂U(α) , U(dα) = I + dU(0) = I + A(0)dα. A(0) = ∂α α=0 Obviously,

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

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

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

= U(α0 ) + A(0)U(α0 )dα = [I + A(0)] U(α0 )dα. h i Hence, A(α0 ) = A(0)U(α0 ) and F α0 +dα (y) = I + A(0) F α0 (y)dα.  2  d 2 Define now the linear operator H = 21 dx 2 − x + 1 . It is known that HΨn (x) =

1 2



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

(6)

From (4) and (6) we have Z ∞ X ∂ α ∂ α j =j = nΨn (y) Ψn (x)f (x)dx, F (y) {F F } (y) ∂α ∂α R α=0 α=0 n=0 α

HF (x) = α

∞ X

nΨn (y)

Z

Ψn (x)f (x)dx.

R

n=0 α

(x) Therefore, j ∂F∂α(x) = HF α (y), ∂F this equaF α (x) = −jH∂α. The solutionh of  2 i 2 1 d  −jα 2 dx2 −x +1 tion is given by F α (x) = e−jαH F and F α = e−jαH = e . Obviously, F α+dα = F dα F α ' (I + dF α ) exp [−jαH] =   ∂F α = I+ dα exp (−jαH) = (I − jHdα) exp (−jαH) , ∂α

where the operator F

dα

1 = (I − jHdα) = I − j 2

24



 d2 2 − x + 1 dα dx2

(7)

is called the infinitesimal Fourier transform or the generator of the fractional Fourier transforms [17,18]. Let us introduce operators (Mx f ) (x) := xf (x) and (My F ) (y) := yF (y). Using the Fourier transform (3), the first of ones may be written as Mx =  d F. Obviously, x2 = Mx2 = −F 1 F −1 j dy

F dα = I − j

1 2



d2 + F −1 dx2

d2 dy 2



d2 dy 2

F. Then 

 F + 1 dα.

Discretization of x-domain with the interval discretization ∆x is equal to the periodization of y-domain  2   2    2  d2 d d d −1 −1 +F F + 1 −→ D∆x +F F + 1. P2π/∆x dx2 dy 2 dx2 dy 2 Discretization of y-domain with the interval discretization ∆y is equal to the periodization of x-domain  2    2  d d −1 D∆x +F P2π/∆x F + 1 −→ 2 dx dy 2   2   2  d d −1 + F P D F + 1. −→ P2π/∆y D∆x 2π/∆x ∆y dx2 dy 2 An approximation for the second derivative can be given by the second order central difference operator d2 f (x) ≈ f (n 1) − 2f (n) + F (n⊕1), dx2 N N

d2 F (y) ≈ F (k 1) − 2F (k) + F (k⊕1), dy 2 N N

where N = 2π/∆x∆y. On the other hand,  2    d −1 F −1 F (y) F ≈ F F (k 1) − 2F (k) + F (k⊕ 1) F= dy 2 N N     2π −j 2π n j 2π n N N = f (n)e − 2f (n) + f (n)e n−1 . = 2f (n) cos N  2  d 2 These allow one to give the approximation for H = 21 dx as follows: 2 − x + 1   2  1 d 2 Hf (x) = − x + 1 f (x) ≈ 2 dx2      1 2π f (n 1) − 2f (n) + f (n⊕1) + 2f (n) cos n − 1 + f (n) = ≈ 2 N N N     1 2π = − cos n − 3/2 f (n) + f (n 1) + f (n⊕1) . N 2 N N

25

In the N -diagonal basis we have 

    dα F f (x) ≈    



f (0) f (1) f (2) f (3) .. . f (N − 1)



     + j∆α×   

  −1/2 1/2 . . . 1/2 f (0)  1/2 cos(1Ω) − 3/2   1/2 . . .    f (1)   .  f (2)  1/2 cos(2Ω) − 3/2 1/2 . .    × . ,  . 1/2 cos(3Ω) − 3/2 1/2 .    f (3)      .. ..  .  . 1/2  . . . 1/2 f (N − 1) 1/2 . . . 1/2 −1/2 (8) 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∆α×     2π 1 × cos k − 3/2 F αi (k) + F αi (k 1) + F αi (k⊕ + 1) . N 2 N N

(9)

It is easy to see that this algorithm requires only 2M N multiplications and 3M N additions vs. M N (2 + log2 N ) multiplications and M N log 2 N additions  in the classical algorithm. In (8), we used O(h2 ) approximation

d2 dx2 f

(k) ≈

(f (k − 1) − 2f (k) + f (k + 1)) . More fine approximations O(h2k ) also can be used [19].

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.

26

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)

27