=Paper=
{{Paper
|id=Vol-1452/paper3
|storemode=property
|title=Fast Infinitesimal Fourier Transform for Signal and Image Processing via Multiparametric and Fractional Fourier Transforms
|pdfUrl=https://ceur-ws.org/Vol-1452/paper3.pdf
|volume=Vol-1452
|dblpUrl=https://dblp.org/rec/conf/aist/OsthaimerLM15
}}
==Fast Infinitesimal Fourier Transform for Signal and Image Processing via Multiparametric and Fractional Fourier Transforms==
https://ceur-ws.org/Vol-1452/paper3.pdf
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
katya@capricat.com,
2
Ural Federal University, pr. Mira, 19,Yekaterinburg, 620002,
Russian Federation
vlabunets05@yahoo.com,
3
SPA Automatics, named after Academician N.A. Semikhatov,
Mamina Sibiryaka, 145, Yekaterinburg, 620002, Russian Federation
stmart2608@gmail.com
Abstract. The fractional Fourier transforms (FrFTs) is one-parametric
family of unitary transformations {F α }2π
α=0 . FrFTs found a lot of applica-
tions in signal and image processing. The identical and classical Fourier
transformations are both the special cases of the FrFTs. They corre-
spond 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 dis-
crete 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 re-
quires M N (2 + log2 N ) multiplications and M N log2 N additions. This
paper develops a new numerical algorithm, which requires 2M N multi-
plications and 3M N additions and which is based on the infinitesimal
Fourier transform.
Keywords: Fast fractional Fourier transform, infinitesimal Fourier trans-
form, Schrödinger 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 decom-
position 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 trans-
2π
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 one-
a ⊕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 spe-
cial 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ödinger 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)|.|Ψ
0 1 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
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
o F 1 = UΛU−1 .
The families F ( 0
α ,...,αN−1 ) 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 )
F 4 = I. Hence, F (a0 ,...,aN −1 ) F (b0 ,...,bN −1 = F 4 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π
N N
0, 1, ..., N −1. So, the ranges of (α0 , . . . , αN −1 ) and α are (Z/2πZ) = [0, 2π] =
N
[−π, π] and Z/2πZ = [0, 2π] = [−π, π], respectively.
21
3 Canonical FrFT
The continuous Fourier transform is a unitary operator F that maps square-
integrable functions on square-integrable ones and is represented on these func-
tions 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) =
n 2
dn x2 2
(−1) ex dx ne , then for n ∈ N0 , functions Ψn (x) = √ n1 √ Hn (x)e−x /2 are
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 eigen-
function. 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,
∞
X X∞
2 2 e−jαn Hn (x)Hn (y)
K α (x, y) := e−jαn Ψn (x) Ψn (y) = e−(x +y ) √ =
n=0 n=0
2n n! π
( ) ( )
1 2xye−jα − e−2jα x2 + y 2 x2 + y 2
=√ √ · exp exp − ,
π 1 − e−2jα 1 − e−2jα 2
(5)
where K α (x, y) is the kernel of the FrFT. In the last step we used the Mehler
formula [19]
∞
( )
X e−jαn Hn (x)Hn (y) 1 2xye−jα − e−2jα x2 + y 2
√ =√ √ exp .
n=0
2n n! π π 1 − e−2jα 1 − e−2jα
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 Z
e− 2 ( 2 α̂−α)
π
j 2
=p f (x) exp x + y 2 cos α − 2xy dx.
2π |sin α| R 2 sin α
There exist several algorithms for fast calculation of spectrum of the frac-
tional Fourier transform F α (y). But all of them are based on the following trans-
form of the FrFT:
j Z h i
e− 2 ( 2 α̂−α) ejy 2 sin α
π 2 cos α
x2
F α (y) = (F α f ) (y) = p f (x)ej 2 cot α e−jxy dx =
2π |sin α| R
= Aα (y) · F {f (x) · Bα (x)} (y),
j
−
e 2 ( π α̂−α
2 ) ejy 2 sin α
2 cos α
x 2
where Aα (y) = √ , Bα (x) = ej 2 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 addi-
tions).
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(α)
A(0) = , U(dα) = I + dU(0) = I + A(0)dα.
∂α α=0
Obviously,
∂U(α)
U(α0 + dα) = U(α0 ) + dU(α0 ) = U(α0 ) + dα =
∂α α0
= U(α0 ) + A(α0 )dα.
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
Define now the linear operator H = 21 dxd 2
2 − x + 1 . It is known that
1 d2 2
HΨn (x) = − x + 1 Ψn (x) = nΨn (x). (6)
2 dx2
From (4) and (6) we have
X∞ Z
∂ α ∂ α
j F (y) =j {F F } (y) = nΨn (y) Ψn (x)f (x)dx,
∂α α=0 ∂α α=0 n=0 R
∞
X Z
α
HF (x) = nΨn (y) Ψn (x)f (x)dx.
n=0 R
α α
Therefore, j ∂F∂α(x) = HF α (y), ∂F (x)
F α (x) = −jH∂α. The solutionh of this equa-
2 i
−jα 1 d 2
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
dα 1 d2 2
F = (I − jHdα) = I − j − x + 1 dα (7)
2 dx2
24
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 d2
F −1 j dy F. Obviously, x2 = Mx2 = −F 1 dy 2 F. Then
1 d2 d2
F dα = I − j + F −1 F + 1 dα.
2 dx2 dy 2
Discretization of x-domain with the interval discretization ∆x is equal to the
periodization of y-domain
2 2 2
d2 −1 d d −1 d
+F F + 1 −→ D∆x +F P2π/∆x F + 1.
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 −1 d
D∆x 2
+F P2π/∆x F + 1 −→
dx dy 2
2 2
d −1 d
−→ P2π/∆y D∆x + F P D
2π/∆x ∆y F + 1.
dx2 dy 2
An approximation for the second derivative can be given by the second order
central difference operator
d2 d2
f (x) ≈ f (n 1) − 2f (n) + F (n⊕1), F (y) ≈ F (k 1) − 2F (k) + F (k⊕1),
dx2 N N dy 2 N N
where N = 2π/∆x∆y. On the other hand,
2
d
F −1 F (y) F ≈ F −1
F (k 1) − 2F (k) + F (k⊕ 1) F=
dy 2 N N
−j 2π n j 2π n 2π
= f (n)e N − 2f (n) + f (n)e N = 2f (n) cos n−1 .
N
2
These allow one to give the approximation for H = 21 dx d 2
2 − x + 1 as follows:
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
2π 1
= − cos n − 3/2 f (n) + f (n 1) + f (n⊕1) .
N 2 N N
25
In the N -diagonal basis we have
f (0)
f (1)
f (2)
dα
F f (x) ≈ f (3) + j∆α×
..
.
f (N − 1)
−1/2 1/2 . . . 1/2 f (0)
1/2 cos(1Ω) − 3/2 1/2 . . .
f (1)
. 1/2 cos(2Ω) − 3/2 1/2 . . f (2)
× . cos(3Ω) − 3/2 1/2 . ,
. 1/2 f (3)
.. ..
. . . 1/2 . 1/2 .
1/2 . . . 1/2 −1/2 f (N − 1)
(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) . (9)
N 2 N N
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
d2
in the classical algorithm. In (8), we used O(h2 ) approximation 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 multi-
plications 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 pre-
sented 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-07-
12168, 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 func-
tional transforms. Proc. Nat. Acad. Sci. 23, 158–164 (1937)
4. Kober, H.: Wurzeln aus der Hankel- und Fourier und anderen stetigen Transforma-
tionen. 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 Fil-
ter 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 Sym-
posium. 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