*199*
*41*
*2MB*

*English*
*Pages (160 Seiten)) : color illlustration, Fotografien
[160]*
*Year 2014*

- Author / Uploaded
- Ell
- Todd A.; Le Bihan
- Nicolas; Sangwine
- Stephen J

*Table of contents : Cover......Page 3Title Page......Page 5Copyright......Page 6Contents......Page 7Nomenclature......Page 11Preface......Page 13Introduction......Page 15 1.1. Definitions......Page 25 1.2. Properties......Page 26 1.3.1. Exponential of a pure quaternion......Page 31 1.3.2. Exponential of a full quaternion......Page 33 1.3.3. Logarithm of a quaternion......Page 34 1.4.1. Polar forms......Page 35 1.4.2. The Cj-pair notation......Page 39 1.4.3. R and C matrix representations......Page 41 1.6. Subfields......Page 42 2.1. Euclidean geometry (3D and 4D)......Page 45 2.1.2. 3D rotations......Page 46 2.1.4. 3D dilations......Page 48 2.1.6. 4D rotations......Page 49 2.2. Spherical geometry......Page 50 2.3. Projective space (3D)......Page 52 2.3.1. Systems of linear quaternion functions......Page 55 2.3.2. Projective transformations......Page 57Chapter 3. Quaternion Fourier Transforms......Page 59 3.1.1. Definitions......Page 62 3.1.2. Basic transform pairs......Page 64 3.1.3. Decompositions......Page 66 3.1.4. Inter-relationships between definitions......Page 69 3.1.5. Convolution and correlation theorems......Page 71 3.2.1. Definitions......Page 72 3.2.2. Basic transform pairs......Page 76 3.2.3. Decompositions......Page 78 3.2.4. Inter-relationships between definitions......Page 79 3.3.1. Coding......Page 81 3.3.3. Verification of transforms......Page 86 4.1.1. Classical grayscale image convolution filters......Page 91 4.1.3. Quaternion convolution......Page 94 4.1.4. Quaternion image spectrum......Page 97 4.2. Generalized correlation......Page 103 4.2.1. Classical correlation and phase correlation......Page 105 4.2.2. Quaternion correlation......Page 110 4.2.3. Quaternion phase correlation......Page 112 4.3.1. Important properties of 1D QFT of a complex signal z(t)......Page 115 4.3.2. Hilbert transform and right-sided quaternion spectrum......Page 120 4.3.3. The quaternion signal associated with a complex signal......Page 122 4.3.4. Instantaneous amplitude and phase......Page 125 4.3.5. The instantaneous frequency of a complex signal......Page 126 4.3.6. Examples......Page 128 4.3.7. The quaternion Wigner-Ville distribution of a complex signal......Page 133 4.3.9. The mean frequency formula......Page 137Bibliography......Page 141Index......Page 147Supplemental Images......Page 153*

W478-Sangwine.qxp_Layout 1 10/04/2014 16:25 Page 1

FOCUS SERIES in DIGITAL SIGNAL AND IMAGE PROCESSING

The fourth and final chapter is dedicated to the illustration of the use of QFT to process color images and complex improper signals. The concepts presented in this chapter are illustrated on simulated and real images and signals.

Todd A. Ell is an Engineering Fellow at UTC Aerospace Systems, Burnsville, MN, USA. He is also a Visiting Fellow at the University of Essex, Colchester, UK. His interests include the study and application of hypercomplex algebras to dynamic systems analysis. Nicolas Le Bihan is currently a Chargé de Recherche at the CNRS working in the Department of Images & Signals (DIS) at the GIPSA-Lab in Grenoble, France. His research interests include polarized signal processing, statistical signal processing on groups and manifolds, geometric (Berry) phases, random processes on non-commutative algebraic structures and applications in physics and geophysics. Stephen J. Sangwine is a Senior Lecturer with the School of Computer Science and Electronic Engineering, University of Essex, Colchester, UK. His interests include linear vector filtering and transforms of vector signals and images, color image processing, and digital hardware design.

www.iste.co.uk

Z(7ib8e8-CBEHIB(

Quaternion Fourier Transforms for Signal and Image Processing

Following the Introduction, Chapter 1 introduces the quaternion algebra H and presents some properties which will be of use in the subsequent chapters. Chapter 2 gives an overview of the geometric transformations which can be represented using quaternions. Chapter 3 provides the definition and properties of QFT. The signals and images considered are those with vector-valued samples/pixels.

Todd A. Ell Nicolas Le Bihan Stephen J. Sangwine

This book presents the state of the art, together with the most recent research results, in the use of Quaternion Fourier Transforms (QFT) for the processing of color images and complex valued signals. It is based on the work of the authors in this area since the 1990s and presents the mathematical concepts, computational issues and applications on images and signals. The book, together with the MATLAB toolbox developed by two of the authors (QTFM, http://qtfm.sourceforge.net/), allows the reader to make use of the presented concepts and experiment with them in practice through the examples provided in the book.

FOCUS DIGITAL SIGNAL AND IMAGE PROCESSING SERIES

Quaternion Fourier Transforms for Signal and Image Processing Todd A. Ell, Nicolas Le Bihan and Stephen J. Sangwine

Quaternion Fourier Transforms for Signal and Image Processing

FOCUS SERIES Series Editor Francis Castanié

Quaternion Fourier Transforms for Signal and Image Processing

Todd A. Ell Nicolas Le Bihan Stephen J. Sangwine

First published 2014 in Great Britain and the United States by ISTE Ltd and John Wiley & Sons, Inc.

Apart from any fair dealing for the purposes of research or private study, or criticism or review, as permitted under the Copyright, Designs and Patents Act 1988, this publication may only be reproduced, stored or transmitted, in any form or by any means, with the prior permission in writing of the publishers, or in the case of reprographic reproduction in accordance with the terms and licenses issued by the CLA. Enquiries concerning reproduction outside these terms should be sent to the publishers at the undermentioned address: ISTE Ltd 27-37 St George’s Road London SW19 4EU UK

John Wiley & Sons, Inc. 111 River Street Hoboken, NJ 07030 USA

www.iste.co.uk

www.wiley.com

© ISTE Ltd 2014 The rights of Todd A. Ell, Nicolas Le Bihan and Stephen J. Sangwine to be identified as the author of this work have been asserted by them in accordance with the Copyright, Designs and Patents Act 1988. Library of Congress Control Number: 2014934161 British Library Cataloguing-in-Publication Data A CIP record for this book is available from the British Library ISSN 2051-2481 (Print) ISSN 2051-249X (Online) ISBN 978-1-84821-478-1

Printed and bound in Great Britain by CPI Group (UK) Ltd., Croydon, Surrey CR0 4YY

Contents

N OMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

P REFACE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

I NTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xiii

C HAPTER 1. Q UATERNION A LGEBRA

. . . . . . . . . . . . . . . . . . . . .

1.1. Deﬁnitions . . . . . . . . . . . . . . . . . . 1.2. Properties . . . . . . . . . . . . . . . . . . . 1.3. Exponential and logarithm of a quaternion 1.3.1. Exponential of a pure quaternion . . . . 1.3.2. Exponential of a full quaternion . . . . 1.3.3. Logarithm of a quaternion . . . . . . . . 1.4. Representations . . . . . . . . . . . . . . . . 1.4.1. Polar forms . . . . . . . . . . . . . . . . 1.4.2. The Cj -pair notation . . . . . . . . . . . 1.4.3. R and C matrix representations . . . . . 1.5. Powers of a quaternion . . . . . . . . . . . . 1.6. Subﬁelds . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . .

1 2 7 7 9 10 11 11 15 17 18 18

C HAPTER 2. G EOMETRIC A PPLICATIONS . . . . . . . . . . . . . . . . . . .

21

2.1. Euclidean geometry (3D and 4D) 2.1.1. 3D reﬂections . . . . . . . . . 2.1.2. 3D rotations . . . . . . . . . . 2.1.3. 3D shears . . . . . . . . . . . 2.1.4. 3D dilations . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

. . . . .

. . . . . . . . . . . .

1

. . . . .

. . . . .

21 22 22 24 24

vi

Quaternion Fourier Transforms for Signal and Image Processing

2.1.5. 4D reﬂections . . . . . . . . . . . . . . . 2.1.6. 4D rotations . . . . . . . . . . . . . . . . 2.2. Spherical geometry . . . . . . . . . . . . . . 2.3. Projective space (3D) . . . . . . . . . . . . 2.3.1. Systems of linear quaternion functions 2.3.2. Projective transformations . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

25 25 26 28 31 33

. . . . . . . . . . . .

35

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

38 38 40 42 45 47 48 48 52 54 55 57 57 62 62

. . . . . . . . . . . . . . .

67

C HAPTER 3. Q UATERNION F OURIER T RANSFORMS 3.1. 1D quaternion Fourier transforms . . . . . 3.1.1. Deﬁnitions . . . . . . . . . . . . . . . . 3.1.2. Basic transform pairs . . . . . . . . . . 3.1.3. Decompositions . . . . . . . . . . . . . 3.1.4. Inter-relationships between deﬁnitions . 3.1.5. Convolution and correlation theorems . 3.2. 2D quaternion Fourier transforms . . . . . 3.2.1. Deﬁnitions . . . . . . . . . . . . . . . . 3.2.2. Basic transform pairs . . . . . . . . . . 3.2.3. Decompositions . . . . . . . . . . . . . 3.2.4. Inter-relationships between deﬁnitions . 3.3. Computational aspects . . . . . . . . . . . . 3.3.1. Coding . . . . . . . . . . . . . . . . . . 3.3.2. Veriﬁcation . . . . . . . . . . . . . . . . 3.3.3. Veriﬁcation of transforms . . . . . . . . C HAPTER 4. S IGNAL AND I MAGE P ROCESSING

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

4.1. Generalized convolution . . . . . . . . . . . . . . . . . . . . . 4.1.1. Classical grayscale image convolution ﬁlters . . . . . . . 4.1.2. Color images as quaternion arrays . . . . . . . . . . . . . 4.1.3. Quaternion convolution . . . . . . . . . . . . . . . . . . . 4.1.4. Quaternion image spectrum . . . . . . . . . . . . . . . . . 4.2. Generalized correlation . . . . . . . . . . . . . . . . . . . . . 4.2.1. Classical correlation and phase correlation . . . . . . . . 4.2.2. Quaternion correlation . . . . . . . . . . . . . . . . . . . . 4.2.3. Quaternion phase correlation . . . . . . . . . . . . . . . . 4.3. Instantaneous phase and amplitude of complex signals . . . 4.3.1. Important properties of 1D QFT of a complex signal z(t) 4.3.2. Hilbert transform and right-sided quaternion spectrum . 4.3.3. The quaternion signal associated with a complex signal . 4.3.4. Instantaneous amplitude and phase . . . . . . . . . . . . . 4.3.5. The instantaneous frequency of a complex signal . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. 67 . 67 . 70 . 70 . 73 . 79 . 81 . 86 . 88 . 91 . 91 . 96 . 98 . 101 . 102

Contents

4.3.6. Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.7. The quaternion Wigner-Ville distribution of a complex signal 4.3.8. Time marginal . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.9. The mean frequency formula . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

vii

104 109 113 113

B IBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 I NDEX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Nomenclature

Symbol R C H V(H) Cμ I (.) (.) i, j, k μ i j k

MR (.) MC (.) [[.]] [.] S(.) V(.) (q1 , q2 ) AB q z qμ q −1 p, q μ×η μ⊥η

Meaning Set of real numbers Set of complex numbers Set of quaternions Set of pure quaternions Subﬁeld of H Complex root of −1 Real part Imaginary part Quaternion basis elements Pure quaternion i-imaginary part j-imaginary part k-imaginary part Real 4 × 4 matrix representation Complex 2 × 2 matrix representation Alternate real 4 × 4 matrix representation Real 4 × 1 vector representation Scalar part Vector part Cj -pair notation of q Arc of great circle between A and B Quaternion conjugate Complex conjugate Involution Inverse Inner product Cross product Orthogonality

Section section 1.1 section 1.1 section 1.1 section 1.1 section 1.6 (See note below.) section 1.1 section 1.1 section 1.1 section 1.1 section 1.1 section 1.1 section 1.1 section 1.4 section 1.4 section 2.3 section 2.3 section 1.1 section 1.1 section 1.4 section 2.2 section 1.2 section 1.4 section 1.2 section 1.2 section 1.2 section 1.2 section 1.2

Equation

[1.2] [1.2] [1.4] [1.4] [1.4] [1.69] [1.70] [2.20] [2.20] [1.3] [1.3] [1.18] [1.54] [1.23] [1.26] [1.10] [1.9]

x

Quaternion Fourier Transforms for Signal and Image Processing

q p {1, i, j, k} q |q| q˜ f ∗g f g F {.} Pη [.] Rq [.] Sα,β,μ [.] Dμ,α [.] L1 (G; K)

Bicomplex product Canonical basis in H Norm Modulus Unit quaternion |˜ q| = 1 Convolution Correlation Fourier transform Reﬂection operator Rotation operator Shear operator Dilation operator Space of absolutely integrable K-valued functions taking arguments in G L2 (G; K) Space of square integrable K-valued functions taking arguments in G sgn(.) Signum function p.v.(.) Principal value of an integral

section 1.4 section 1.1 section 1.2 section 1.2 section 1.2 section 3 section 3 section 3 section 2.1 section 2.1 section 2.1 section 2.1

[1.65] [1.14] [1.16] [1.30] [3.1] [3.1] [3.1] [2.1] [2.3] [2.9] [2.10]

section 3.1 section 4.3 section 4.3 section 4.3

[4.31] [4.32]

N OTE.– The complex root of −1 which is usually denoted i, or, in engineering texts j, is denoted throughout this book by a capital letter I, in order to avoid any confusion with the ﬁrst of the three quaternion roots of −1, all three of which are denoted throughout in bold font like this: i, j, k.

Preface

This book aims to present the state of the art, together with the most recent research results in the use of quaternion Fourier transforms (QFTs) for the processing of color images and complex-valued signals. It is based on the work of the authors in this area since the 1990s and presents the mathematical concepts, computational issues and some applications to signals and images. The book, together with the MATLAB® toolbox developed by the authors, [SAN 13b] allows the readers to make use of the presented concepts and experiment with them in practice through the examples provided. Todd A. E LL Minneapolis, MN, USA Nicolas le B IHAN Melbourne, Australia Stephen J. S ANGWINE Colchester, UK April 2014

Introduction

This book covers a topic that combines two branches of mathematical theory to provide practical tools for the analysis and processing of signals (or images) with three- or four-dimensional samples (or pixels). The two branches of mathematics are not recent developments, but their combination has occurred only within the last 25–30 years, and mostly since just before the millennium. I.1. Fourier analysis Fourier analysis was, in 1822, with Joseph Fourier’s development of techniques, the ﬁrst to analyze mathematical functions into sinusoidal components. In signal and image processing, Fourier’s ideas underpin the two fundamental representations of a signal: one in the time (or image) domain where the signal (or image) is represented by samples (or pixels) with amplitudes and the other in the frequency domain where the signal (or image) is represented by sinusoidal frequency components, each with an amplitude and a phase. Mathematically, these concepts are not limited to time and frequency: one can use Fourier analysis on a function of any variable, resulting in a representation in terms of sinusoidal functions of that variable. However, this book is concerned with signal and image processing, and we will therefore use the terms time and frequency rather than more general concepts. It should be understood throughout that when we talk of images, the concept of time is replaced by the two spatial coordinates that deﬁne pixel position within an image. Today, Fourier analysis is classically taught to mathematicians, scientists and engineers in several related ways, each applicable to a speciﬁc subset of mathematical functions or signals: 1) Fourier series analysis [SNE 61] in which continuous periodic functions of time, with inﬁnite duration, are represented as sums of cosine and sine functions, each with inﬁnite duration;

xiv

Quaternion Fourier Transforms for Signal and Image Processing

– Fourier integrals or transforms [BRA 00, ROB 68] in which continuous (but aperiodic) functions of time are represented as continuous functions of frequency (or vice versa); 2) Discrete Fourier transforms in which signals deﬁned at discrete intervals in time are represented in the frequency domain by cosine and sine functions. This topic is broken down into: – discrete-time Fourier transforms, in which discrete-time signals of limited duration are represented as continuous frequency-domain distributions; – discrete Fourier transforms, in which discrete-time, discretized (that is digital) signals of ﬁnite duration are represented by a ﬁnite-length array of digital frequency coefﬁcients. (These are usually computed numerically using the fast Fourier transform (FFT)). The key to all of the above ideas is the representation of a signal using complex exponentials, often known as harmonic analysis, although this term has a somewhat wider meaning in mathematics than its usage in signal and image processing. The complex exponential with angular frequency ω and phase φ: f (t) = A exp(ωt + φ) = A (cos(ωt + φ) + I sin(ωt + φ)) has cosine and sine components in its real and imaginary parts, respectively. Since, in this book, we are concerned with signals that have three- or four-dimensional samples, it is helpful to consider classical Fourier analysis in terms of complex exponentials rather than in terms of separate cosines and sines. Figure I.1 shows a real-valued signal (on the left-hand side of the plot, with time increasing away from the viewer). The signal is a sawtooth waveform reconstructed from its ﬁrst ﬁve non-zero harmonics, which are plotted in the center of the ﬁgure as helices. (The horizontal spacing between the helices is introduced simply to make them clearer: there is no mathematical signiﬁcance to it). The ﬁve helices on the left are the positive frequency complex exponentials and the ﬁve helices on the right are the negative frequencies. Note that the positive and negative frequency exponentials have opposite directions of rotation. The real parts of the harmonics are projected onto the right-hand side of the ﬁgure (these sum to give the reconstructed waveform on the left) and the imaginary parts of the harmonics are projected onto the base of the ﬁgure (these cancel out because the exponentials occur in complex conjugate pairs at positive and negative frequencies, a symmetry due to the original signal being real-valued). In general, with a complex signal analyzed into complex exponentials in the same way, there would be no symmetry between the positive and negative frequency exponentials. This case is a useful model for what follows in this book, where we consider signals and images with three- and four-dimensional samples. Figure I.2 shows a complex signal constructed by bandlimiting a random complex signal.

Introduction

xv

Time is plotted on the right, increasing to the right, and at each time instant the signal has a complex value. The signal evolution over time traces out a path in the complex plane, and the ﬁgure renders this path as a three-dimensional view by plotting the signal values, in effect, on a stack of 2,000 transparent complex planes perpendicular to the time axis. The real and imaginary parts of the signal are also plotted on the base of the axes, and on the rear plane of the axes. Analysis of a complex signal into positive and negative frequency complex exponentials is not conceptually different from the real case depicted in Figure I.1: each complex exponential will have an amplitude and phase, and their sum will reconstruct the original signal. 300

200

Amplitude

100

0

−100

−200

−300 2000 1500 1000 500 Time

0

Figure I.1. Analysis of a real signal into complex exponential harmonics

The time and frequency domain representations of a signal are not mutually exclusive: the ﬁeld of time-frequency analysis [FLA 98] is concerned with intermediate representations that combine aspects of time and frequency. The need for intermediate representations arises due to the variation of frequency content in a signal over time. This is not an easy concept to understand, but it follows from the uncertainty principle or Gabor limit: a signal cannot be bandlimited (i.e. with frequency content limited to a ﬁnite range of frequencies) and simultaneously be of limited time duration. A pure sinusoidal signal with unlimited duration (inﬁnite

xvi

Quaternion Fourier Transforms for Signal and Image Processing

extent) can be represented in the frequency domain as an impulse (that is a function with zero value everywhere except at one frequency point). Conversely, an impulse in the time domain has inﬁnite bandwidth. However, a signal that contains a speciﬁc frequency for a limited time requires a time-frequency representation. Examples of such signals occur widely in the real world: speech and music contain frequencies that are present for a short time (one note played on a musical instrument, for example, which lasts for the duration of the note, plus some reverberation time afterward). An in-depth discussion of these ideas is outside the scope of this book, but is assumed to be understood; although much of the contents of the book relates to Fourier transforms, the quaternion approach can easily be applied to time-frequency concepts, such as fractional and short-time Fourier transforms, by combining quaternion transform formulations with existing knowledge from classical signal processing.

0.3 0.2

ℑ(f(t))

0.1 0 −0.1 −0.2 −0.3 −0.4 −0.2

2000

−0.1

1500

0

1000

0.1

500 0.2

ℜ(f(t))

0

t

Figure I.2. A bandlimited complex signal showing real and imaginary parts projected onto the base and rear of the grid box

I.2. Quaternions In this book, we are concerned with signals and images that have vector-valued samples (that is samples with three or more components), and their processing using

Introduction

xvii

Fourier transforms based on four-dimensional hypercomplex numbers (quaternions). In Chapter 4 (section 4.3), we show that quaternion Fourier transforms also have applications for the processing of complex signals, exploiting the symmetry properties of a quaternion Fourier transform that are missing from a complex Fourier transform. A vector-valued signal (in three dimensions, for example) evolves over time and traces out a path in three-dimensional space. To render a plot of such a signal requires four dimensions, and therefore we cannot produce a graphical representation like the one in Figure I.2. Decomposition of a vector-valued signal into harmonic components requires a Fourier transform in an algebra with dimension higher than 2, and this is the motivation for the use of quaternions, which, as we will see, are the next available higher-dimensional algebra after the complex numbers. Quaternions followed the work of Fourier just over 20 years later, Sir William Rowan Hamilton in 1843 to generalize the complex numbers to three dimensions, was forced to resort to four dimensions in order to obtain what we now call a normed division algebra, that is, an algebra where the norm of a product equals the product of the norms, and where every element of the algebra (except zero) has a multiplicative inverse [WAR 97]. Hamilton opened a door in mathematics to hypercomplex algebras in general [STI 10, Chapter 20], [KAN 89], leading to the octonions [CON 03, BAE 02] in less than a year, and the Clifford algebras about 30 years later [LOU 01, POR 95]. I.3. Quaternion Fourier transforms Quaternion Fourier transforms, the subject of this book, are a generalization of the classical Fourier transform to process signals or images with three- or four-dimensional samples. Such signals arise very naturally in the physical world from the three dimensions of physical space. Quite independently, for very different (physiological) reasons connected with the trichromatic nature of human color vision [MCI 98], color images have three components per pixel. The fourth dimension of the quaternions plays a role in at least two ways: the frequency-domain representation of a signal with three-dimensional samples requires four dimensions (just as in the complex case, two dimensions are required in the frequency domain, even if the original signal has one-dimensional samples). But more generally, the four dimensions of the quaternions can be used to represent a most general set of geometric operations in three dimensions using homogeneous coordinates, which are explored in a later chapter (see section 2.3) and in [SAN 13a]. Of course, generalizations to higher dimensions are possible, and there is a wide range of work on Clifford Fourier transforms, which is outside the scope of this book (we refer the readers to a recent volume for further details [HIT 13], and in particular the historical introduction contained within [BRA 13]).

xviii

Quaternion Fourier Transforms for Signal and Image Processing

I.4. Signal and image processing Fourier transforms are a fundamental tool in signal and image processing. They convert a signal or image from a representation based on sample or pixel amplitudes into a representation based on the amplitudes and phases of sinusoids. The latter representation is said to be in the frequency domain, and the original signal is said to be in the time domain for a signal which is a time series, or in the image domain for an image captured with a camera or scanner. Of course, signals may be encountered that are not time series, for example, measurements of some physical quantity made at (regular) intervals in space; in this book, we will use the terminology of time series for simplicity, since the processing of other signals is mathematically no different. The Fourier transformation is invertible, which means that the original signal or image may be recovered from the frequency domain representation. More interestingly, the frequency domain representation may be modiﬁed before inversion of the transform, so that the recovered signal or image is a modiﬁed version of the original, for example, with some frequencies or bands of frequencies suppressed, attenuated or ampliﬁed. In some applications, inversion of the transform is not needed: the processing performed in the frequency domain directly yields information that can be immediately utilized. An example is computer vision, where a decision based on analysis of an image may result in an action without any need to construct an image from the processed frequency domain representation. At a more detailed level, another example includes correlation, where the signal or image is processed in the frequency domain to yield information about the location of a known object within an image (the same applies in signal processing to ﬁnd a known signal occurring within a longer, noisy signal). The classical Fourier transform is inherently based on complex numbers. This is obvious from the fact that the frequency domain representation must represent both the amplitude and the phase of each frequency present in the signal or image. The symmetry of the transform means that the signal may be complex without any modiﬁcation of the transform. (There are some specialized variants of the Fourier transform that handle only real signals, for example the Hartley transform [BRA 86]). Given a signal with three components (representing, for example, acceleration in three mutually perpendicular directions), how can a frequency domain representation be calculated? The question is very similar if one considers a color image: is it possible to construct a holistic frequency-domain representation of the entire image? Obviously one can compute separate classical (i.e. complex) Fourier transforms of the three components in both of these cases, but one then has three separate frequency-domain representations, each representing one aspect of the original image (the frequency content of one of the color or luminance/chrominance channels). Processing of separate representations is sometimes known as marginal processing, for reasons connected with techniques in the hand computation of marginal distributions in statistics [TRU 53, section 1.22]. It is axiomatic in this book

Introduction

xix

that marginal processing is not the best way to handle signals and images with more than two components per sample, but we will attempt to justify this belief throughout the book, by showing how holistic approaches with quaternions yield better results. There is another reason for using a quaternion Fourier transform in some applications, and it provided the motivation for the earliest published work on quaternion Fourier transforms (in the ﬁeld of nuclear magnetic resonance (NMR)). When a two-dimensional signal is captured (that is samples are measured over a two-dimensional grid, like an image), it is sometimes necessary to regard the two dimensions of the sampling grid as independent time-like axes. Processing such a signal with a classical two-dimensional complex transform mixes the two dimensions, whereas a suitably formulated quaternion transform does not. This is because it is possible to associate each of the time-like dimensions with a different dimension of the four-dimensional quaternion space, thus keeping the frequency-domain representations of the ﬁrst and second time-like axes apart. There were two independent (as far as we are aware) early formulations of quaternion Fourier transforms, by Ernst [ERN 87, section 6.4.2] and Delsuc [DEL 88, equation 20], which are almost equivalent (they differ in the relative placement of the exponentials and the signal, and in the signs, inside the exponentials): F (ω1 , ω2 ) =

∞

∞

−∞

−∞

f (t1 , t2 )eiω1 t1 ejω2 t2 dt1 dt2

[I.1]

Note that the two-dimensional signal f (t1 , t2 ), is scalar-valued (i.e. not quaternion-valued as in most of the cases discussed in this book). The two time-like axes t1 and t2 are treated using separate quaternion roots of −1 (i and j), and therefore they are not mixed. This may also be regarded as being due to the orthogonality of the imaginary parts of the two exponentials. Similar considerations motivated Thomas Bülow [BÜL 99, BÜL 01] when he processed grayscale images using a quaternion Fourier transform. By using a transform with samples of dimension greater than 2, Bülow was able to study symmetries present in certain images in a way that is not possible with the two-dimensional complex Fourier transform. I.5. Other hypercomplex algebras Of course, there are alternatives to quaternions for the construction and computation of Fourier transforms for the type of applications suggested so far in this introductory chapter, and it is worth brieﬂy reviewing them here in order to provide a full picture. As already noted, the quaternion algebra is not the only hypercomplex algebra. In fact, the quaternion algebra is a speciﬁc example of a more general class of hypercomplex algebras discovered by William Kingdon Clifford in 1876 (about 30

xx

Quaternion Fourier Transforms for Signal and Image Processing

years after Hamilton’s discovery of quaternions) and known as Clifford algebras. The Clifford algebras include the complex numbers and quaternions, but not the octonions, curiously. However, as already mentioned, the quaternions share with the real numbers and the complex numbers a very speciﬁc property that sets them aside from all other Clifford algebras: every non-zero quaternion, q = 0, has a multiplicative inverse such that q q −1 = q −1 q = 1. Furthermore, the quaternion algebra is normed. This means that it is possible to deﬁne a norm (representing the squared length of the quaternion in four dimensions) such that the norm of a product of two quaternions equals the product of the norms of the two quaternions taken separately: pq = p q . This is discussed in section 1.2 (see [1.14]). Hypercomplex algebras in general have other troublesome properties. Many contain values that are idempotent or nilpotent. An idempotent value q squares to give itself: q 2 = q; and a nilpotent value squares to give zero: q 2 = 0. Such values obviously have the potential to cause problems in numerical algorithms [ALF 07], and the choice of the quaternion algebra avoids these problems entirely, because there are no nilpotent or idempotent quaternions other than 0 and 1, respectively. The one property of the quaternions that cannot be avoided is that multiplication of quaternions is not commutative. This means that pq gives a different result in general from qp for two arbitrary quaternions p and q. The reason for this can be stated quite simply – the vector (or cross) product in three dimensions is not commutative, and the product of two quaternions contains a vector product. It is important to understand that this is an inherent property of three-dimensional geometry and is not speciﬁc to the quaternions. This is again discussed in Chapter 1. Non-commutative multiplication can be avoided by choosing a different hypercomplex algebra, but since any hypercomplex algebra which is commutative contains divisors of zero (a consequence of the Frobenius [DIC 14, section 11] ), any attempt to avoid non-commutative multiplication will inevitably lead to other problems, which may well be more troublesome than non-commutativity. Non-commutative multiplication also occurs in linear algebra, of course, where the product of matrices is dependent on ordering; so it should not cause undue concern to anyone contemplating using quaternions. I.6. Practical application The ideas and concepts in this book are realisable in practice in several different ways, particularly using software. I.6.1. Software libraries The library [SAN 13b] permits experimentation with transforms and other algorithms operating on three- and four-dimensional data in MATLAB®. Since

Introduction

xxi

MATLAB® can generate C and C++ code (with some restrictions on supported language features), quaternion code can, in principle, be used to generate code for stand-alone applications, subject to licensing1. Alternatively, code can be custom-written, using a quaternion library such as the Boost library for C++, which contains some quaternion functions in the Math toolkit [BRI 13]. This is discussed again in section 3.3, particularly with respect to the use of decompositions into complex transforms to avoid the need to code elaborate algorithms directly in quaternion code. Both the QFTM and Boost libraries adopt the approach of directly coding quaternion operations, that is they represent quaternions as quadruplets of real (or complex) values, and provide elementary functions to add and multiply quaternions, implementing the famous rules for ijk given in section 1.1, directly in code. I.6.2. Matrix representations An alternative to the use of quaternion libraries is possible, using matrix representations, which we discuss here, and will return to with a very practical application (to veriﬁcation) in section 3.3.3.1. Hypercomplex algebras (with the exception of the octonions [CON 03, BAE 02], which are not associative) have matrix representations. What this means is that for a given algebra, there exists a matrix algebra with real or complex elements that is equivalent to the given hypercomplex algebra, in the sense that multiplication (and addition of course) of the matrix representations is equivalent to multiplication in the hypercomplex algebra. There are also other equivalences, for example the norm of a hypercomplex value may be equivalent to the determinant of the matrix representation. The matrix algebra using the given representation is said to be isomorphic to the hypercomplex algebra. Matrix representations for the quaternions are discussed in section 1.4.3, but we discuss the ramiﬁcations here. Given the existence of a matrix representation of quaternions, it is theoretically possible to substitute matrix representations for quaternions, both in algebraic manipulation and in computer coding (the same is true for other hypercomplex algebras except the octonions). Doing so can be a useful technique in theoretical development because it can reduce a hypercomplex problem to a problem involving real (or complex) matrices, and thus provide a deep insight into the relation between the hypercomplex case and the well understood real and complex cases. However, there are disadvantages of using a matrix representation compared to a direct quaternion approach, in practice: 1 QTFM is licensed under the GNU General Public License.

xxii

Quaternion Fourier Transforms for Signal and Image Processing

1) Computation with matrices is numerically inefﬁcient, and it requires four times as much memory as a direct quaternion representation storing only four values. A quaternion product requires 16 multiplications and 12 additions, whereas the equivalent 4 × 4 matrix product requires 64 multiplications and 48 additions, an increase by a factor of 4. This disadvantage effectively rules out the use of matrix representations for implementation, even for computer simulations (as run-time is four times less using a quaternion library that directly calculates quaternion products than using a general matrix package and a matrix representation for each quaternion). 2) The result of a sequence of arithmetic operations in matrix form may not be an accurate representation of a quaternion matrix. This disadvantage again rules out the use of matrix representations for implementation (a quaternion library will yield more accurate results). 3) The matrix representation provides little geometric insight. As will be shown in Chapter 2, the quaternion algebra provides a remarkably intuitive link with the geometry of three or four dimensions. It is possible to manipulate quaternion symbols algebraically in order to derive expressions for geometric operations. This is the central idea in geometric algebras [SOM 01]. Geometric algebras, as might be expected from the preceding text, are not division algebras; so there is a price to be paid for their additional geometric utility in the form of divisors of zero, which make them less attractive for applications in digital signal and image processing. The matrix representation is certainly useful, and it is helpful in any quaternion library to have the ability to convert between a direct quaternion representation and the matrix representation. In the QTFM library [SAN 13b], for example, functions called adjoint2 and unadjoint are provided to perform the conversion, even for matrices of quaternions (the adjoint in the latter case is a block matrix with each block representing one quaternion). I.7. Overview of the remaining chapters The rest of the book is divided into four chapters. Chapter 1 covers the quaternion algebra, and provides the mathematical deﬁnitions and concepts necessary for the later chapters. Chapter 2 presents the geometric applications of quaternions, and provides the ideas necessary to understand how quaternions can be used to represent both three- and four-dimensional values and geometric operations applied to them. Chapter 3 gives a detailed and comprehensive account of quaternion Fourier transforms, including their deﬁnitions, operator formulas and how they may be computed. Chapter 4 shows how quaternion Fourier transforms can be applied in signal and image processing.

2 The “adjoint” terminology was taken from Zhang’s 1997 paper [ZHA 97].

1 Quaternion Algebra

This chapter introduces the quaternion algebra H and presents some properties that will be useful in later chapters.

1.1. Deﬁnitions Quaternions are one of the four existing normed division algebras over the real numbers. Classically denoted by H in honor of Sir W.R. Hamilton who discovered them in 1843, they form a non-commutative algebra. A quaternion q ∈ H is a fourdimensional (4D) hypercomplex number and has a Cartesian form given by: q = a + ib + jc + kd

[1.1]

where a, b, c, d ∈ R are called its components. The three imaginary units i, j, k are square roots of −1 and are related through the famous1 relations: i2 = j 2 = k2 = ijk = −1 ij = −ji = k ki = −ik = j

[1.2]

jk = −kj = i A quaternion q ∈ H can be decomposed into a scalar part S(q) and a vector part V(q): q = S(q) + V(q)

[1.3]

1 These relations deﬁning the three imaginary units of an element of H were carved by Sir W. R. Hamilton on a stone of the Broome bridge in Dublin on 16 October 1843.

2

Quaternion Fourier Transforms for Signal and Image Processing

where S(q) = a and V(q) = q − S(q) = ib + jc + kd. Obviously, S(q) ∈ R and we will also refer to it as the real part of q, i.e. (q) = a. Now, q ∈ H will be called a pure quaternion if its real part is null, i.e. if S(q) = 0. The set of pure quaternions will be denoted as V(H), while the set of quaternions with a null vector part are trivially identiﬁed with elements of R, i.e. S(H) ≡ R. If q has a null vector part, V(q) = 0, then q is simply an element of R. To identify the different imaginary components of a quaternion q = a + ib + jc + kd, we will use the following notations: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩

i (q)

=b

j (q)

=c

k (q)

=d

[1.4]

so that a quaternion q ∈ H can be written as: q = (q) + i

i (q)

+j

j (q)

+k

k (q)

[1.5]

The Cartesian notation for a quaternion q ∈ H is, in fact, its expression in a speciﬁc 4D basis of the algebra H, namely in the basis {1, i, j, k}. Recall that, as an algebra, H possesses a vector space structure that allows the expression of any of its elements in terms of its components in a basis of H. The basis {1, i, j, k} is the most common and popular basis to express a quaternion. However, we may encounter some other bases later on in this book, leading to alternate notations for a quaternion q ∈ H. Before introducing these notations, we ﬁrst review some remarkable properties of quaternions. 1.2. Properties Here, we list some of the properties of quaternions that will be used throughout the book. From the algebra structure of H, the sum of two quaternions is trivial. Given two quaternions q and p, we have: q + p = S(q) + S(q) + V(p) + V(q)

[1.6]

Expressing the two quaternions in their Cartesian forms, q = a + ib + jc + kd and p = e + if + jg + kh, their sum is: q + p = (a + e) + i(b + f ) + j(c + g) + k(d + h)

[1.7]

Quaternion Algebra

3

and their product takes the form: qp

= (a + ib + jc + kd) (e + if + jg + kh) = ae − (bf + cg + dh) + a (if + jg + kh) + e (ib + jc + kd) + i (ch − dg) + j (df − bh) + k (bg − cf )

[1.8]

Using the scalar/vector notation, this product takes the following form: qp = S(q) S(p) − V(q) , V(p) + S(q) V(p) + S(p) V(q) + V(q) × V(p)[1.9] where ., . is the scalar product and × is the vector cross product. These are understood in the classical sense of the three-dimensional (3D) vector cross and inner products, which means that: V(q) , V(p) = bf + cg + dh

[1.10]

which is scalar valued, i.e. V(q) , V(p) ∈ R, and that: V(q) × V(p) = i(ch − dg) + j(df − bh) + k(bg − cf )

[1.11]

where the result is a pure quaternion, i.e. (V(q) × V(p)) ∈ V(H). A very noticeable property is that the product of two quaternions is not commutative so that in general: qp = pq

[1.12]

This can be inferred from the presence of the non-commutative cross product in [1.9]. Note, however, that the product of quaternions is associative so that for any three quaternions q, p, r ∈ H, the following is true: (qp) r = q (pr)

[1.13]

The norm of a quaternion q is deﬁned as: q = a2 + b2 + c2 + d2

[1.14]

A quaternion q ∈ H with q = 1 is said to be a unit quaternion. As previously mentioned, H is one of the four existing normed division algebras. As a result, given any two quaternions p, q ∈ H, then: qp = q

p

[1.15]

4

Quaternion Fourier Transforms for Signal and Image Processing

It can also be easily checked that qp = pq . A related quantity that will be used in the following is the modulus of a quaternion. It is deﬁned as the length of the quaternion in Euclidean 4D space. The modulus of q ∈ H is denoted by |q| and is expressed as: |q| =

a2 + b2 + c2 + d2 = q

1 2

[1.16]

Obviously, |q| ∈ R+ and |q| = 0 if and only if q = 0. Like the norm, the modulus of a product of two quaternions p and q has the following property: |pq| = |p| |q| = |qp|

[1.17]

Just as with the complex numbers, the conjugate of a quaternion q is obtained by negating its imaginary part. However, in H the imaginary part is 3D and consists of the entire vector part V(q). Denoted by q, the conjugate of q is thus deﬁned as: q = a − ib − jc − kd = S(q) − V(q) by:

[1.18]

It follows that the scalar and vector parts of any quaternion q ∈ H can be obtained 1 2

S(q) =

(q + q) and

V(q) =

1 2

(q − q)

[1.19]

Conjugation in H has the following property: q=q

[1.20]

In contrast to the complex case, conjugation is not an involution2 but an antiinvolution, such that for q, p ∈ H: qp = p q

[1.21]

that is, the order of the factors in a quaternion product is reversed by the conjugation operator. Note that the modulus (and also the norm) of a quaternion q can be expressed using the conjugate of q as: |q| =

qq

and

q = qq

2 An involution f : H → H is such that for q, p ∈ H: f (f (q)) = q f (p + q) = f (p) + f (q) f (pq) = f (p)f (q)

[1.22]

Quaternion Algebra

5

It is also possible to deﬁne involutions over H. Involutions are deﬁned with respect to a pure unit quaternion μ (μ ∈ V(H) and |μ| = 1). The most general case (together with many properties) is presented in [ELL 07d]. As special cases, one can choose μ as one of the standard basis elements of H, i.e. i, j or k. Involutions with respect to these three unit pure quaternions can be called canonical. Given a quaternion q, its three canonical involutions are: qi qj qk

= a + ib − jc − kd = a − ib + jc − kd = a − ib − jc + kd

= −iqi = −jqj = −kqk

[1.23]

Clearly, q and its three canonical involutions allow us to recover the four components a, b, c, d of q by linear combination. Now, the most general deﬁnition for involution is: q μ = −μqμ

[1.24]

where μ ∈ V(H) and |μ| = 1. Involutions in H possess many properties (see [ELL 07d] for details) among which, for q, p ∈ H and μ ∈ V(H) and μ2 = −1 (i, j and k are possible choices for μ): ⎧ ⎪ qp μ ⎪ ⎨ μ qμ ⎪ ⎪ ⎩q i j

= q μp μ [1.25]

=q =q

k

As H is a division algebra, any non-null quaternion possesses an inverse. The inverse of a given quaternion q ∈ H is given by: q −1 =

q q

[1.26]

where it can be easily checked that qq −1 = 1 because of [1.22]. Note that for a pure unit quaternion μ, ( μ = 1 and S(μ) = 0), the following holds: μ−1 = −μ. Now that we have introduced the inverse of a quaternion, we are ready to look at the ratio of two quaternions p and q. Ratios must be handled with care in H (indeed, in any non-commutative algebra), and it is preferable to avoid the p/q notation when possible, as it is ambiguous, since p/q can be interpreted as the product of p by q −1 ;

6

Quaternion Fourier Transforms for Signal and Image Processing

the notation p/q does not specify the order of the product so that it leaves the possibility for pq −1 and q −1 p. The ambiguity arises from the fact that in general: pq −1 = q −1 p

[1.27]

The above non-equality arises from the fact that: pq −1 =

pq q

while q −1 p =

qp q

[1.28]

Thus, it is important to consider ratios as products with the inverse and to take care of the order of the product. Now, the modulus of a ratio is an interesting quantity to consider, as it does not suffer from the order in the multiplication, due to the property of the modulus given in [1.17]. It thus follows that: q −1 p = pq −1 =

|p| |q|

[1.29]

It can be useful to write a given quaternion q ∈ H as a product of a scalar positive number (its modulus) and a unit quaternion. This can be done in the following way: q = |q| q˜ = |q|

a b c d +i +j +k |q| |q| |q| |q|

[1.30]

where we used the notation q˜ for the unit modulus version of q, i.e. q˜ = q/ |q| so that |˜ q | = 1. Note that the decomposition of q into the product of its modulus and its normalized version is unique. The normalized version of q, denoted by q˜, is also sometimes called a versor. Finally, it must be emphasized that if q is a pure quaternion, i.e. q ∈ V(H), then it is uniquely written as q = |q| μ where we have denoted q˜ = μ to highlight the fact that it is a pure unit quaternion. In [1.5], we introduced the Cartesian form of a quaternion q ∈ H, in which it is expressed using the sum of a real part (q), an i−imaginary part i (q), a j−imaginary part j (q) and a k−imaginary part k (q). This expression is a special case of the expansion of a quaternion over a 4D basis. The speciﬁc basis used in [1.5] is {1, i, j, k}. This is the classical basis used by most authors. Now, it is possible to use a different basis and it turns out that there is an inﬁnite amount of choices for a basis in H. Given two pure unit quaternions μ and η, i.e. μ2 = η 2 = −1, which are orthogonal to each other: μ⊥η, the set {1, μ, η, μη} is a basis of H. The inﬁnite number of possibilities arises from the inﬁnite number of possible choices for the

Quaternion Algebra

7

pair of orthogonal pure unit quaternions μ and η. Note that the orthogonality between pure unit quaternions should be understood as: S(μη) = μ, η = 0

[1.31]

which is the scalar part of the quaternion product of two pure quaternions (see the expression for the quaternion product in [1.9]). Also note that the orthogonality condition does not mean that the quaternion product between μ and η is equal to zero. On the contrary, as can be seen from [1.9], we have: μη = V(μη) = μ × η

[1.32]

which is indeed a pure unit quaternion orthogonal to both μ and η. This makes μη the third pure unit quaternion that, together with 1, μ and η, forms the quaternion basis. Quaternion bases play an important role in the quaternion notation as well as in the computation of quaternion Fourier transforms (QFT) with arbitrary axis, as will be presented in Chapter 3. The list of quaternion properties presented in this section is not intended to be exhaustive. More properties can be found in various references [CON 03, KAN 89, WAR 97, KUI 02, HAN 06]. 1.3. Exponential and logarithm of a quaternion Among the functions with quaternion-valued arguments that will be considered in the sequel, two will be of major importance: the exponential and logarithm functions. The former will be central to the study and use of QFTs in subsequent chapters. 1.3.1. Exponential of a pure quaternion The exponential function exp : V(H) → H can be deﬁned through its power series expansion; thus, for a given (non-null) pure quaternion ξ ∈ V(H) written as ξ = |ξ| ξ with ξ a pure unit quaternion3 (i.e. ξ ∈ V(H) and ξ 2 = −1) and |ξ| ∈ R+ , we can write: eξ =

+∞

+∞

n

ξn |ξ| ξ n = n! n! n=0 n=0

[1.33]

3 Just like any quaternion, a pure quaternion ξ ∈ V(H) can be uniquely written as ξ = |ξ| ξ. ξ is unit quaternion (in this case, pure) and |ξ| is its modulus.

8

Quaternion Fourier Transforms for Signal and Image Processing

where we make use of the fact that |ξ| and ξ commute as any quaternion commutes with a real number. Now, just like the familiar complex imaginary number I does, a pure unit quaternion behaves, when it is raised to an integer power of n, as: ξn =

(−1)k (−1)k ξ

if if

n = 2k n = 2k + 1

[1.34]

With this in mind, the exponential eξ simply is: eξ =

+∞ n=0

(−1)k

2k

+∞

2k+1

|ξ| |ξ| +ξ (−1)k (2k)! (2k + 1)! n=0

[1.35]

In this equation, the right-hand side terms are the series expansions of the classical cos and sin functions. We have established that for a pure quaternion ξ = |ξ| ξ, the exponential of ξ is: e|ξ|ξ = cos |ξ| + ξ sin |ξ|

[1.36]

This shows that the exponential of a pure quaternion can be expressed easily in terms of cosine and sine functions just as in the complex case. The difference lies in the axis ξ which is a unit pure quaternion, while the argument is the modulus of ξ. Thus, the exponential of a pure quaternion is a full quaternion, with real/scalar part cos |ξ| and vector part ξ sin |ξ|. It must also be noticed that the exponential of a pure quaternion is always of unit modulus: eξ = 1,

∀ξ ∈ V(H)

[1.37]

which is easily veriﬁed using [1.36] as well as property [1.34] with n = 2. Thus, we conclude that the exponential of a pure quaternion is a full unit quaternion. Now, the reciprocal property is interesting: any full unit quaternion can be expressed as the exponential of a pure quaternion4. This is a remarkable property that will be tackled in section 1.3.3 for the expression of the logarithm of a quaternion, and especially for the polar form and Euler formula over H (sections 1.4.1 and 1.4.1.1). 4 A full unit quaternion q is uniquely expressed as the exponential of a pure quaternion, i.e. q = exp(ξ) with ξ = |ξ| ξ, iff |ξ| ∈ [0, 2π[.

Quaternion Algebra

9

Another important property, which differs drastically from the complex exponential case, is that the product of two exponentials of pure quaternions5 is not an exponential with argument equal to the sum of the arguments of the original exponentials. This means that for α, β ∈ R+ and any two distinct pure unit quaternions μ, ν, we have: eνα eμβ = eνα+μβ

[1.38]

One should note that in the special case where μ = ν, the equality is fulﬁlled, meaning that for α, β ∈ R+ and μ a pure unit quaternion, one has exp (μα) exp (μβ) = exp (μ (α + β)). In fact, this is a well-known property of exponential functions over non-commutative algebras (for example, the exponentials of matrices) and a consequence of the Baker–Campbell–Hausdorff formula (see, for example, [GIL 08] for illustrations of this formula). The exponential function of a pure quaternion will be of use when considering polar representations of quaternions as well as when deﬁning QFTs. 1.3.2. Exponential of a full quaternion We have already introduced the exponential of a pure quaternion in section 1.3.1 for use later in the deﬁnition of Fourier transforms. Here, we consider the exponential of full quaternions, i.e. the function exp : H → H. Just as in section 1.3.1, the exponential function is directly deﬁned through its series expansion that is absolutely convergent. Thus, for a quaternion q ∈ H, its exponential is given by: eq =

+∞

qn n! n=0

[1.39]

Now, recalling that any quaternion can be written as q = S(q) + V(q), it follows directly that: eq = eS(q) eV(q)

[1.40]

leading to the special case of the exponential of a pure quaternion if S(q) = 0. Now, it is also possible to expand the eV(q) part into cos and sin, as V(q) ∈ V(H), leading to the following expression for the exponential of q: eq = eS(q) cos |V(q)| + V(q) sin |V(q)| 5 This generalizes to the exponential of full quaternions.

[1.41]

10

Quaternion Fourier Transforms for Signal and Image Processing

with V(q) = |V(q)| V(q), where V(q) is the normalized vector part of q, i.e. the versor associated with V(q). Note that due to the fact that eS(q) is real-valued, it commutes with the exponential of the vector part V(q). Finally, note that given a quaternion q ∈ H and a complex number z ∈ C, their exponentials eq and ez are isomorphic, provided that |q| = |z| and S(q) / |V(q)| = (z)/ (z). 1.3.3. Logarithm of a quaternion The logarithm of a quaternion q is simply deﬁned as the inverse of the exponential function, so that for q, p ∈ H, we have: p = ln q

[1.42]

which means that ep = q. However, it is possible to obtain an expression for the logarithm of q in terms of its elements. First, we recall that q ∈ H can be expressed as q = |q| q˜, where q˜ is the normalized version of q. Now, it follows directly that the logarithm of q is: ln q = ln |q| + ln q˜

[1.43]

Now, the second term on the right-hand side of the equation can be found to have an interesting literal expression. There is no doubt that ln |q| ∈ R; now, we are going to show that ln q is a pure quaternion. First, remember from section 1.3.1 that, as q is a full unit quaternion, it can be expressed as the exponential of a pure (non-unit) quaternion ξ : q˜ = eξ = e|ξ|ξ = eφq μq

[1.44]

where we used the following notation: φq = |ξ| μq = ξ

∈ R+ ∈ V(H)

[1.45]

This notation was chosen as φq and μq will be identiﬁed later in section 1.4. Now, with this notation, and by direct substitution of [1.44] into [1.43], we have: ln q = ln |q| + μq φq

[1.46]

where we use the fact that μq φq = φq μq as φq is real-valued. This expression of the logarithm of a quaternion can be seen as a generalization of the well-known logarithm

Quaternion Algebra

11

of a complex number. In particular, the fact that we ﬁxed the range of values to be taken by the argument φq between 0 and 2π ensures that the logarithm function with a quaternion argument is not multivalued. 1.4. Representations Here we describe several existing representations of quaternions that will be used throughout the book. 1.4.1. Polar forms 1.4.1.1. Euler formula In addition to the Cartesian form of q given in [1.1], there are several other representations for quaternions that have been introduced since their discovery. One of the most important notations, which was introduced by Hamilton himself and is called the polar form of a quaternion q, is the quaternion equivalent of what Richard Feynman called the “jewel” of complex numbers, namely the Euler formula. It encapsulates the link between geometry and algebra. For a quaternion q ∈ H, it reads as: q = |q| eμq φq = |q| cos φq + μq sin φq

[1.47]

where |q| is the modulus, μq is called the axis and φq is the phase/argument of q. The elements of the polar form are given in terms of the Cartesian components as: ⎧ |q| ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨μ

q

⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩φq

=

√

a2 + b2 + c2 + d2

bi + cj + dk =√ b2 + c2 + d2 √ b2 + c2 + d2 = arctan a

[1.48]

where μq ∈ V(H) and μq = 1. From this notation, it follows that unit quaternions have a modulus equal to 1, and pure quaternions (i.e. with a = 0) have a phase equal to π/2. The polar form allows us to link geometrical concepts such as rotations in 3D and 4D to the algebra of quaternions. This will be elaborated upon in Chapter 2. It is interesting to look at the special case of quaternions with unit modulus. In such a case, the Euler form simply reads as: eμφ = cos φ + μ sin φ

[1.49]

12

Quaternion Fourier Transforms for Signal and Image Processing

still with μ ∈ V(H) and |μ| =1. Note that cos φ ∈ R and μ sin φ is the vector part of this quaternion. As a result, taking its conjugate consists of simply negating the vector part, giving the following identity: cos φ − μ sin φ = e−μφ

[1.50]

where one can see that conjugating a unit quaternion consists of either reversing its axis or negating its phase. Note that this could have been guessed from the fact that for unit quaternions q ∈ H, i.e. |q| = 1, then q = q −1 . Now, due to expressions [1.49] and [1.50], we can express the sin and cos functions in the following way: cos φ sin φ

= 12 eμφ + e−μφ = − 12 μ eμφ − e−μφ

[1.51]

which will prove to be useful later in the study of QFTs in Chapter 3. 1.4.1.2. The Euler angle parameterization polar form Another polar form that was used in [BÜL 01] has a direct connection to Euler angles [ALT 86]. As there are 12 possible conventions when using Euler angles, there are equivalently 12 possible polar forms. Here, we illustrate the polar form with the XY Z convention as used in [BÜL 01]. With this convention, a quaternion q ∈ H can be expressed as: q = |q| eηi eκj eφk

[1.52]

The three angles of q, i.e. η ∈ [0, 2π), κ ∈ [0, π), φ ∈ [0, 2π), can be identiﬁed as three phases, which are related to Euler angles when |q| = 1. 1.4.1.3. The Cayley–Dickson form It can also be useful to consider a quaternion as a pair of complex numbers with speciﬁc multiplication rules. This is the idea behind the Cayley–Dickson (CD) form of a quaternion q, which can be obtained using the doubling procedure [KAN 89]. Thus, a quaternion q ∈ H has a CD form that reads: q = z1 + z2 j

[1.53]

where z1 = a + ib ∈ Ci and z2 = c + id ∈ Ci . Obviously, there is an inﬁnite number of CD forms, depending on the unit pure quaternion that is used to split the quaternion into two different planes. This will be detailed in section 1.4.1.4. All the quaternion properties could be rephrased by this notation. However, we will not do so as it is not of use in the following. As an example, we give the expression of the quaternion conjugation in the CD notation, which is: q = z1 − z2 j where denotes the classical complex conjugation.

[1.54]

Quaternion Algebra

13

1.4.1.4. The ortho-split or symplectic form A more general representation exists, in the spirit of the CD form, that allows the interpretation of a quaternion in terms of two complex numbers found in two non-intersecting two-dimensional (2D) planes in 4D space6. Consider a basis in H: {1, μ, ν, μν}, where μ2 = ν 2 = −1 and S(μν) = 0, i.e. μ ⊥ ν. μ and ν are pure unit quaternions (square roots of −1). As a result, μν ⊥ μ and μν ⊥ ν. Then, any quaternion q ∈ H can be written as: q = (a + b μ) + (c + d μ)ν

[1.55]

where the two above-mentioned planes are spanned by {1, μ} and {ν, μν}. This representation of a quaternion is sometimes referred to as its ortho-split representation or symplectic form7 as in [ELL 07c]. It is also known as the decomposition of q into its simplex part and perplex part. Using the quantities deﬁned in [1.55] the simplex part of q is (a + b μ) and the perplex part of q is (c + d μ). Later in this book, we will sometimes make use of the notation: q = qs + qp ν

[1.56]

with notations qs ∈ Cμ for the simplex part of q ∈ H and qp ∈ Cμ for its perplex part, both understood with respect to the axis ν satisfying μ⊥ν. Using this symplectic notation, a quaternion q can be written as: q = q+ + q−

[1.57]

where: q+ q−

= =

1 2 1 2

(q + μ q ν) (q − μ q ν)

[1.58]

which can be seen as a generalization of the expressions in [1.19] using a quaternion q and its conjugate. We also mention here the swap-rule for the symplectic notation, which will be of use in the even-odd split study of the QFT in section 3.1.3. If, in a symplectic 6 The two planes are called non-intersecting even though it is not strictly true as they intersect at the origin. However, as this is the only point in 4D space where the planes intersect, we will keep the terminology “non-intersecting”. 7 Note that the terminology symplectic used here is different from the classically used term “symplectic” used in differential geometry, topology or group theory where it designates a special non-degenerate 2-form.

14

Quaternion Fourier Transforms for Signal and Image Processing

decomposition, one needs to have the imaginary axis ν to be placed on the left, then one can simply write q as: q = (a + b μ) + (c − d μ)ν = qs + νqp

[1.59]

where the complex conjugation is used on qp despite its quaternionic nature as it is isomorphic to a complex number, i.e. qp ∈ Cμ . Note ﬁnally that choosing a basis in H to express a quaternion q induces a choice of symplectic representation. The symplectic form is thus a generalization of the CD form presented in section 1.4.1.3. 1.4.1.5. The polar Cayley–Dickson form Recently, a new representation was introduced for elements of H, called the polar Cayley–Dickson representation [SAN 10]. Given a quaternion q ∈ H, its polar CD form is: q = ρq eΦq j

[1.60]

with ρq ∈ C the complex modulus of q and Φq ∈ C its complex phase. This representation of a quaternion is unique. Its construction is given in detail in [SAN 10]. Here, we present the main expressions of the polar CD form, without details of the sign ambiguity that is known to exist during the construction of the polar CD form of a quaternion from its Cartesian form. For a complete discussion of this sign issue, one should refer to [SAN 10]. Now, given a quaternion q ∈ H, with CD form q = z1 + z2 j, its complex modulus and phase ρq and Φq are obtained by: ⎧ ⎪ ⎪ ⎨ρq

= z1

⎪ ⎪ ⎩ Φq

= − ln

|q| |z1 | z1 q j z1

[1.61]

Note that in the expression for Φq , the argument of the logarithm is a product between a complex number and a quaternion, meaning that the order is important and that q should be multiplied from the left by z1 . Note also that the product z1 q is rather 2 special in that if we replace q by its CD form, we get z1 q = |z1 | + z1 z2 j, which is a degenerate quaternion having a real part, an j part and an k part, but its i part is null. As a result, taking the logarithm of such a degenerate quaternion leads to a quaternion with only a j and a k part. The negation and right multiplication by j ﬁnally lead to the fact that Φq is a complex number from Ci . It is also interesting to note that any ortho-split form as introduced in section 1.4.1.3 will lead to another possible polar CD form based on the axis used for the split. An illustration of the usefulness of the polar CD form is made in section 4.3.4.

Quaternion Algebra

15

1.4.2. The Cj -pair notation In the style of the ortho-split notation or CD form, it is possible to think of a quaternion as an ordered pair of complex numbers and to manipulate them as such. Here, we mention this notation, called the Cj -pair associated with a quaternion, as it highlights a very special product, namely the bicomplex product of two quaternions. This product will be of use when studying convolution for the one-dimensional (1D) QFT in Chapter 4. Consider the quaternion q = a+ib+jc+kd ∈ H and let us write it as q = q1 +iq2 . In this way, q can be considered as the pair q = (q1 , q2 ) with q1 = a + jc and q2 = b + jd being elements of Cj . Obviously, it is a special case of an ortho-split decomposition, with axis i, where the quaternion q is simply understood as a pair of complex numbers. Now, with this notation, we could rewrite anything related to quaternion algebra. For example, the conjugate of q is q = (q1 , −q2 ). More interestingly, if we consider the two quaternions q = (q1 , q2 ) and p = (p1 , p2 ), then their (quaternion) product is: qp = (q1 p1 − q2 p2 , q2 p1 + q1 p2 )

[1.62]

and the inverse of q is simply: q −1 =

as:

q1

2

|q|

,

−q2 2

|q|

[1.63]

Now, the involutions of q with respect to the canonical axis of H can be expressed ⎧ i ⎪ ⎨q qj ⎪ ⎩ k q

= (q1 , −q2 ) = (q1 , q2 ) = (q1 , −q2 )

[1.64]

Indeed, all the quaternion operations can be handled with this notation, but we will make special use in the following of one particular operation that comes in a very useful form when using the complex pair notation: the bicomplex product denoted by . Given two quaternions in their complex pair forms q = (q1 , q2 ), p = (p1 , p2 ) ∈ H, their bicomplex product is: q

p = (q1 p1 − q2 p2 , q2 p1 + q1 p2 )

[1.65]

16

Quaternion Fourier Transforms for Signal and Image Processing

where we can see the difference with the standard quaternion product in [1.62]. A very remarkable property of the bicomplex product of two quaternions is that it is commutative. This can be easily seen from the expression of p q: p

q = (p1 q1 − p2 q2 , p2 q1 + p1 q2 )

[1.66]

and by remembering that p1,2 and q1,2 are complex numbers and thus pi qj = qj pi for any pair i, j = 1, 2. This leads to the following interesting property: q

p=p

q

∀p, q ∈ H

[1.67]

Note that the commutativity of this product is associated with the Cj -pairing notation. If we use a different way to divide a quaternion into two complex numbers (a Cμ pairing with μ a square root of −1), then we need to carefully look at how to deﬁne an associated commutative product. Rather than providing a lengthy study of the commutative product, we give here some properties of this special product over H, for any quaternion q = (q1 , q2 ) = (a + jc, b + jd), any Ci number z = (z)+i i (z) = ( (z), i (z)), any Cj number w = (w)+j j (w) = (w, 0) and Ck number s = (s)+k k (s) = ( (s), j k (s)). Then the following properties hold: ⎧ q ⎪ ⎪ ⎪ ⎪ ⎪ q ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ q ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ q ⎪ ⎪ ⎪ ⎨q ⎪ q ⎪ ⎪ ⎪ ⎪ ⎪ i ⎪ ⎪ ⎪ ⎪ ⎪ k ⎪ ⎪ ⎪ ⎪ i ⎪ ⎪ ⎪ ⎪ ⎪i ⎪ ⎪ ⎪ ⎪ ⎪ j ⎪ ⎪ ⎩ k

w z s i j k j j k i j k

= (q1 , q2 ) (w, 0) = (q1 w , q2 w) = (q1 , q2 ) ( (z), i (z)) = (q1 (z) − q2 i (z) , q2 (z) + q1 i (z)) = (q1 , q2 ) ( (s), j k (s)) = (q1 (s) − jq2 k (s) , q2 (s) + jq1 k (s)) = (q1 , q2 ) (0, 1) = (−q2 , q1 ) = (q1 , q2 ) (j, 0) = (q1 j, q2 j) = (q1 , q2 ) (0, j) = (−q2 j, q1 j) =j i=k = j k = −i = k i = −j = −1 = −1 =1

[1.68]

where all the given results can be checked by direct calculation using the bicomplex product expression given in [1.65].

Quaternion Algebra

17

As previously mentioned, the bicomplex product will not be encountered many times in this book, but it will be of use in section 4.3 when studying the convolution theorem for the QFT of complex-valued signals. The interested reader is referred to [PRI 91, ROC 04] for more materials on bicomplex numbers. 1.4.3. R and C matrix representations In addition to scalar representations, there exist two matrix representations over R and C that are isomorphic to H. Given a quaternion q ∈ H, its real matrix representation, denoted as MR (q) ∈ R4×4 , is given by8: ⎡

b a d −c

⎤ c d −d c ⎥ ⎥ a −b ⎦ b a

1 0 0 0

⎡ ⎤ 0 0 0 0 ⎢ −1 1 0 0 ⎥ ⎥ + b⎢ ⎣ 0 0 1 0 ⎦ 0 0 1 0

1 0 0 0

0 0 0 1

⎤ 0 0 ⎥ ⎥ −1 ⎦ 0

0 ⎢ 0 +c⎢ ⎣ −1 0

⎡ ⎤ 0 1 0 0 ⎢ 0 0 0 1 ⎥ ⎥ + d⎢ ⎣ 0 0 0 0 ⎦ −1 0 0 −1

0 0 1 0

0 −1 0 0

⎤ 1 0 ⎥ ⎥ 0 ⎦ 0

a ⎢ −b ⎢ MR (q) = ⎣ −c −d ⎡ ⎢ a⎢ ⎣

=

⎡

[1.69]

and for any quaternion q ∈ H, its complex matrix representation, denoted as MC (q) ∈ C2×2 , is of the form: MC (q) = =a

z1 −z2 z2 z1 1 0

0 1

+b

i 0 0 −i

+c

0 −1 1 0

+d

0 −i i 0

[1.70]

Comparing the real and complex matrix expressions for q with the Cartesian form, one gets a direct identiﬁcation of the imaginary units over H, i.e. i, j and k with their real 4 × 4 and complex 2 × 2 matrix forms. As a real vector space, quaternions are spanned by these four real or complex matrices. For both of these matrix representations, the quaternion product is identiﬁed with the matrix product over R4×4 and C2×2 . From a practical point of view, it is thus possible to manipulate 8 The transpose of the matrix representation given here is also valid.

18

Quaternion Fourier Transforms for Signal and Image Processing

quaternions through their matrix representations. However, as noted in section I.6.2, this is not optimal from a computational point of view as it requires much more computation and storage than necessary. Matrix representations of quaternions have been found to be of interest in the study of matrices with quaternion entries (for example [ZHA 97]), and in many more application ﬁelds. 1.5. Powers of a quaternion In many situations, we need to multiply quaternions together. Due to the fact that H is an algebra, the product of two quaternions is still a quaternion. Now, if we multiply a quaternion q ∈ H, say n times, by itself, it is interesting to look at how the components of the nth power of q can be expressed using the components of q. First, note that a quaternion commutes with itself. This may seem trivial, but it is interesting to note that the right and left multiplications are equivalent in this special case. Now, assuming that q has the polar form given in [1.47], it follows directly that: n

n

q n = |q| enμq φq = |q|

cos(nφq ) + μq sin(nφq )

[1.71]

This is simply a generalization of the de Moivre formula for quaternions [KAN 89]. The demonstration can be carried out easily by using the series expansion of the exponential and the property of the nth power of the pure unit imaginary μq (see [1.34]). Apart from its simplicity, this formula has useful consequences in the use of quaternions to represent rotations. Anticipating results from Chapter 2, multiplying9 n times by the same quaternion q will be shown to be equivalent to performing a rotation around the axis μq by an angle nφq /2. 1.6. Subﬁelds As previously mentioned, quaternions form a 4D algebra over R. A quaternion q with no vector part, V(q) = 0, is thus a real number, i.e. q ∈ R if V(q) = 0. Similarly, any quaternion with two of its three imaginary components equal to zero is homomorphic to a complex number. This means that for q ∈ H of the form: q = a + wμ

[1.72]

with μ = i, j, k, q is equivalent to a complex number. This is to say that numbers taking values in R ⊕ μR, with μ2 = −1, form an algebra homomorphic to complex numbers. q is said to be an element of Cμ . Obviously, we can see that such a construction is the restriction of H to a 2D plane spanned by {1, μ}, just as the 9 The exact way to perform rotations using quaternion multiplication is presented in detail in sections 2.1.2 and 2.1.6.

Quaternion Algebra

19

complex plane is classically known to be spanned by {1, I}. Now, the above-mentioned μ was required to be one of the three “classical” imaginary units from H. In fact, it may be any root of −1 in H so that the subﬁeld (equivalent to the complex ﬁeld) is now a different 2D plane inside the 4D space of H. Obviously, there are inﬁnitely numerous such subspaces that can be deﬁned within H. Arbitrarily choosing one such subﬁeld yields a new basis for the expression of any quaternion q ∈ H through the CD form given in section 1.4.1.4. This shows also that H contains two copies of C as subﬁelds.

2 Geometric Applications

The multiplication rules for the quaternion operators i, j and k provide a double meaning to quaternion numbers: on the one hand, they are geometric forms (e.g.vectors in four-dimensional (4D) space) and, on the other hand, they are transformation operators (e.g. rotation operators in three-dimensional (3D) space). Their usefulness as a tool in solving real-world problems is related to the various thought processes, or mental models, used to describe how they work. The geometric forms of signal and image samples upon which the models are based are inherently linked into the quaternion Fourier transform of those samples. There are three models: 3D and 4D vectors, spherical geometry and 3D projective space.

2.1. Euclidean geometry (3D and 4D) The ability of quaternions to represent geometric transformations has been known right from their discovery and it has been studied and applied in many areas [ALT 86, CON 03, HAN 06, KUI 02, WAR 97, SHO 85]. It is standard knowledge that 3D geometry can be well-handled using quaternions. However, 4D geometry is also a major concept that is inherent to quaternions. A very good reference on 4D geometry using quaternions is Coxeter’s paper [COX 46]. Material in the following sections can be found in [COX 46] together with [ELL 07a]. We introduce here the major geometrical transformations that will be of use in the later chapters. When considering 3D and 4D transformations, it will be assumed that full quaternions (i.e. elements of H) represent points or vectors in 4D space and pure quaternions (i.e. elements of V(H)) represent points or vectors in 3D space.

22

Quaternion Fourier Transforms for Signal and Image Processing

2.1.1. 3D reﬂections Given a 3D vector p, represented as a pure quaternion, the reﬂection of p across the plane deﬁned by its normal vector η (i.e. η ∈ V(H) and |η| = 1 so that it is a pure unit quaternion) is given by: Pη [p] = ηpη

[2.1]

Reﬂections are the simplest isometric transformations and the basic bricks for building up more elaborate isometric transformations, as is well known (see for example [COX 74, chapter III]). Finally, it is interesting to note that reﬂections are their own inverse, meaning that if the same reﬂection is applied twice, one returns to the original vector: Pη [Pη [p]] = η (ηpη) η = η 2 p η 2 = p

[2.2]

which is true due to the fact that η is a pure unit quaternion, i.e. η 2 = −1. 2.1.2. 3D rotations Given a 3D vector p, represented as a pure quaternion, the rotations of p in 3D space are encoded using one unit quaternion q (i.e. q ∈ H and |q| = 1) and are given by: Rq [p] = qpq

[2.3]

where the operator Rq [.] denotes the rotation with respect to q. As q is a unit quaternion, it can be expressed in its polar (Euler) form as q = eμq φq , with μq its axis and φq its angle/phase, as explained in section 1.4.1.1. The rotation operator Rq [.] consists of a rotation around the axis μq , with an angle of 2φq (twice the angle of q). Note that, as q is a pure quaternion, then q = q −1 , allowing the alternate notation for the rotation operator: Rq [p] = qpq −1 . It is noticeable that the inverse rotation operator R−1 q [.] is given by: R−1 q [p] = qpq

[2.4]

which is equivalent to a direct rotation encoded with q = q −1 . This means that the axis is the same but that the angle is exactly the opposite of the one in the rotation Rq [.].

Geometric Applications

23

A very good property of quaternions with respect to rotations is that, when composing such transformations, the resulting transformation can be expressed in terms of the product of the quaternions that represent the individual rotations. Indeed, the composition of two rotations, ﬁrst Rq1 [.] and then Rq2 [.], applied to a vector p takes the form: Rq2 [Rq1 [p]] = q2 (q1 pq1 ) q2 = q2 q1 p(q2 q1 )

[2.5]

which is equivalent to a rotation operator Rq2 q1 [.] encoded using the quaternion q2 q1 . Note that q2 q1 is a unit full quaternion as it is the product of two unit full quaternions. Obviously, the composition of rotations can be handled easily due to the fact that the group of unit quaternions (isomorphic to SU (2)) is the double universal covering of the 3D space rotation group SO(3) [ALT 86]. Generalizing the idea of composition, we can see that composing an arbitrary number of rotations using n different quaternions q1 , q2 , . . . , qn is equivalent to performing one rotation with the operator Rqn qn−1 ...q2 q1 [.]. A special case of composition is applying the same rotation n times. Applying the operator Rq [.] n times will look like: n

Rq [Rq [ . . . Rq [Rq [ p]] . . . ]] = q n p (q)

[2.6]

n-times

As we would suspect, this is a rotation of p around the axis μq with an angle nφq /2. The fact that the axis is kept the same during the successive rotation operations and that only the angle is affected is a consequence of the de Moivre formula given in [1.71]. Finally, it must be noted that, as is well known in Euclidean geometry, rotations can be expressed using reﬂections. Using the previously introduced reﬂection operator, we can see that two successive reﬂections with different axis/plane represent a rotation in the following way: Pη2 Pη1 [p] = η 2 (η 1 pη 1 ) η 2 = η 2 η 1 p (η 2 η 1 ) = Rη2 η1 [p]

[2.7]

where the unit quaternion that encodes the rotation is q = η 2 η 1 . We made use of the property shared by pure quaternion: η = −η together with the anti-involution property of the quaternion conjugation. Note that as η 1 and η 2 are pure unit quaternions, then q = η 2 η 1 is a unit full quaternion, which is necessary to have a 3D rotation. It is also necessary to have η 1 = η 2 in order to avoid degenerate cases. Also note that if η 1 ⊥ η 2 , then their product is a pure quaternion, i.e. η 1 η 2 = η 3 such that {η 1 , η 2 , η 3 } is a triad. The rotation obtained by succesive reﬂections with

24

Quaternion Fourier Transforms for Signal and Image Processing

axis η 1 and η 2 is then a special case of the rotation operator obtained in equation [2.7], where the angle of rotation is equal to π. Using the rotation operator notation from equation [2.3], it means that q = η 3 and the angle of q is φq = π/2. This involves that cos φq = 0 and thus that q = η 3 is a pure unit quaternion. The rotation operator consists of a rotation around η 3 with an angle of π. 2.1.3. 3D shears As introduced in [ELL 07a], we present here the quaternion formalism for the shear operator. In particular, two types of shears are studied here: axial-shears and beam-shears. The simplest shears are the axial-shears. They are deﬁned for a 3D vector p (again considered as a pure quaternion) using two parameters, a scalar α (the shear factor) and a pure unit quaternion μ3 (the shear axis), and take the form: Sα,μ3 [p] = p 1 −

α α μ − μ2 pμ1 2 3 2

[2.8]

with μ1 ⊥μ2 ⊥μ3 and μ1 μ2 = μ3 , so that {1, μ1 , μ2 , μ3 } form a quaternion basis. This means that the three pure unit quaternions form a right-handed triad. The shear operator reduces to the identity operator when α = 0 and reversing α inverts the shear operation. Now, the next level for shears is the beam-shear, that was introduced in [CHE 00]. This shear operator is parameterized by two scalar factors and a pure unit quaternion, part of a triad. The beam-shear operator reads: Sα,β,μ3 [p] = p

1 α + μ3 2 2

−

α μ p μ1 + 2 2

1 β + μ2 2 2

p−

β μ p μ3 2 1

[2.9]

where again μ1 , μ2 , μ3 form a right-handed triad. 2.1.4. 3D dilations The quaternion formalism for dilations was ﬁrst introduced in [ELL 07a]. Two types of dilations are presented here: axial- and radial-dilations. Considering a 3D vector p, the axial-dilation operator, denoted by Da , takes the form: a Dμ,α [p] =

1+α 1−α p+ μpμ 2 2

[2.10]

Geometric Applications

25

where μ is the axis of dilation/compression. The factor α controls the behavior of the transformation, leading to a compression if |α| < 1, a dilation if |α| > 1 and identity operator for α = 1. Negative values of α result in a reversal of the vector p. For the radial dilation case, the operator is denoted by D r and takes the form: r Dμ,α [p] =

2+α α p+ μpμ 2 2

[2.11]

This operator expands space outward from an invariant line deﬁned by μ. Here, dilations occur for α > 0 and compression occurs for α < 0, while α = 0 deﬁnes the identity operator. 2.1.5. 4D reﬂections As mentioned previously, quaternions have the ability to encode 4D geometric transformations [COX 46]. Here, we introduce the 4D reﬂection operator for a 4D vector p, encoded as a full quaternion, i.e. p ∈ H. The reﬂection of p across the hyperplane deﬁned by its normal 4D vector q (a unit full quaternion, i.e. q ∈ H and |q| = 1) is given by: Pq [p] = −q p q

[2.12]

Just like in the 3D case in section 2.1.1, the 4D reﬂection operator is its own inverse since: Pq [Pq [p]] = q(q p q) q = q q p q q = p

[2.13]

where we used the fact that q is a unit quaternion, i.e. qq = qq = 1. 2.1.6. 4D rotations The last Euclidean transformation is the 4D rotation of a 4D vector p (encoded as a full quaternion, i.e. p ∈ H). The 4D rotation operator is encoded using a pair of full unit quaternions q and r, i.e. q, r ∈ H and |q| = |r| = 1, such that S(q) = S(r) [COX 46]. The rotation operator takes the form: Rq,r [p] = q p r

[2.14]

where the rotation Rq,r [.] stands for a rotation about a plane (to be identiﬁed later) with angle 2φ if S(q) = S(r) = cos φ. Note that the scalar parts of q and r are cosines due to the fact that |q| = |r| = 1.

26

Quaternion Fourier Transforms for Signal and Image Processing

Now, just as in the 3D case, a rotation can be expressed as a combination of two reﬂections. In the 4D case, it arises in the following way when considering the product of two reﬂections across planes with h1 and h2 , respectively, as normal vector (with the requirement that h1 = h2 ): Ph2 [Ph1 [p]] = h2 (h1 p h1 ) h2 = h2 h1 p h1 h2

[2.15]

As h1 and h2 are pure unit quaternions, their product is a full unit quaternion. The quaternions h2 h1 and h1 h2 are different (multiplication order); however, they have their scalar part in common, which is not affected by the order of quaternion multiplication, so that S h2 h1 = S h1 h2 . This, in addition to the unit norm of the two quaternion products, is the requirement for the operator to be a 4D rotation. In fact, this rotation operator performs a rotation of angle φ with cos φ = S h2 h1 = S h1 h2 about a plane deﬁned by S h2 h1 p = S h1 h2 p = 0. We can also identify the obtained quaternions with the elements of the 4D rotation operator in section 2.14 like: h2 h1 = q and h1 h2 = r. The interested readers can refer to [COX 46, COX 74] for more on 4D Euclidean geometry. 2.2. Spherical geometry Geometry on a sphere resembles geometry on a plane. Thus, if A and B are two points on a sphere, then in analogy to a segment on a plane one can consider AB = B − A to be the distance between the points A and B, or the arc of a great circle running from point A to point B. The idea of the difference of two points as an arc joining them is completely analogous to the difference between two points in a plane being the vector joining them. Hence, two arcs will be equal if and only if they lie on the same great circle, and their lengths and directions are equal. Two equal arcs can be superimposed by moving one along the great circle on which they are located. One difference between planar and spherical geometry is that the addition of arcs is in general non-commutative. The addition of arcs is commutative only if they lie on the same great circle, in which case the addition reduces to a simple arithmetic sum. Branets [BRA 74], following Hamilton, considered the unit quaternion λ that can be represented in the form of a quotient of two equal-length pure quaternions or vectors in 3D-space, a and b, i.e., λ = b a−1 = cosφ + μsinφ.

[2.16]

The reciprocal of a is deﬁned through [1.26]. It can be shown that the following conditions hold: 1) the angle φ from the vector a to b satisﬁes φ ≤ π;

Geometric Applications

27

2) the unit vector μ is normal to the plane containing a and b; 3) the vectors a, b and μ span all of R3 . In other words, if we associate the point A on a sphere with the vector a, and the point B with the vector b, then the angle of the arc AB is in fact φ given in [2.16]. This situation is depicted in Figure 2.1. This implies that the quaternion product of two normalized quaternions corresponds to the geometric sum of the arcs of a sphere of radius |a| = |b|. So we have at our disposal the equivalence between arcs and unit quaternions, denoted as AB ⇐⇒ λb,a , which allows the correspondence AB + BC = AC ⇐⇒ λc,b λb,a = λc,a where we have included a third point C with its vector c.

Figure 2.1. Spherical geometry as unit length quaternion

Conversely, any unit quaternion λ deﬁnes an arc of a great circle, whose plane is normal to μ and whose arc length is φ; the position of the arc on its great circle is arbitrary, i.e. the arc may be slid along the circle. Conjugating the quaternion λ reverses the direction of the arc along the same great circle. Negating the quaternion −λ corresponds to the arc from −a to b. The special cases of an arc of length zero corresponds to a quaternion λ = 1 + 0μ so that the plane of the arc is arbitrary. The case where the arc length is π corresponds to λ = −1 + 0μ, as again the direction of the arc is arbitrary when moving a point

28

Quaternion Fourier Transforms for Signal and Image Processing

to its diametrically opposite point across the sphere – any direction will do. These associations provide a mechanism for encoding spherical triangles. The geometric sum of the three arcs AB, BC and CA forming a closed spherical triangle, when mapped to their unit quaternions, is: λa,c λc,b λb,a = (a c−1 )(c b−1 )(b a−1 ) = a (c−1 c)(b−1 b) a−1 = 1 which is the arc of zero length. In general, unit quaternions which represent arcs that are cyclically arranged edges of a closed polygon have a quaternion product equal to 1. 2.3. Projective space (3D) Full quaternions, those with non-zero scalar part, can be used to represent points in homogeneous coordinates. The quaternion identity: q =s+v=s 1+

v s

[2.17]

where s = S(q) and v = V(q). In this form, q can be used to represent a point located at the end of the vector p = v/s from the origin with weight s. Thus, in weighted-point form: q = s[1 + p] Hence, the set of quaternions can be associated with the real projective space P3 ∼ = R . This interpretation was discussed by Joly [JOL 05] and MacFarlane [MAC 06] in 1905 and 1906, respectively. Under this interpretation, a unit-weight point at the origin is denoted as q = 1[1 + 0i + 0j + 0k], and a weightless point at inﬁnity in the (x, y, z) direction is denoted as q = [0 + xi + yj + zk], which can be seen by letting s → 0 in [2.17]. The resulting pure quaternion can also be viewed as a translation vector when added to a weighted-point. 4

∼ V(H) and the algebra of the set of (weighted) The algebra of the set of vectors V = points P ∼ = H, both described using quaternion algebra, are different in subtle ways. For example, the direct sum of two vectors, v1 and v2 , follows the parallelogram law, whereas the direct sum of two (weighted) points, p1 and p2 , follows the principle of center-of-mass. This can be seen in the following equation: q1 + q2 = (s1 + s2 )[1 + (s1 p1 + s2 p2 ) / (s1 + s2 )] weight

center-of-mass

Geometric Applications

29

Likewise, scaling a vector v ∈ V by α ∈ R changes its length but leaves its direction unchanged, whereas scaling a point p ∈ P leaves its position unchanged but changes its weight since αq = αs[1 + p]. The set of vectors is closed with respect to quaternion addition and subtraction; however, the set of weighted-points is not closed. For example, the difference between two equal-weight points, q1 and q2 , is the weightless point (i.e. the vector) that joins them as: q1 − q2 = s [1 + p1 ] − s [1 + p2 ] = 0 + (p1 − p2 ) Table 2.1 provides a summary of how various quaternion algebraic operations map elements between point and vector subsets. Mapping αV → V αP → P V ±V →V P ±V →P P +P →P P − P → V or P VV → V or P PP → V or P

Projective space interpretation Scales vector length, direction unchanged. Scales point weight, position unchanged. Parallelogram law. Translates point in direction of vector. Principle of center-of-mass. V iff weights are equal. V iff vectors are perpendicular. V iff p1 , p2 = 1.

Table 2.1. Summary of point (P ∼ = H) and vector (V ∼ = V(H)) mappings with α ∈ R

Quaternion conjugation, by deﬁnition, negates the vector part of a quaternion. ¯ = −v. In Consequently, conjugation and negation of a vector are the same, i.e. v contrast, the conjugate of a weighted-point results in the point’s position being reﬂected across the origin, whereas negation changes the sign of the point’s weight, but its position remains the same, i.e.: s [1 + p] = s [1 − p] = −s [1 + p2 ] . Not withstanding these differences, some algebraic combinations behave similarly for vectors and points. For example, given two vectors, v1 and v2 , their linearly weighted sum: vα = (1 − α) v1 + (α) v2 traces a line from v1 to v2 as α varies from 0 to 1, i.e. it draws a line diagonally across the parallelogram deﬁned by v1 and v2 . Likewise, given two points, p1 and p2 , their linearly weighted sum: qα = (1 − α) s1 [1 + p1 ] + (α) s2 [1 + p2 ] = s1 [1 + p1 ] + s2 [1 + p2 ]

30

Quaternion Fourier Transforms for Signal and Image Processing

also traces the same path across the diagonal of the parallelogram deﬁned by points p1 and p2 . However, for a constant rate of change in α the weighed points are not in general equally spaced, as it is in the vector case, due to the relative weight of the two points. The points would be equally spaced only if their weights are equal. See Figure 2.2 for an example comparison where p2 has a weight three times that of p1 .

Figure 2.2. Linear weighted sum of points (top) versus vectors (bottom)

Geometric Applications

31

Even though points and vectors are disjoint subsets, they share a common framework when viewed as a projective space. A vector represents a point at inﬁnity along its direction so that geometrically all the points at inﬁnity have been included into a single 3-space model. We may ask, “what advantage does this projective space model of separating points and vectors grant us?” We saw in section 2.1 the use of various linear quaternion functions that map vectors into vectors by way of rotations, reﬂections, shears and dilations. This set of linear operations can now be expanded into projections (including translations). To understand this, we must digress and discuss the most general form of linear quaternion equation. 2.3.1. Systems of linear quaternion functions Real linear functions take the monomial form: f (x) = mx, where m, x ∈ R. Linear combinations, i.e. direct sums or compositions, of such functions can always be reduced to this same form. In contrast, quaternion linear functions have the multinomial form: P

f (q) =

mp q np

[2.18]

p=1

where all factors are quaternion-valued, i.e., q, mp , np ∈ H. We would expect under functional compositions that the number of terms in the summation can become arbitrarily large due to quaternion multiplication being non-commutative. However, this general multinomial can always be reduced to its quaternary canonical form: f (q) = A q + B q i + C q j + D q k,

[2.19]

where A, B, C, D ∈ H. The linear sum of two such functions is given by the component-wise addition. Let f1 (q) = A1 q + B1 q i + C1 q j + D1 q k and f2 (q) = A2 q + B2 q i + C2 q j + D2 q k then: f1 (q) + f2 (q) = A3 q + B3 q i + C3 q j + D3 q k where A3 = A1 + A2 , B3 = B1 + B2 , C3 = C1 + C2 and D3 = D1 + D2 .

32

Quaternion Fourier Transforms for Signal and Image Processing

The composition of two linear functions, f2 (f1 (q)), is given by: f2 (f1 (q)) = A3 q + B3 qi + C3 qj + D3 qk where: A3

= A2 A1 − B2 B1 − C2 C1 − D2 D1 ,

B3

= A2 B1 + B2 A1 − C2 D1 + D2 C1 ,

C3

= A2 C1 + B2 D1 + C2 A1 − D2 B1 ,

D3

= A2 D1 − B2 C1 + C2 B1 + D2 A1 .

Careful examination of this composition rule reveals that it has the same structure as the standard quaternion multiplication. The composition is, of course, not commutative: f2 (f1 (q)) = f1 (f2 (q)). This is because geometrical operations do not commute (the composition of two rotations, for example, gives different results, in general, according to the order in which the rotations are carried out). The equivalence between invertible functions in quaternary canonic form and the general linear group of 4 × 4 invertible real matrices, GL4 (R), will be used later. Let p = p0 + p1 i + p2 j + p3 k and q = q0 + q1 i + q2 j + q3 k be two arbitrary quaternions. Using standard hypercomplex operator product rules, the components of the quaternion product pq = q0 + q1 i + q2 j + q3 k = q are: q0

= p0 q0 − p1 q1 − p2 q2 − p3 q3 ,

q1

= p1 q0 + p0 q1 − p3 q2 + p2 q3 ,

q2

= p2 q0 + p0 q2 + p3 q1 − p1 q3 ,

q3

= p3 q0 + p0 q3 − p2 q1 + p1 q2 .

By gathering terms into matrix-vector notation: ⎡ ⎤ ⎡ ⎤⎡ ⎤ q0 p0 −p1 −p2 −p3 q0 ⎢q1 ⎥ ⎢ p1 p0 −p3 p2 ⎥ ⎢ q1 ⎥ ⎢ ⎥=⎢ ⎥⎢ ⎥ ⎣q2 ⎦ ⎣ p2 p3 p0 −p1 ⎦ ⎣ q2 ⎦ q3 p3 −p2 p1 p0 q3

[2.20]

we can deﬁne an equivalent matrix-vector form for the quaternion product. Let [[p]] and [q] denote the matrix and vector equivalences, respectively. Since the quaternion product is bilinear, we can instead place the components of p into the vector and the components of q into the matrix, but in doing so the lower right 3 × 3 submatrix is

Geometric Applications

33

transposed from the one given above. Denote this transmuted form of the matrix as † [[q]] . Therefore, any quaternion product has two equivalent matrix-vector forms: †

pq ↔ [[p]] [q] = [[q]] [p] . Hence, any triple-product pqr may be arbitrarily reordered in matrix form as †

pqr ↔ [[p]] [[q]] [r] = [[p]] [[r]] [q] The linear quaternary equation q = Aq + Bqi + Cqj + Dqk can be written in matrix-vector form as †

†

[q ] = [[A]] + [[B]] [[i]] + [[C]] [[j]] + [[D]] [[k]]

†

[q] = [[Q]] [q]

[2.21]

†

Now, both a matrix [[p]] and a transmuted matrix [[q]] contain only four degrees of freedom: one for each component of the quaternion. Careful examination of the four terms in the above matrix equation reveals that these matrices are independent, giving the total of 16 degrees of freedom necessary for an arbitrary matrix. It is straightforward to start with an arbitrary matrix [[Q]] and solve for the elements of the four quaternions {A, B, C, D}. For details of this analysis, see [ELL 07b]. 2.3.2. Projective transformations Having established the connection between GL4 (R), the general linear group of 4 × 4 invertible real matrices, we can now discuss the various types of geometrical transformation possible using a linear quaternion function. We draw here on the classic text of Bruce Meserve [MES 83] which clearly sets out the possibilities (although not using a quaternion formalism as given here). The most general transformation is projective. The afﬁne transformations are a subset of the projective transformations with the property that straight lines are preserved (that is, remain straight) under the transformation. Reﬂection, scaling and shear are examples of afﬁne transformations. Similarity transformations are a subset of the afﬁne transformations with the additional property that angles are preserved (thus geometric ﬁgures are transformed into similar geometric ﬁgures). Uniform scaling (that is, scaling along all axes by the same scale factor) has this property. Euclidean transformations are a subset of the similarity transformations. Translation and rotation are examples. They preserve angles as well as straight lines. Two additional subsets of the afﬁne transformations “cut across” the other subsets. The isometries preserve the distance between points. The linear transformations are those that obey superposition. In geometric terms, it means we may decompose a vector or point into the sum of two or more vectors or points and obtain the same

34

Quaternion Fourier Transforms for Signal and Image Processing

results by transforming the decomposed parts separately as by transforming the undecomposed vectors or points. Now a very signiﬁcant point is that all these transformations are linear when expressed as operations on a projective space, even if non-linear in the 3D Euclidean space of the signal samples or image pixels. All of these linear transforms can be encoded using real 4 × 4 matrices as shown in Figure 2.3, which were shown to be equivalent to linear quaternion functions in section 2.3.1.

Figure 2.3. Projective maps using real 4 × 4 matrix

3 Quaternion Fourier Transforms

Quaternion Fourier transforms provide elegance and expressive power in the analysis of vectorvalued signals and images. There is, however, a cost – an overwhelming number of transform deﬁnitions. This chapter provides some of the possible quaternion Fourier transforms deﬁnitions and their properties which appear best suited to vector-valued signals.

The classical (complex) Fourier transform of a one-dimensional (1D) complexvalued function f (t) : R → C can be deﬁned as: F{f }(ω) = F (ω) =

∞ −∞

f (t) e−Iωt dt,

[3.1]

where I is the complex square root of −1. The inverse transform is deﬁned by conjugating the exponential kernel and adding a scale factor as in F −1 {F }(ω) = f (t) =

1 2π

∞ −∞

F (ω) e+Iωt dω.

[3.2]

Symbolically, we write the signal and its transform pair as f (t) F (ω). The extension of the complex Fourier transform to two dimensions is straightforward. Given a signal f (x, y): R2 → C, the two-dimensional (2D) complex Fourier transform is F{f }(ω, ν) = F (ω, ν) =

∞

∞

−∞

−∞

f (x, y)e−I(ωx+νy) dx dy

36

Quaternion Fourier Transforms for Signal and Image Processing

The index space of the transformed function is of course 2D, like the original function. As in the 1D case, the inverse transform is deﬁned by conjugating the exponential kernel and adding a scale factor as in F −1 {F }(x, y) = f (x, y) =

1 (2π)2

∞

∞

−∞

−∞

F (ω, ν)e+I(ωx+νy) dω dν.

A sufﬁcient (though not necessary) condition that a function f (t) has a complex Fourier transform is that the function should be absolutely integrable, i.e.,

R

|f (t)| dt < ∞

[3.3]

The set of such functions is denoted by L1 (R; C). To be able to recover the original function from the inverse Fourier transform, the function must be of bounded variation, which means that f (t) can be expressed as the difference of two bounded monotonically increasing functions. Alternately, if the function is square integrable, i.e., 2

R

|f (t)| dt < ∞,

[3.4]

then both the Fourier transform and its inverse exist. Square integrable functions are denoted by L2 (R; C). To generalize the Fourier transform to quaternion-valued signals, the approach taken in this work is to replace the single complex square root of −1 = I 2 with one of many possible quaternion square roots of −1 = μ2 . In the quaternion domain, μ is taken from the set of all pure unit length quaternions, namely μ ∈ {q ∈ V(H) : q = 1}. This provides an additional degree of design freedom, beyond the complex Fourier transform, when treating problems in vector signal processing and can be chosen to advantage for each application. The ability to choose the root of −1 is particularly useful when dealing with 2D signals where each integral may be associated with its own root, or not. For purposes of later comparison, Table 3.1 contains the basic properties of the complex Fourier transform used in image and signal processing. Table 3.1 contains two operator formulas, which are the building blocks of image and signal processing: convolution and correlation. The formal deﬁnition of a convolution of two functions is (f1 ∗ f2 )(t) =

+∞ −∞

f1 (τ )f2 (t − τ ) dτ.

Quaternion Fourier Transforms √1 2π

+∞ −∞

eIωt F (ω) dω = f (t)

Time-domain Time reversal Complex conjugation Reversed conjugation

f (t) f (−t) f (t) f (−t)

F (ω) =

√1 2π

+∞ −∞

F (ω) F (−ω) F (−ω) F (ω)

e−Iωt f (t) dt Frequency-domain Frequency reversal Reversed conjugation Complex conjugation

Purely real f (t) ∈ R Purely imaginary f (t) ∈ IR Even, symmetric f (t) = f (−t) Odd, anti-symmetric f (t) = −f (−t)

F (ω) = F (−ω) Even, symmetric F (ω) = −F (−ω) Odd, anti-symmetric F (ω) ∈ R Purely real F (ω) ∈ IR Purely imaginary

Time shift Modulation Time scale

F (ω) e−Iωτ F (ω − ν) 1 F(ω ) |α| α F (αω)

Linearity Multiplication Convolution Correlation Delta function Shifted delta

f (t − τ ) f (t) eIνt f (αt) 1 f ( αt ) |α| αf1 (t) + βf2 (t) √ f1 (t) f2 (t) 2π (f1 ∗ f2 )(t) (f1 f2 )(t) √ δ(t) √2π δ(t − τ ) 2π eIνt 1 dn f (t) dtn n

nth derivative

t f (t)

37

Modulation Frequency shift Frequency scale

αF1 (ω) + βF2 (ω) F √1 (ω) ∗ F2 (ω) √2π F1∗(ω) F2 (ω) 2π F1 (ω)F2 (ω)

Convolution Multiplication

1 e√−Iωτ √2π δ(ω − ν) 2π δ(ω)

Shifted delta Delta function

(Iω)n F (ω) dn I n dω n F (ω)

nth derivative

N OTE.– f (t) : R → C, α, β ∈ R. Table 3.1. 1D Complex Fourier transform pairs

Likewise, the cross-correlation of two functions is deﬁned as:

(f1 f2 )(t) =

+∞ −∞

f1 (τ )f2 (t + τ ) dτ.

The correlation is superﬁcially similar to convolution, but it has a different physical interpretation. The generalization from complex to quaternion algebra comes at a cost: neither the correlation nor the convolution are commutative, and which function is shifted (f1 or f2 ) changes the results. In the complex domain, these differences are immaterial, but with quaternions it means that the formulas given in Table 3.1 do not work: instead we must ﬁnd the correct quaternion generalizations, as discussed in section 3.1.5.

38

Quaternion Fourier Transforms for Signal and Image Processing

3.1. 1D quaternion Fourier transforms In the complex case, the sign of the exponents of [3.1] and [3.2] must differ. Which of the two you choose to be positive does not matter, so long as you keep to the rule. If someone else has used the opposite choice, the Fourier pair calculated from a given function will be the complex conjugate of that given by your choice. This simple rule does not extend to the quaternion case. Still further deﬁnitions for the quaternion Fourier transform (QFT) arise from the fact that quaternion multiplication is non-commutative; hence, the placement of the exponential kernel in the transform deﬁnition becomes meaningful. This situation gives rise to four unique 1D QFTs. 3.1.1. Deﬁnitions The 1D right-sided 1 QFTs of a quaternion-valued signal f (t) : R → H are: R F∓μ {f }(ω) = F R (ω) = κ−

∞ −∞

f (t) e∓μωt dt

[3.5]

where the −μ version is referred to as the forward transform and the +μ version is the reverse transform, which is distinct from their inverse transforms given, respectively, by −R F∓μ {F R }(t) = f (t) = κ+

∞ −∞

F (ω) e±μωt dω.

[3.6]

In both equations, μ ∈ H is a pure unit quaternion so that μ2 = −1. The scale −1 factors must satisfy: κ+ κ− = (2π) . If κ− = κ+ , then the transform is a unitary transform with the consequence that the transform preserves the energy of the original signal f (t) when it is mapped to the frequency domain2. Likewise, the left-sided QFTs of a quaternion-valued signal, f (t) : R → H are: L F∓μ {f }(ω) = F L (ω) = κ−

∞ −∞

e∓μωt f (t) dt,

[3.7]

1 In this book, we use the term right-sided to refer to the placement of the exponential relative to the signal (the exponential is on the right for a right-sided transform). This is consistent with usage in the QTFM library [SAN 13b], where the transform handedness is selected by a parameter with values “L” or “R”. Care is needed to check whether this convention is used when reading other quaternion literature, since the convention may differ. 2 The scale factor requirement also holds for complex Fourier transforms.

Quaternion Fourier Transforms

39

and their inverses are given by: −L F∓μ {F L }(ω) = f (t) = κ+

∞ −∞

e±μωt F L (ω) dω.

[3.8]

When the context makes it clear, or when either a left- or a right-sided transform may be used, the superscript L/R will be dropped. Not all quaternion-valued functions can be Fourier-transformed. They are transformable if they fulﬁll certain conditions, known as the Dirichlet conditions. We need not use ﬁrst principles to discover which functions are transformable. Instead, we can borrow from the classic complex case as follows. Split f (t) via symplectic decomposition (equation [1.55]) into simplex and perplex parts with respect to μ, the axis of the transform, as: f (t) = fs (t) + fp (t) μ2 , where μ2 ⊥ μ and fs , fp ∈ Cμ . The left-sided QFT of f (t) becomes: F L (ω) = κ−

R

e−μωt fs (t) dt + κ−

R

e−μωt fp (t) dtμ2 .

[3.9]

so that: F L (ω) = Fs (ω) + Fp (ω)μ2 ,

[3.10]

where: Fs (ω) = κ−

R

e−μωt fs (t) dt,

and

Fp (ω) = κ−

R

e−μωt fp (t) dt.

Likewise, the right-sided QFT of f (t) becomes: F R (ω) = κ− = κ−

R

R

fs (t)e−μωt dt + κ− fs (t)e−μωt dt + κ−

R

R

fp (t)μ2 e−μωt dt

[3.11]

fp (t)e+μωt dt μ2

so that: F R (ω) = Fs (ω) + Fp (−ω)μ2 .

[3.12]

40

Quaternion Fourier Transforms for Signal and Image Processing

Due to the equivalence between Cμ and the standard complex numbers, and the fact that the QFT can be decomposed into a sum of complex subﬁeld functions, the existence and invertibility of the QFT is inherited from the complex case. So, one sufﬁcient (though not necessary) condition is that if f (t) is absolutely integrable (i.e. f (t) ∈ L1 (R; H)) and if it is of bounded variation on every ﬁnite interval, then F L (ω) (and F R (ω)) exists and f (t) can be recovered from the inverse Fourier transform relationship at each point at which f (t) is continuous. Since Fs (ω), Fp (ω) ∈ Cμ , they can be computed using the complex Fourier transforms of fs (t) and fp (t). This is done by replacing μ by I in fs and fp , taking their complex Fourier transforms, and then back-substituting μ for I in the results. This is also the basis of a numerical implementation method using complex FFTs, which is discussed in section 3.3.1.3. 3.1.2. Basic transform pairs There are several basic relationships that are of use in manipulating Fourier pairs, summarized in Table 3.2. For the most part, their proofs are elementary. The art of determining the Fourier transform of a function is in using these relationships. Linearity: of the transform integrals implies αf1 (t) + βf2 (t)

αF1 (ω) + βF2 (ω),

where α, β ∈ Cμ and fn (t) operators.

[3.13]

Fn (ω). Hence, all four 1D QFT deﬁnitions are linear

Scaling: either time or frequency results in the following pairs: f (αt)

1 F |α|

ω α

and

F (αω)

1 f |α|

t α

,

[3.14]

where α ∈ R. Shifting: The Fourier pairs for shifts, in either time f (t−τ ) or frequency F (ω−ν), are unique for each Fourier deﬁnition, namely: f (t − τ )

R F−μ (ω) e−μωτ

and

R F−μ (ω − ν)

f (t) eμνt

[3.15]

f (t − τ )

L e−μωτ F−μ (ω)

and

L F−μ (ω − ν)

eμνt f (t)

[3.16]

or

Quaternion Fourier Transforms

41

depending on the chirality of the transform used. As with the complex Fourier transforms, shifts and modulation operations are inter-related between the time and frequency domains, but the handedness of the transform must be taken into account. √1 2π Right-sided: √12π

Left-sided:

+∞ −∞ +∞ −∞

eμωt F (ω) dω = f (t) F (ω)e

Time-domain Time reversal Quaternion conjugation Reversed conjugation Real Pure μ-imaginary Even, symmetric Odd, anti-symmetric Simplex Pure perplex

μωt

F (ω) =

dω = f (t)

F (ω) =

f (t) f (−t)

F (ω) F (−ω)

√1 2π √1 2π

+∞ −∞ +∞ −∞

Frequency-domain Frequency reversal

f (t)

See section 3.1.4 See section 3.1.4

f (t) ∈ R f (t) ∈ μR f (t) = f (−t) f (t) = −f (−t) f (t) ∈ Cμ2 f (t) ∈ Cμ2 μ

F (ω) = F (−ω) Even, symmetric F (ω) = −F (−ω) Odd, anti-symmetric F (ω) ∈ R Real F (ω) ∈ μR Pure μ-imaginary F (−ω) = F (ω) F (−ω) = −F (ω)

f (t − τ ) f (t − τ ) Modulation (left-sided transform) eμνt f (t) Modulation (right-sided transform) f (t) eμνt Time scale f (αt) 1 f ( αt ) |α|

e−μωτ F (ω) F (ω) e−μωτ F (ω − ν) F (ω − ν) 1 F(ω ) |α| α F (αω)

Linearity Convolution Correlation

αF1 (ω) + βF2 (ω)

αf1 (t) + βf2 (t) (f1 ∗ f2 )(t) (f1 f2 )(t) √ δ(t) √2π δ(t − τ ) 2π eμνt 1

nth derivative nth derivative (Left-sided transform) (Right-sided transform)

f (t)e−μωt dt

f (−t)

Time shift

Delta function Shifted delta

e−μωt f (t) dt

dn f (t) dtn dn f (t) n dt n

t f (t) tn f (t)

(Left-sided transform) (Right-sided transform) Frequency shift Frequency shift Frequency scale See section 3.1.5 See section 3.1.5

1 e√−μωτ √2π δ(ω − ν) 2π δ(ω) (μω)n F (ω) F (ω) (μω)n dn μn dω n F (ω) dn F (ω) μn dω n

Shifted delta Delta function (Left-sided transform) (Right-sided transform) nth derivative

N OTE.– f (t) : R → H, μ2 ⊥ μ, and α, β ∈ R. Table 3.2. 1D forward quaternion Fourier transform pairs

nth derivative

42

Quaternion Fourier Transforms for Signal and Image Processing

Derivatives: In a similar manner, derivatives, with respect to either time or frequency, have a chirality as: dn f (t) dtn

R F−μ (ω) (μω)n

and

dn R F (ω) μn dω n −μ

tn f (t)

[3.17]

dn f (t) dtn

L (μω)n F−μ (ω)

and μn

dn L F (ω) dω n −μ

tn f (t).

[3.18]

or

Impulse functions: The QFT of the Dirac delta function δ(t) mimics the complex transform exactly and is identical for all QFT forms: √ δ(t) 2π

1

and

1

√

2π δ(ω).

[3.19]

Likewise, the shifted delta function is common across all QFT forms as: √ δ(t − τ ) 2π

e−μωτ

and

eμνt

√

2π δ(ω − ν).

[3.20]

In the last equation above, if the reverse transforms are used, then the sign of μ is changed. 3.1.3. Decompositions In the time domain, complex-valued functions f (t) may happen to have one or more special symmetries. It might be purely real or purely imaginary, or it might be even, f (t) = f (−t), or odd, f (t) = −f (−t). These time-domain symmetries lead to symmetries in the frequency domain. Put another way, any function f (t) can always be decomposed into even-odd or real-imaginary or both parts. Decomposition of 1D QFTs is helpful in simplifying QFT pairs and studying the relationships among the various deﬁnitions. 3.1.3.1. Even-odd split In general, the function f (t) is quaternion-valued; denoting by fs (t) and fp (t) its simplex and perplex parts with respect to the transform axis μ, we have3: f (t) = fs (t) + fp (t)μ, where fs (t), fp (t) ∈ Cμ2 , μ2 ⊥ μ,

[3.21]

3 It should be noted that μ2 ⊥ μ allows for design freedom. If μ2 is ﬁxed, then the transform axis μ can be steered as needed in the plane perpendicular to μ2 . Conversely, if μ is ﬁxed, then μ2 can be steered in the plane perpendicular to μ.

Quaternion Fourier Transforms

43

and its (right-handed) QFT is F (ω) = Fs (ω) + Fp (ω)μ, where Fs (ω), Fp (ω) ∈ Cμ2

[3.22]

where S(ω) =

κ−

Fp (ω) = −κ−

R

R

[fs (t) cos ωt + fp (t) sin ωt] dt,

[3.23]

[fs (t) sin ωt − fp (t) cos ωt] dt.

[3.24]

Likewise, using the QFT inversion formula; with F (ω) = Fs (ω) + Fp (ω)μ, we obtain: fs (t) = κ+ fp (t) = κ+

R

R

[Fs (ω) cos ωt − Fp (ω) sin ωt] dω,

[3.25]

[Fs (ω) sin ωt + Fp (ω) cos ωt] dω

[3.26]

Now under the special condition that f (t) is strictly μ2 -complex, i.e. fp (t) = 0, then the simplex and perplex parts of F (ω) are given by: Fs (ω) = κ−

R

fs (t) cos ωt dt and Fp (ω) = −κ−

R

fs (t) sin ωt dt

[3.27]

From which we conclude that Fs (ω) is even and Fp (ω) is odd, i.e. Fs (−ω) = Fs (ω) and Fp (−ω) = −Fp (ω).

[3.28]

Therefore: R R (ω), if f (t) ∈ C , μ ⊥ μ. F−μ (−ω) = F−μ μ2 2

[3.29]

R R (ω), then f (t) = 0, i.e. f (t) ∈ C . Conversely, if F−μ (−ω) = F−μ p μ2

Starting instead with f (t) being entirely perplex, i.e. fs (t) = 0, so that f (t) = fp (t)μ, then the simplex and perplex parts of F (ω) are given by Fs (ω) = κ−

R

fp (t) sin ωt dt and Fp (ω) = κ−

R

fp (t) cos ωt dt

[3.30]

44

Quaternion Fourier Transforms for Signal and Image Processing

Therefore, Fs (ω) is odd and Fp (ω) is even, i.e. Fs (−ω) = −Fs (ω) and Fp (−ω) = Fp (ω).

[3.31]

Therefore: R R (ω), if f (t) ∈ C , μ ⊥ μ. F−μ (−ω) = −F−μ μ2 2

[3.32]

R R Conversely, if F−μ (−ω) = −F−μ (ω), then fs (t) = 0, i.e. f (t) ∈ Cμ2 .

Imposing further restrictions on μ2 -complex f (t), if we assume f (t) is an even function of time, i.e. f (−t) = f (t), then Fp (ω) = 0. Conversely, if the QFT is a μ2 complex simplex function, i.e. Fp (ω) = 0, then f (−t) = f (t), so that f (t) is even. However, if we assume instead f (t) is μ2 -complex and an odd function of time, i.e. f (−t) = −f (t), then Fs (ω) = 0. Conversely, if the QFT of a μ2 -complex function is perplex, i.e. Fs (ω) = 0, then f (−t) = −f (t), so that f (t) is odd. as:

What we have shown is that if f (t) is μ2 -complex and split into even and odd parts f (t) = fe (t) + fo (t)

[3.33]

where fe (t) is even with respect to time and fo (t) is odd, then: F (ω) = Fse (ω) + Fpo (ω)μ

[3.34]

where Fse (ω) = F R {fe } and Fpo (ω)μ = F R {fo } .

[3.35]

To switch to left-sided QFT’s one instead starts with an alternate simplex-perplex decomposition, namely, let f (t) = fs (t) + μfp (t) where the μ factor is on the left of fp (t) instead of on the right 4. The above analysis mimics the symmetry properties of the complex Fourier transform, but we speak of μ2 -complex even and odd functions instead of real even and odd functions of time. 4 Note this fp is the conjugate of fp from [3.21] by the swap rule.

Quaternion Fourier Transforms

45

It is desirable that the complex Fourier transform properties be contained within the QFT as a special case. Indeed this is so, for if we restrict the function to the μcomplex subﬁeld of H, then the classic complex properties carry over, namely, if f (t) is strictly real-valued then: F (ω) = F (−ω) and if f (t) is purely μ-imaginary, i.e. f (t) ∈ μR, F (ω) = −F (−ω). Likewise, if f (t) is even, symmetric, i.e. f (t) = f (−t), then F (ω) ∈ R, but if f (t) is odd, anti-symmetric f (t) = −f (−t), then F (ω) ∈ μR. Proof of this is simple; since both f (t) and the kernel of the transform eμωt are from the same complex subﬁeld, Cμ = R + μR, they are isomorphic to the standard complex numbers. 3.1.3.2. Orthogonal 2D planes split Every q ∈ H can be split as q = q+ + q− with q± = 12 (q ± μ1 qμ2 ) as in [1.58] where μ1 μ2 = μ3 [HIT 07]. In terms of components, if q = r0 +μ1 r1 +μ2 r2 +μ3 r3 , then q± = {(r0 ± r3 ) + μ1 (r1 ± r2 )}

1 ± μ3 1 ± μ3 = {(r0 ± r3 ) − μ2 (r1 ± r2 )} . 2 2

3.1.4. Inter-relationships between deﬁnitions The left- and right-handed transforms given in [3.5] and [3.7] are related, as we might expect, and the operations that play a role in these relationships are: conjugation, time/frequency reversal and symplectic decomposition. The operation with the simplest effect to consider is conjugation. In the classical complex case (when complex signals are being considered), conjugation and time/frequency reversal are connected: conjugation in one domain corresponds to conjugation and reversal in the other. For example, F {f (t)} = F{f (−t)}. In the quaternion case, however, there is a difference, and the reason is the different behavior of the quaternion conjugate compared to the complex conjugate, when applied to a product. In the complex case, the conjugate of a product is equal to the product of the conjugates: (z1 z2 ) = z1 z2 . In the quaternion case, however, order must be reversed: q1 q2 = q2 q1 . The difference between the left- and right-sided

46

Quaternion Fourier Transforms for Signal and Image Processing

transforms is the order of the signal and exponential; hence, it is not surprising that change of handedness and conjugation are connected. Figure 3.1 shows how transforms of opposite handedness are related with respect to conjugation. Algebraically, the relationships can be expressed as: R L {f } F±μ f = F∓μ

and

L R {f }. F±μ f = F∓μ

[3.36]

(For simplicity, the ﬁgure omits the permutations of sign in the transform axis μ). F R (ω) ✉ FR ✒ f (t) ✉ ❅ ❅ ❅ ❘ F L❅ ❅

F R (ω) ✉ ❅ ❅ L ❘F ❅ ❅ ❅

❅✉f (t)

✒F R ❅✉ L

F (ω)

✉ F L (ω)

Figure 3.1. Graph showing conjugation relationships between left- and right-handed one-sided QFTs. The dotted lines represent conjugation. The solid directed lines represent transforms

Time/frequency reversal is slightly more complicated5. In the classical complex case, changing the sign of the time variable (time reversal) corresponds to changing the sign of the frequency variable (frequency reversal). Thus, if we have: F{f (t)} = F (ω), then: F{f (−t)} = F (−ω). In the quaternion case, however, the symplectic decomposition plays a part. Figure 3.2 shows how transforms of opposite handedness are related by time/frequency reversal and symplectic decomposition. Note that application of the same transform twice to a signal f (t) results in a timereversed signal f (−t) (ignoring scale factors). It is immaterial whether the transform chosen is left- or right-handed. The relationship between the two transforms is shown 5 Computationally time/frequency reversal is not trivial and we discuss it in section 3.3.3.2.

Quaternion Fourier Transforms

47

by the vertical dashed line that represents a partial frequency reversal. The symplectic decomposition of the right-handed QFT was given by [3.22] as: F R (ω) = Fs (ω) + Fp (ω)μ,

[3.37]

where Fs (ω), Fp (ω) ∈ Cμ2 . The left-handed QFT using the same symplectic decomposition is found to be: F L (ω) = Fs (ω) + Fp (−ω)μ,

[3.38]

so that the simplex parts Fs (ω) are the same in both cases but the perplex parts Fp (±ω) of the two transforms are reversed in frequency. We have not considered the dual case, where the frequency-domain representation is reversed in frequency, but it should be obvious that the effect in the time domain would be a partial time reversal with only the perplex component reversed.

FR ✒ f (t)

✉ ❅

See [3.37] F R (ω) ✉ ❅ ❅ R ❅ ❘F ❅ ❅ ❅✉f (−t)

❅ ❘ ❅ F L❅ ❅ ❅✉ F L (ω)

✒F L See [3.38]

Figure 3.2. Graph showing time/frequency reversal relationships between leftand right-handed one-sided QFTs. The dotted lines represent time/frequency reversal (horizontally) or partial time/frequency reversal (vertically). The solid directed lines represent transforms

3.1.5. Convolution and correlation theorems Theorems for convolution and correlation as shown in Table 3.1 for the complex Fourier transform are possible in the quaternion case, but the formulas are different for each formulation of transform. The results in the previous section show that leftand right-sided transforms are closely related, and conjugation plays a key role because of the interaction between conjugation and the re-ordering of a quaternion

48

Quaternion Fourier Transforms for Signal and Image Processing

product, which are steps that often occur in the derivation of formulas involving quaternion products. Conjugation of a quaternion exponential requires a change of sign in the root of −1, and of course this corresponds to changing a forward transform into a reverse or vice versa. For these reasons, it is more sensible to construct such formulas using a mix of transforms, rather than conﬁning the formula to a single transform, because this can yield a simpler formula. A secondary consideration is that when the forward and inverse transforms have different scale factors (as they do in some numerical implementations), combining forward and inverse transforms of the same handedness becomes awkward because of the different scale factors: it is simpler to combine transforms of the same direction (forward and inverse) of different handedness, in order to avoid such difﬁculties. We show how a correlation formula may be derived in section 4.2.2. 3.2. 2D quaternion Fourier transforms The extension of the 1D QFT to two dimensions gives rise to a multitude of possibilities. The obvious choice of altering the transform kernel from e−μωt to e−μ(ωx+νy) for functions f (x, y) : R → H provides the same four left-right and forward-reverse transforms, but there are other choices. One concept left unexplored is the implication of the exponential product of two quaternions, i.e. in general ep eq = ep+q , as noted in [1.38]. Hence, e−(μ1 ωx+μ2 νy) = e−μ1 ωx e−μ2 νy . As soon as the exponent is split like this, we have the ability to “sandwich” the function f (x, y) between the exponents. We also have the ability to use a different root of −1 for each exponent. Table 3.3 details these eight forward forms. Doubling this number to account for the reverse transforms yields a total of 16 different QFTs. Only a limited number of most commonly used forms will be explored in this book.

Single-axis Dual-axis Factored

Left e−μ(ωx+νy) f (.) e−(μ1 ωx+μ2 νy) f (.) e−μ1 ωx e−μ2 νy f (.)

Right f (.) e−μ(ωx+νy) f (.) e−(μ1 ωx+μ2 νy) f (.) e−μ1 ωx e−μ2 νy

Sandwich e−μωx f (.) e−μνy — −μ1 ωx e f (.) e−μ2 νy

Table 3.3. Forward transform kernel deﬁnitions for f : R2 → H

3.2.1. Deﬁnitions The 2D single-axis, right-sided QFTs (forward and reverse) are deﬁned as: R F∓μ {f } (ω, ν) = F R (ω, ν) =

f (x, y) e∓μ(ωx+νy) dx dy, R2

[3.39]

Quaternion Fourier Transforms

49

with their inverses being formed with the conjugate of the exponential kernel as: −R F∓μ F R (x, y) = f (x, y) =

F R (ω, ν) e±μ(ωx+νy) dω dν,

[3.40]

R2

The 2D single-axis, left-sided QFTs (forward and reverse) are deﬁned as: L F∓μ {f } (ω, ν) = F L (ω, ν) = κ2∓

e∓μ(ωx+νy) f (x, y) dx dy,

[3.41]

R2

and their inverses are derived from the conjugate of the kernels as: −L F∓μ F L (x, y) = f (x, y) = κ2±

e±μ(ωx+νy) F L (ω, ν) dω dν.

[3.42]

R2

The 2D factored, sandwich QFTs (forward and reverse) are deﬁned as: S F∓(μ {f } (ω, ν) = F L (ω, ν) = κ2∓ 1 ,μ2 )

e∓μ1 ωx f (x, y) e∓μ2 νy dx dy,

[3.43]

R2

and their inverses are derived from the conjugate of the kernels as: −S F∓(μ

1 ,μ2 )

F L (x, y) = f (x, y) = κ2±

e±μ1 ωx F S (ω, ν) e±μ2 νy dω dν. [3.44] R2

The remaining ﬁve transform deﬁnitions follow a similar pattern, with kernel placements as shown in Table 3.3. Discovering which 2D quaternion-valued functions can be Fourier-transformed follows the same logic of the 1D case in section 3.1.1; we can borrow from the classic complex case. Split f (x, y) via symplectic decomposition into simplex and perplex parts with respect to μ2 as: f (x, y) = fs (x, y) + fp (x, y) μ2 , where fs , fp ∈ Cμ ,

50

Quaternion Fourier Transforms for Signal and Image Processing

so that the right-sided QFT of [3.39] becomes: F R (ω, ν) = =

R2

R2

fs (x, y) e−μ(ωx+νy) dx dy + fs (x, y) e−μ(ωx+νy) dx dy +

R2

R2

fp (x, y) μ2 e−μ(ωx+νy) dx dy, fp (x, y) eμ(ωx+νy) dx dy μ2

= Fs (ω, ν) + Fp (−ω, −ν)μ2 , In the same manner, [3.41] yields for the left-sided QFT: F L (ω, ν) =

R2

e−μ(ωx+νy) fs (x, y) dx dy +

R2

e−μ(ωx+νy) fp (x, y) dx dy μ2

= Fs (ω, ν) + Fp (ω, ν)μ2 From the fact that the left- and right-sided 2D QFTs can be decomposed into a sum of complex subﬁeld functions, fs,p (x, y), the existence and invertibility of these 2D QFTs is inherited from the 2D complex case. This is true even for the sandwich form of the QFT of [3.43], but showing this takes more work. Split f (x, y) via symplectic decomposition into simplex and perplex parts with respect to μ2 as: f (x, y) = fs (x, y) + fp (x, y) μ2 , where fs , fp ∈ Cμ1 , so that [3.43] becomes: F (ω, ν) = Fs (ω, ν) + Fp (ω, ν) μ2

[3.45]

where Fs (ω, ν) = Fp (ω, ν) =

R2

R2

e−μ1 ωx fs (x, y)e−μ2 νy dx dy, and e−μ1 ωx fp (x, y)e−μ2 νy dx dy.

The derivation begins with a 2D QFT equivalent to a 2D complex Fourier transform. Let: Fs,μ1 (ω, ν) =

R2

e−μ1 ωx fs (x, y)e−μ1 νy dx dy

[3.46]

Quaternion Fourier Transforms

51

so that Fs,μ1 (ω, ν) ∈ Cμ1 because both kernels of the QFT are the same. This is different from the deﬁnition of [3.44] on purpose. Expand the following expression, using Euler’s cosine equation [1.51] in the last step, as: 1 2

Fs,μ1 (ω, ν) + Fs,μ1 (ω, −ν) =

1 2

=

1 2

=

R2

R2

R2

e−μ1 ωx fs (x, y)e−μ1 νy dx dy +

1 2

R2

e−μ1 ωx fs (x, y)e+μ1 νy dx dy

e−μ1 ωx fs (x, y) e−μ1 νy + e+μ1 νy dx dy

e−μ1 ωx fs (x, y) cos(νy) dx dy

Likewise, using the sine function of [1.51], another expression: 1 2

Fs,μ1 (ω, ν) − Fs,μ1 (ω, −ν) =

1 2

=

1 2

R2

R2

=

R2

e−μ1 ωx fs (x, y)e−μ1 νy dx dy −

1 2

R2

e−μ1 ωx fs (x, y)e+μ1 νy dx dy

e−μ1 ωx fs (x, y) e−μ1 νy − e+μ1 νy dx dy e−μ1 ωx fs (x, y) sin(νy) dx dy {−μ1 }

Combining these two expressions, with the second post-multiplied by −μ3 (= μ2 μ1 ), yields: 1 2

Fs,μ1 (ω, ν) + Fs,μ1 (ω, −ν) + =

R2

1 2

Fs,μ1 (ω, ν) − Fs,μ1 (ω, −ν) {−μ3 }

e−μ1 ωx fs (x, y)e−μ2 νy dx dy = Fs (ω, ν)

which is the simplex part of the 2D QFT we are looking for from [3.45]. So that, after simplifying the combined expressions, Fs (ω, ν) = Fs,μ1 (ω, ν)

1 − μ3 1 + μ3 + Fs,μ1 (ω, −ν) 2 2

[3.47]

Substituting fp (x, y) for fs (x, y) in [3.46] provides similar expressions for the perplex part of the 2D QFT, namely: Fp (ω, ν) = Fp,μ1 (ω, ν)

1 − μ3 1 + μ3 + Fp,μ1 (ω, −ν) 2 2

[3.48]

52

Quaternion Fourier Transforms for Signal and Image Processing

which gives us the desired result: F (ω, ν) = =

Fs (ω, ν) + Fp (ω, ν) μ2 1 − μ3 1 + μ3 + Fs,μ1 (ω, −ν) 2 2 μ + μ1 μ − μ1 + Fp,μ1 (ω, ν) 2 + Fp,μ1 (ω, −ν) 2 . 2 2 Fs,μ1 (ω, ν)

3.2.2. Basic transform pairs Following is a list of theorems that can be proved using the deﬁning equations of section 3.2.1; it is assumed that all functions considered actually have Fourier transforms. For the sake of brevity, we consider only the forward versions of each form; the reverse form pairs will be obvious. Table 3.4 summarizes the results that follow. Linearity: of the transform integrals implies: αf1 (x, y) + βf2 (x, y)

αF1 (ω, ν) + βF2 (ω, ν),

[3.49]

where α, β ∈ R and fn (x, y) Fn (ω, ν). Hence, all eight 2D QFT deﬁnitions shown in Table 3.3 are linear operators. Scaling: Consider the function f (αx, βy), which is obtained by replacing (x, y) by (αx, βy), where α, β ∈ R. The QFT of this function is given as: FμS1 ,μ2 {f (αx, βy)} =

1 FS |α| |β|

ω ν , α β

.

[3.50]

Change of Coordinates: 2D functions can undergo a change of coordinates in the function’s domain, or index space. In general, a change of coordinates is a mapping: (x, y) → (αx + βy, γx + δy) where α, β, γ, δ ∈ R and αδ − βγ = 0. So, consider the function f (αx + βy, γx + δy), which is obtained by a change of coordinates of the function f (x, y). The QFT of this function is given as: FμS1 ,μ2 {f (αx + βy, γx + δy)} = 1 αδ − βγ

1 + μ3 ˆ + 1 − μ3 F (ω δˆ + ν γˆ , ν α ˆ [3.51] F (ω δˆ − ν γˆ , ν α ˆ + ω β) ˆ − ω β) 2 2

where μ3 = μ1 μ2 , α ˆ=

α ˆ αδ−βγ , β

=

β ˆ αδ−βγ , γ

=

γ ˆ αδ−βγ , δ

=

δ αδ−βγ .

Quaternion Fourier Transforms

53

Shifts: If a function is shifted in both directions by a constant (δx , δy ), the spectral amplitude remains the same but a linear term is added to the respective “phase” of each axis as: FμS1 ,μ2 {f (x − δx , y − δy )} = e−μ1 ωδx F S (ω, ν)e−μ2 νδy .

[3.52]

Derivatives: Taking the nth partial derivative of a function with respect to x, we conclude that: FμS1 ,μ2

∂n f (x, y) ∂xn

= (μ1 ω)n F S (ω, ν).

[3.53]

Likewise, the mth partial derivative with respect to y gives: FμS1 ,μ2

∂n f (x, y) ∂y n

= F S (ω, ν) (μ2 ν)n

[3.54]

and together they yield: FμS1 ,μ2

∂ n+m f (x, y) ∂y n ∂y m

= (μ1 ω)n F S (ω, ν) (μ2 ν)m .

[3.55]

These formulas for the partial derivatives do not guarantee that the transforms exist; but if they do, the formulas will be as given. Impulse functions: In applications we use functions that are singular, e.g., the Dirac delta function. To obtain the following results, we use this deﬁnition for the delta function: with f (x, y) an arbitrary function, continuous at a given point (xo , yo ), δ(x, y) is such that:

R2

δ(x − xo , y − yo )f (x, y) dx dy = f (xo , yo ).

It is assumed here that δ(x, y) is a real-valued function so as to commute with f (x, y). The QFT of the 2D Dirac delta function is identical for all QFT forms, namely, δ(x, y) 2π

1

and

1

2π δ(ω, ν).

[3.56]

However, the shifted delta function is only common across all factored and singleaxis 2D QFT forms as: δ(x − xo , y − yo ) 2π

e−μ1 ωxo e−μ2 νyo

54

Quaternion Fourier Transforms for Signal and Image Processing

and eμ1 ωo x eμ2 ν o y

2π δ(ω − ωo , ν − ν o ).

Whereas for the dual-axis (un-factored) form: δ(x − xo , y − yo ) 2π

e−(μ1 ωxo +μ2 νyo )

and eμ1 ωo x+μ2 ν o y

2π δ(ω − ωo , ν − ν o ).

3.2.3. Decompositions Any 2D signal can be split into two different ways with regard to symmetry: either into fully even and odd components along the x- and y-axes as6: f (x, y) = fee (x, y) + foe (x, y) + feo (x, y) + foo (x, y), where feo denotes the part of f that is even with respect to x and odd with respect to y, and so forth, or the signal may be split into partial even and odd components as f (x, y) = fe (x, y) + fo (x, y)

[3.57]

= fe (x, y) + fo (x, y)

[3.58]

where fe denotes the part of f that is even with respect to x and no symmetry is assumed with respect to y, and fo denotes the part of f that is odd with respect to y and no symmetry is assumed with respect to x, etc. The partial even-odd split is related to the full even-odd split as: fe (x, y) = 12 {f (x, y) + f (−x, y)} = fee (x, y) + feo (x, y), fo (x, y) = 12 {f (x, y) − f (−x, y)} = foe (x, y) + foo (x, y), fe (x, y) = 12 {f (x, y) + f (x, −y)} = fee (x, y) + feo (x, y), fo (x, y) = 12 {f (x, y) − f (x, −y)} = foe (x, y) + foo (x, y), 6 Recall that a function is even along the x-axis if f (x) = f (−x) and odd if f (x) = −f (−x). Hence, an arbitrary signal can be split using: fe (x) = 12 {f (x) + f (−x)} and fo (x) = 1 {f (x) − f (−x)}, so that f (x) = fe (x) + fo (x). 2

Quaternion Fourier Transforms

55

while the full even-odd split is given by: fee (x, y) = 14 {f (x, y) + f (−x, y) + f (x, −y) + f (−x, −y)}, feo (x, y) = 14 {f (x, y) + f (−x, y) − f (x, −y) − f (−x, −y)}, foe (x, y) = 14 {f (x, y) − f (−x, y) + f (x, −y) − f (−x, −y)}, foo (x, y) = 14 {f (x, y) − f (−x, y) − f (x, −y) + f (−x, −y)}. Using the 2D sandwich QFT [3.43], it can be shown that: S S S S F S (ω, ν) = Fee (ω, ν) − μ1 Foe (ω, ν) − Feo (ω, ν)μ2 + μ1 Foo (ω, ν)μ2

where S S Fee (ω, ν) = F−(μ {fee } (ω, ν) , 1 ,μ2 ) S S Foe (ω, ν) = F−(μ {foe } (ω, ν) , 1 ,μ2 ) S S Feo (ω, ν) = F−(μ {feo } (ω, ν) , 1 ,μ2 ) S S Foo (ω, ν) = F−(μ {foo } (ω, ν) . 1 ,μ2 )

3.2.4. Inter-relationships between deﬁnitions Using the symplectic decomposition of f (x, y) with respect to μ1 , i.e. f (x, y) = fs (x, y) + fp (x, y)μ2 , where fs,p ∈ Cμ1 the right-sided QFT is given by: F R (ω, ν) = Fs (ω, ν) + Fp (ω, ν)μ2 ,

[3.59]

where Fs (ω, ν), Fp (ω, ν) ∈ Cμ1 . The left-sided QFT, using the same symplectic decomposition for f (x, y), is found to be: F L (ω, ν) = Fs (ω, ν) + Fp (−ω, −ν)μ2 ,

[3.60]

so the simplex parts of the left-side and right-side QFTs are the same, while the perplex parts of the two transforms are reversed in frequency.

56

Quaternion Fourier Transforms for Signal and Image Processing 1 2π

R2

eμ1 ωx F (ω, ν)eμ2 νy dω dν

Spatial-domain Spatial reversal

f (x, y) f (−x, y) f (x, −y)

1 2π

R2

e−μ1 ωx f (x, y)e−μ2 νy dx dy

F (ω, ν) F (−ω, ν) F (ω, −ν)

Frequency-domain Frequency reversal

f (x − xo , y − yo ) e f (x, y) eμ2 ν o y Spatial scale f (αx, βy) 1 x y f( α , β) |α| |β| Coordinate change f (αx + βy, γx + δy)

e−μ1 ωxo F (ω, ν) e−μ2 νyo F (ω − ωo , ν − ν o ) Frequency shift 1 F(ω , ν) |α| |β| α β F (αω, βν) Frequency scale 1+μ3 1 ˆ ˆ { 2 F (ω δ − ν γˆ , ν α ˆ + ω β) αδ−βγ 1−μ3 ˆ ˆ + 2 F (ω δ + ν γˆ , ν α ˆ − ω β)}

Linearity

αF1 (ω, ν) + βF2 (ω, ν)

Spatial shift

μ1 ωo x

αf1 (x, y) + βf2 (x, y)

Delta function Shifted delta

δ(x, y) 2π δ(x − xo , y − yo ) 2π eμ1 ωo x eμ2 ν o y 1

nth derivative

∂n f (x, y) ∂xn ∂n f (x, y) ∂y n

∂ n+m f (x, y) ∂y n ∂y m

1 e−μ1 ωxo e−μ2 νyo 2π δ(ω − ωo , ν − ν o ) Shifted delta 2π δ(ω, ν) Delta function (μ1 ω)n F (ω, ν) F (ω, ν) (μ2 ν)n

(μ1 ω)n F (ω, ν) (μ2 ν)m

N OTE.– f (t) : R2 → H, α, β, γ, δ ∈ R, μ3 = μ1 μ2 β γ α δ α ˆ = αδ−βγ , βˆ = αδ−βγ , γˆ = αδ−βγ , δˆ = αδ−βγ . Table 3.4. 2D factored, sandwich quaternion Fourier transform pairs

If the signal f (x, y) is split into the partial even and odd parts with respect to the y-axis, i.e. f (x, y) = fe (x, y) + fo (x, y), then we ﬁnd: F S (ω, ν) = FeL (ω, ν) − μ1 FoL (ω, ν)μ2

[3.61]

where L FeL (ω, ν) = F−μ {fe } (ω, ν) L FoL (ω, ν) = F−μ {fo } (ω, ν)

so that the sandwich form can be constructed from the left-handed form. See [YEH 08] for details of the proof. Likewise, if the same signal f (x, y) is split into the partial even and odd parts with respect to the x-axis, i.e., f (x, y) = fe (x, y) + fo (x, y), then we ﬁnd: F S (ω, ν) = FeR (ω, ν) − μ1 FoR (ω, ν)μ2

[3.62]

Quaternion Fourier Transforms

57

where R FeR (ω, ν) = F−μ {fe } (ω, ν) R FoR (ω, ν) = F−μ {fo } (ω, ν)

so that the sandwich form can also be constructed from the right-handed form. 3.3. Computational aspects In this section, we consider the practical and numerical issues that arise in the computation of QFTs, and the process of veriﬁcation that the numerical results computed by the transform are in agreement with the mathematical deﬁnition. Veriﬁcation is also applicable to formulas such as those presented in section 3.1.4, and we discuss the details of time/frequency reversal, because this plays a key role in the inter-relationships between transforms. 3.3.1. Coding Coding of a QFT requires at least some minimal quaternion arithmetic operations to be implemented. These operations are normally provided by a software library, and unless a library already exists for the chosen computing platform or language, it will be necessary to create one. Examples of libraries are the Boost library for C++, which contains some quaternion functions in the Math toolkit [BRI 13], and the QTFM library for MATLAB® [SAN 13b], which implements a wide range of quaternion functions using as far as possible standard MATLAB® notation and conventions. Both of these libraries are open source; so implementation of a quaternion library for another language or platform should not be difﬁcult, since the code of the existing libraries can be inspected for ideas. Even if a library exists, there may be bugs in it that render the results of a quaternion transform computation incorrect under some circumstances, and this fact should be considered when testing an implementation to verify correctness. A further use for a quaternion library is the analysis of the results from computing a transform, whether this is some form of processing (e.g. ﬁltering in the frequency domain) or interpretation (e.g. by plotting the Fourier coefﬁcients or some aspects of them such as their magnitude, axis or angle). In the following sections, we consider four possible ways to compute a QFT numerically: 1) a direct quaternion computation using a discrete QFT implementing directly the desired transform as a single or double summation;

58

Quaternion Fourier Transforms for Signal and Image Processing

2) a direct quaternion FFT algorithm using decimation in time or frequency or another method [BRA 00]; 3) a decomposition into a pair of complex Fourier transforms that can be computed using existing complex FFT library routines; 4) a matrix method that avoids the need for a quaternion library, and thus provides a way to validate a computation based on a quaternion library. Of these methods, the third is probably the best choice for most purposes, since it requires minimal effort to code, and exploits the signiﬁcant work that others have invested in the fast and accurate computation of complex FFTs. We discuss each of the above methods in turn in the following sections. 3.3.1.1. Direct computation of a DFT Any discrete Fourier transform (DFT) (complex, quaternion or otherwise) can be computed directly from the deﬁnition, using a summation (or double summation in the case of a 2D transform). This is often referred to in the literature as a DFT, the word discrete referring to the discretization of time (or position in an image). For example, a 1D left-sided QFT pair as deﬁned in [3.7] could be computed in discrete form using the following pair of summations7: N −1

F [u] =

exp −μ 2π n=0

1 f [n] = N

N −1

nu f [n] N

nu exp +μ 2π F [u] N u=0

[3.63]

Now obviously the conversion from the integral notation to the discrete notation assumes another difference: the integral is over all time, whereas the summation is over a limited range of sample indices. This is a matter that is standard in digital signal and image processing, where the integral notation is often used, even though a discretized form is actually used for numerical computation. The summation is simple enough to be coded in a few lines of code in a high-level language. In MATLAB®, it is possible to avoid a loop and perform a vectorized computation by precomputing the exponential for all index values8, then multiply the 7 Exactly how the scale factor is handled is a matter of choice. It could be applied to the forward transform, or the square root of the scale factor could be applied in both directions. 8 This assumes there is enough memory available to store the entire set of exponential values. For longer transforms, this is not a reasonable assumption, and a more sophisticated computation is recommended in order to store the exponential values more economically.

Quaternion Fourier Transforms

59

resulting matrix by the signal vector to obtain the result. This is how the QTFM function qdft.m is coded, and the following code extract shows the core of the code: F = E * f; The variable E here contains a matrix of dimension N × N containing the values of the exponential function for all possible index values n and u in [3.63], as discussed, for example, by Golub and van Loan [GOL 96, section 4.6.4]. A single matrix product then multiplies the signal by this matrix to yield the frequency-domain coefﬁcients. This works because the quaternion operations for element-wise and matrix multiplication are implemented within the toolbox, and the exponential function is vectorized, just like the standard MATLAB® exponential function. (Of course the classical complex DFT can be computed by exactly the same code, by replacing mu with 1i.) A major disadvantage of this method, whether implemented using loops, or as vectorized code, is that the run-time scales as the square of the length of the signal array, N . This can be easily seen in Figure 3.3, which was obtained by running the qdft function on random quaternion arrays from length 1 to 2,000, and averaging over three measurements. For comparison, it also shows the results from the qfft function (see section 3.3.1.2). 0.3

qdft qfft

0.25

Run−time, seconds

0.2

0.15

0.1

0.05

200

400

600

800 1000 1200 Transform length, N

1400

1600

1800

Figure 3.3. Run-time for the qdft and qfft functions of the QTFM toolbox [SAN 13b] on random quaternion signals

60

Quaternion Fourier Transforms for Signal and Image Processing

As can be seen, even for a short array of only 2,000 samples, the run-time for the qdft function is around a third of a second, whereas the qfft computes in approximately 10 ms, with a very much slower increase in run-time with increasing N . The corresponding transform for a 2D array has run-time that scales as N 4 , assuming a square array of dimension N × N . (None of this is unique to QFTs, but since computing with quaternions takes four times as much memory and 16 times as much time as computing with reals, run-time is an important issue.) Despite the disadvantage of long run-times, and run-times that scale with the square of the signal length, a direct coding of a quaternion DFT is useful as a reference against which the results of faster methods can be compared. As can be seen, the translation of mathematics into code can be so direct that there is little room for error. However, as already explained, it is not difﬁcult to verify the results of a quaternion DFT against a complex Fourier transform implementation. 3.3.1.2. Direct coding of an FFT The fast Fourier transform (FFT) is a family of algorithms for computing the DFT with much reduced run-time scaling and better accuracy. Both advantages arise from a signiﬁcant reduction in the number of arithmetic operations needed. We will not discuss here how this reduction in run-time is achieved, since it is not pertinent to this book, and is well-covered elsewhere (see, for example, [BRA 00]). A DFT takes time proportional to N 2 , whereas an FFT has run-time proportional to N log N . The difference is clearly illustrated in Figure 3.3 where the second trace shows the run-time of the qfft function (close to the horizontal axis). Direct coding of a quaternion FFT algorithm is possible, but it requires signiﬁcant effort, even for a power-of-two FFT (the easiest type of algorithm to implement), and it does not seem worth attempting nowadays. It was used by Sangwine in 1996 [SAN 96] because there was no other fast method capable of working with 512×512 color images, but it required very signiﬁcant effort to code and test, even after a quaternion arithmetic library had been written (in Ada 95). The implementation was limited to array sizes that were powers of 2 because handling other sizes was much more difﬁcult, and it was not necessary to handle other sizes in order to demonstrate the validity of using QFTs on color images. Direct coding of quaternion FFTs was effectively made obsolete by the idea presented in the next section, which is the method recommended for anyone wishing to code a quaternion FFT. 3.3.1.3. Computation using decompositions and complex FFTs In 2000, two of the authors published a paper showing how to decompose a QFT into a pair of complex Fourier transforms, and thus compute a fast quaternion transform without needing to code an explicit quaternion FFT [ELL 00a]. The method is also explained in [ELL 07c, section IV.B]. This technique has a major advantage over an explicit quaternion FFT: it exploits the capabilities of any existing

Quaternion Fourier Transforms

61

fast complex FFT code, and any future enhancements made to the code (subject to recompilation or relinking with the newer version, of course). The most obvious candidate is the excellent FFTW library written at MIT by Matteo Frigo and Steven G. Johnson [FRI 05]9. FFTW is written in C, and it handles arbitrarily sized arrays of data, even with large prime-number lengths. For repeated calculation on multiple data of the same length, it may be ﬁne-tuned to optimize the speed on a particular machine architecture and at a particular address in memory (the latter affects cache hits and misses, which impacts performance, and can be taken into account in the tuning process). The mathematical basis of the decomposition method is described in section 1.4.1.4 where it is shown that any quaternion can be decomposed into two quaternions parallel and perpendicular to a chosen unit pure quaternion. Applying this idea to a particular QFT formulation, as in section 3.1.1 (see [3.9] and [3.11]), it is possible to decompose the signal and the transform into a pair of quaternion transforms each of which is isomorphic to a complex Fourier transform. Note that it is necessary to ﬁnd the decomposition formula for each particular transform. After computing the pair of complex transforms, the complex results are used to construct isomorphic quaternions, which are then combined into the quaternion result. Numerically, a simpler implementation of this idea is to compute a change of basis of the vector part of the quaternion signal, and this is how the technique is actually implemented in the QTFM library [SAN 13b, Function change_basis]. Given a basis, for example, {1, μ1 , μ2 , μ1 μ2 } and a quaternion q = w + xi + yj + zk expressed in the standard basis {1, i, j, k}, the four numerical components of the quaternion expressed in the new basis are simply computed as: w =w x = q, μ1 y = q, μ2 z = q, μ3 After the change of basis, the numerical components of the quaternion signal can be directly utilized as complex signals z1 = w + x I and z2 = y + z I, and the complex transforms computed. Finally, the change of basis has to be inverted after reassembling the complex results into a quaternion result. Inversion requires that the standard basis be expressed in terms of the new basis. This is straightforward because a basis can be expressed as an orthogonal 3 × 3 matrix, which is trivially invertible (note that the change of basis does not affect the scalar part). Therefore, the three unit quaternions needed to invert the change of basis are simply computed from μ1 , μ2 and 9 FFTW is the library used by MATLAB® to compute FFTs. It is available at: http://www.fftw.org.

62

Quaternion Fourier Transforms for Signal and Image Processing

μ1 μ2 by constructing new unit pure quaternions from their x, y and z components, in the same manner as constructing new rows from the columns of the corresponding 3 × 3 matrix. 3.3.2. Veriﬁcation The transforms presented in this book are not just theoretical results. They have been implemented in computer code and applied to problems in signal and image processing. This requires numerical implementations, and numerical implementations should be veriﬁed carefully. In the experience of the authors, numerical implementations also provide an excellent way to verify the mathematics, particularly with hypercomplex formulas where the interactions between non-commutativity and operations such as conjugation provide plentiful sources of error in mathematical derivations. 3.3.3. Veriﬁcation of transforms Here, we consider the problem of verifying that the numerical results produced by a transform implementation are in agreement with the mathematical deﬁnition. This is not a trivial issue: it is simple to check that an implementation of a particular transform inverts correctly – compute the forward transform of a randomly created signal or image f , giving F , and then compute the inverse transform of F , giving g. Then, f and g should be identical to within a suitably small tolerance (computing their difference is the simplest way to check this, followed by ﬁnding the maximum absolute value of the result). However, this process does not verify that the transform computes the Fourier coefﬁcients correctly, according to the mathematical deﬁnition, although it is a useful initial check – if the test described fails (the difference between f and g is not small), then there is obviously an error in the coding of the transform. A most likely error is that f and g could differ signiﬁcantly by a scale factor, which is not serious, and easily corrected. This problem is easily diagnosed by computing the quaternion ratio of the elements in f to those in g (multiply each sample in f by the quaternion inverse of the corresponding sample in g). A scale factor error will manifest as a constant real-valued result for all elements of f and g. If the error is not a scale factor error, the problem of diagnosis is more difﬁcult. This is where it is valuable to code a DFT implementation, as discussed in section 3.3.1.2, because this provides a simple, directly coded implementation that is useful as a reference for a more complex implementation based on decompositions and complex FFTs. However, even a DFT implementation requires a quaternion library, and there could be errors in a library that affect both transform implementations. Given the limited number of hypercomplex libraries in existence, running on different platforms, and written in different languages, even checking one against another is a non-trivial task. Until 2010, there was no known solution for this problem, other than the use of a

Quaternion Fourier Transforms

63

different library for veriﬁcation, but the formulation presented in the next section even overcame that problem. It is also possible to verify results against a classical complex Fourier transform implementation. The latter is possible because a QFT reduces to the classical complex Fourier transform under certain conditions. If these conditions are satisﬁed, a numerical comparison is possible between a quaternion implementation and a classical complex implementation. The most general way to do this is to construct a quaternion signal so that all the samples will have the same axis, and this axis is the same as the transform axis10. Here is the process: construct a complex test signal f (e.g. using a random number generator), such that f = fr + Ifi . Then construct an isomorphic quaternion signal q = fr + μfi . Compute the complex Fourier transform of the complex signal, F = F {f }, and the QFT of the quaternion signal, with axis μ, Q = F {q}. Then, to within a tolerance, we should ﬁnd that Fr = S(Q) and μFi = V(Q). 3.3.3.1. Transform formulation using a matrix exponential A solution to the problem of verifying a transform without requiring a second quaternion library was ﬁrst published in 2012 [SAN 12], and it is of considerable mathematical interest. It is based on a matrix representation of quaternions (as already discussed in sections I.6.2 and 1.4.3), and it can compute a QFT without a quaternion library. A one-sided quaternion DFT pair such as that given in [3.5] can be represented in the following form, which combines matrix representations of quaternions with quaternions represented as four-element column vectors: N −1

F [u] = S

exp −J 2π

nu f [n] N

[3.64]

exp +J 2π

nu F [u] N

[3.65]

n=0 N −1

f [n] = T u=0

where J is a 4 × 4 matrix root of −1 representing a quaternion root of −1 using the representation of [1.69]; f [n] and F [u] are quaternion-valued discrete-time signals with N samples; each sample, indexed by n or u respectively, being a column vector of four elements corresponding to the four components of the corresponding quaternion signal sample; and the two scale factors S and T multiply to give 1/N in the usual way (the distribution of the scale factors between the forward and inverse transforms being arbitrary). 10 It is assumed here that we are discussing a one-sided transform with a single axis. Obviously, matters get a bit more complicated if the transform is two-sided with multiple axes.

64

Quaternion Fourier Transforms for Signal and Image Processing

This formulation depends on two results that are not obvious at ﬁrst glance. The ﬁrst is that Euler’s formula is valid when the root of −1 is replaced by a matrix that squares to a negated identity matrix (we call this a matrix root of −1): exp Jθ = I cos θ + J sin θ where I is an identity matrix, and J2 = −I. For a proof, we refer the reader to [SAN 12], but note that the result is easily veriﬁed using MATLAB® or other mathematical software. The exponential, of course, must be a matrix exponential [ROW 14]. The second result is due to Ickes [ICK 70], who showed that the matrix representation of a quaternion, as discussed in section I.6.2 and 1.4.3 could be combined with a four-element column vector containing the four components of a quaternion to represent the quaternion product as a matrix-vector product (the result of course is a four-element column vector). Surprisingly, Ickes also showed that the matrix does not have to represent the left operand in the product: by transposing the lower-right 3 × 3 part of the matrix, the matrix-vector product yields the result when the quaternion represented by the matrix is the right operand, despite the fact that the matrix remains the left operand in the matrix-vector product. The transform pair in [3.64] can be implemented in MATLAB® code as follows (note the directness of the coding): N = size (f , 2); F = zeros ( size ( f )); for n = 0: N -1 for u = 0: N -1 F (: , u + 1) = F (: , u + 1) + ... expm ( - J .* 2 .* pi .* n .* u ./ N ) * f (: , n + 1); end end Here, F and f are arrays of four rows and N columns, representing a signal of N quaternion samples; and J is the quaternion root of −1 used in the transform, represented as a 4 × 4 matrix according to [1.69]. This code runs on standard MATLAB®, i.e. without needing any quaternion toolbox or library, and it has been veriﬁed to produce the same numerical results as the quaternion functions qdft and qfft in the QTFM library [SAN 13b]. The approach discussed here has been shown also to work for two-sided transforms such as that in [3.43], using the simple approach of representing all the quaternions (left and right exponentials, and the signal samples) as matrices. However, using the matrix-vector approach of Ickes [ICK 70], it should be possible to formulate even a two-sided transform as a matrix-matrix-vector product, placing the two exponential matrices on the left, but using the middle one to represent the right-side exponential in Ickes’ partially transposed manner.

Quaternion Fourier Transforms

65

Finally, note that the matrix formulation also works in the classical complex case, and for other hypercomplex algebras such as Clifford algebras [SAN 12]. 3.3.3.2. Time and frequency reversal As we have seen in section 3.1.4, the various transform formulations are inter-related, and there are two operations that play a key role in the inter-relationships: conjugation and time/frequency reversal. Conjugation is computationally simple, but it is worth expanding a little on time/frequency reversal because veriﬁcation of relationships involving time/frequency reversal is not as simple as it might appear. In the continuous-time case, time and frequency reversal simply require the time or frequency variable to be negated, and this is mathematically and algebraically trivial. In discrete time, however, some care is needed to implement this operation correctly, whether working algebraically, or with a numerical implementation. To understand this point, it is necessary to understand how the frequency-domain coefﬁcients of a discrete-time Fourier transform are stored.

N even

Swap ❄

N odd DC 0

···

1

N 2

−1

❄

❄ Nyq.

Positive frequencies

DC 0

❄

N 2

Negative frequencies N 2

+1

···

N −1

Swap ❄

❄ Positive frequencies 1

···

❄

❄

Negative frequencies

N−1 N +1 2 2

···

N −1

Figure 3.4. Layout of frequency-domain coefﬁcients for a DFT of N samples

Figure 3.4 shows the layout of the samples in the frequency domain for a 1D DFT (classical or quaternion, the layout is the same). The indexing in the ﬁgure shows zerobased indexing, as used in mathematical formulas. In programming languages with one-based indexing, program code must be adapted to use array indices for accessing the transform elements, but zero-based indexing for calculations. In the even-length case, the coefﬁcients are divided into four categories: – the zero-frequency or DC coefﬁcient, stored in the ﬁrst index position with numerical index zero;

66

Quaternion Fourier Transforms for Signal and Image Processing

– the Nyquist frequency coefﬁcient, stored in the ﬁrst index position of the second half of the array with numerical index N/2; – the positive frequency coefﬁcients, stored between the DC and Nyquist coefﬁcients, with the lowest frequency on the left, with index 1, and the highest frequency on the right, with index (N/2) − 1; – the negative frequency coefﬁcients, stored between the Nyquist coefﬁcient and the end of the array, with the lowest frequency coefﬁcients (i.e. closest to zero frequency) at the end, with index N − 1, and the highest frequency next to the Nyquist coefﬁcient, with coefﬁcient (N/2) + 1. In the odd-length case, the Nyquist frequency coefﬁcient is not present. The zerofrequency coefﬁcient is at index zero, the positive frequency coefﬁcients run from index 1 to (N − 1)/2, and the negative frequency coefﬁcients run from index N − 1 down to (N + 1)/2. Time reversal of a signal corresponds to the interchanging of the positive and negative frequency coefﬁcients. Classically, in digital signal processing texts concerned with real signals, time reversal corresponds to conjugation in the frequency domain (because the Fourier coefﬁcients of a real signal have complex conjugate symmetry – the negative frequency coefﬁcients are the conjugates of the corresponding positive frequency coefﬁcients). However, as soon as we move from real signals to complex, it is much simpler to think in terms of exchanging positive and negative frequency coefﬁcients, since this is fundamentally how the correspondence between the two domains operates. The swaps indicated by the arrows in Figure 3.4 show how to achieve time reversal of a signal in either domain. In practice, the swaps can be very simply achieved by reversing the order of the whole array, excluding the DC coefﬁcient, which should stay where it is in index position zero. This operation leaves the Nyquist coefﬁcient (in the even case), unmoved. Alternatively, reversing the whole array, and then applying a right circular shift by one position, will put the DC coefﬁcient back in index position zero. To understand why it is that the same operation works in both domains, note that in time, the sample at index zero must remain in the same position, since when negating the time variable, t = 0 remains at 0. Second, we must remember that the DFT treats the signal as periodically extended to inﬁnity. Thus, the sample stored at the rightmost end of the array is also logically adjacent to the sample with index zero, but on its left. In two dimensions, the above condition applies to both the rows and columns of an image. The simplest way to consider this is to imagine the reversal applied to the rows ﬁrst, then to the columns of the result, and indeed this is a simple way to implement the operation numerically.

4 Signal and Image Processing

This chapter presents some examples of applications of quaternion Fourier transforms (QFTs) to process color images and complex-valued signals. The concepts presented in this chapter will be illustrated on simulated and real images and signals.

4.1. Generalized convolution One of the early motivations for the study of QFTs is that they can be used to describe color image ﬁlters when the image pixels are treated as vectors and thus offer freedom to process color images holistically, rather than as separated color space components. This section illustrates the use of QFTs to process color images through the use of a vector convolution operation. 4.1.1. Classical grayscale image convolution ﬁlters Convolution is a way of describing how systems change an input signal into an output signal. If the system under consideration is a linear time- (or shift-) invariant process, the Fourier domain description of that system is greatly simpliﬁed. Time invariance means that whether we apply an input to the system now or τ seconds from now, the system output will be identical except for a time delay of τ seconds. Hence, the system is time invariant because the output does not depend on the particular time at which the input is applied. A similar description applies to shift-invariant systems. For example, image processing ﬁlters can be described by their point-spread function, which is independent of where, in the input image, the ﬁlter is applied. The classical 2D convolution of two functions, f (x, y) and h(x, y), is given as: (f ∗ h)(x, y) =

∞

∞

−∞

−∞

f (x − x , y − y ) h(x , y ) dx dy

[4.1]

68

Quaternion Fourier Transforms for Signal and Image Processing

where f (x, y) is the input signal and h(x, y) is the system’s point-spread function. The discrete equivalent of this convolution, used in image processing, is given by: N−1 M −1

(f ∗ h)[n, m] =

f [n − i, m − j] h[i, j]

[4.2]

i=0 j=0

where f [n, m] ∈ RN ×M is the input image, a scalar array of grayscale pixels. The ﬁlter convolution mask h[n, m] is typically much smaller in dimensions than the input image. The 2D complex Fourier transform of this equation is given by F {f ∗ h} = F [u, v] H[u, v]

[4.3]

where F [u, v] = F {f } and H[u, v] = F {h}. This is an element-wise product, not a matrix product. This equation, called the Fourier convolution theorem, gives insight into what the ﬁlter mask does to an arbitrary input image. It scales and phase shifts each unique frequency component of the input image, F [u, v], for a speciﬁc frequency index [u, v] by the spectral coefﬁcient of the mask at the same index. In fact, the desired ﬁlter mask can be designed by construction of the mask spectral coefﬁcients and taking the inverse transform of the results. Or, if only one speciﬁc image needs to be ﬁltered, it can be manually shaped in the frequency domain and the result inverse transformed, skipping the construction of the mask altogether. The discrete convolution of [4.2] can be visualized by imagining the mask h[i, j] being shifted across the input image f [n, m] to different offsets; then the superimposed values at this offset are multiplied together, and the products added. The resulting sum at this offset is the value of the output image g[n, m] = (f ∗ h)[n, m] at the pixel location corresponding to the center of the mask. This is shown in Figure 4.1. The classical edge detectors, attributed to Prewitt, Sobel and Kirsch, are examples of 2D convolution ﬁlters deﬁned by [4.3]. The simplest of the three is the Prewitt ﬁlter, which is deﬁned by a real-valued convolution mask over a 3 × 3 neighborhood as: ⎡

1 h=⎣ 0 −1

⎤ 1 1 0 0 ⎦ −1 −1

[4.4]

This mask detects horizontal edges; a transposed version of this mask can be used to detect vertical edges. Both ﬁlters can be combined to yield an edge strength output image. This can be done by applying the second ﬁlter to the output image of the ﬁrst, or

Signal and Image Processing

69

a combined ﬁlter can be created by convolving one mask against the other. Likewise, the classic Sobel ﬁlter has the mask values: ⎡

1 h=⎣ 0 −1

2 0 −2

⎤ 1 0 ⎦ −1

[4.5]

for detecting horizontal edges. Finally, the classic Kirsch ﬁlter has the mask values: ⎡

−3 h = ⎣ −3 −3

5 0 −3

⎤ 5 5 ⎦ −3

[4.6]

and seven other similar masks, each sensitive to vertical or horizontal edges, or edges inclined at 45◦ to the vertical. Each of these edge detection ﬁlters is a high-pass ﬁlter, meaning that low-frequency components of the input image are attenuated and high-frequency components accentuated. This can be seen by computing the complex Fourier transforms of their masks, h[n, m], and examining the resulting spectral coefﬁcients, H[u, v], at speciﬁc frequencies (±u, ±v).

Figure 4.1. Convolution mask operation example

The extension of these ideas to color images could be achieved by repeated application to each color component of the image, as if each component were a grayscale image. This is, in fact, what is commonly done. But this implies that such marginal, component-wise ﬁltering is sufﬁcient under all desired color image ﬁltering operations. It has been shown that this assumption is false, and the rest of this section shows why.

70

Quaternion Fourier Transforms for Signal and Image Processing

4.1.2. Color images as quaternion arrays Color image pixels typically have three components, as such they can be represented in quaternion form using pure quaternions1. For example, a pixel at image coordinates (n, m) in an RGB color space image can be represented as f [n, m] = r[n, m]i + g[n, m]j + b[n, m]k, where r[n, m] is the red component and g[n, m] and b[n, m] are the green and blue components of the pixel, respectively. The same approach may be used with other color space models, for example, a luminance-chrominance color space, such as Y Cb Cr [PAL 98]. Using the 3D pure imaginary quaternion is neither coincidental nor arbitrary. As we saw in Chapter 2, a full quaternion may be interpreted as the ratio of two vectors, the quantity that multiplies one vector to give another. This geometric interpretation is exploited in the design of color ﬁlters. In classical image ﬁltering, we typically have to handle the problem of color space closure. Image pixel values have a ﬁnite range, typically positive and scaled to the interval [0, 1]. If a ﬁltering operation results in a pixel value outside this range, then it must be brought back into the range, ideally without distortion artifacts. One says the range of legitimate pixel values is not closed with respect to the ﬁltering operation. For grayscale images, closure is typically forced by clipping the pixel value at the range boundaries. For color image ﬁlters, the color space is a 3D bounded volume, e.g., the RGB color space is typically mapped to a unit edge-length cube with one corner (the black color) at the origin. The color space closure problem becomes more difﬁcult, especially if we wish to avoid visual distortions in the color image output. To reduce these distortions, the origin is moved to the center of the color cube. This means that all pixel values are now vectors pointing away from mid-gray, the center of the cube. Now closure can be enforced by clipping output pixel values to the value where the pixel vector passes through the cube’s surface. This three-space clipping eliminates hue distortions by holding the orientation of the pixel vector constant during the clipping operation, i.e. only the length of the vector is altered [SAN 04]. 4.1.3. Quaternion convolution There are three possible quaternion convolution deﬁnitions available: left-, rightand bi-convolution deﬁned, respectively, as: N −1 M −1

(hL ◦ f )[n, m] =

hL [i, j] f [n − i, m − j]

[4.7]

i=0 j=0

1 Some color spaces are four-dimensional but because of the tri-stimulus response of the human visual system most are three-dimensional.

Signal and Image Processing

71

N−1 M −1

(f ◦ hR )[n, m] =

f [n − i, m − j] hR [i, j]

[4.8]

i=0 j=0

and N −1 M −1

(hL ≺ f

hR )[n, m] =

hL [i, j] f [n − i, m − j] hR [i, j]

[4.9]

i=0 j=0

where f [n, m] is the input image and hL [n, m] and hR [n, m] are the convolution masks. The bi-convolution operation is shown in Figure 4.2 where the ﬁlter masks sandwich the input image.

Figure 4.2. Bi-convolution mask operation example

When deﬁning quaternion ﬁlter masks for color images, care must be taken to ensure that the convolution equations map pure quaternions into pure quaternions. In general, the product of a vector (pure quaternion) with another vector is a full quaternion; hence, there are constraints on the convolution masks. The simplest example of ensuring this constraint is in the design of color edge detection ﬁlters. This was ﬁrst done in [SAN 98c] with a generalization of the Prewitt ﬁlter shown in [4.4]. To this design was added two additional ﬁlters, generalizations of the Sobel and Kirsch edge detection ﬁlters [SAN 00], all based on vector rotation operations. These ﬁlters are a generalization of the classical grayscale edge detection ﬁlters of Prewitt, Sobel and Kirsch, as described in section 4.1.1. They employ a vector rotation in color space about the gray line of the RGB color space on which all achromatic pixel values lie; hence, they require the use of the bi-convolution equation [4.9].

72

Quaternion Fourier Transforms for Signal and Image Processing

The Prewitt color ﬁlter is deﬁned by two quaternion-valued bi-convolution masks over a 3 × 3 neighborhood as ⎡

1 α⎣ 0 q

1 0 q

⎤⎡ ⎤⎡ ⎤ 1 1 1 1 0 ⎦⎣ f ⎦⎣ 0 0 0 ⎦ q q q q

[4.10]

where q = exp(μπ/2) and α = 1/6 is a scale factor to account for the addition of the six non-zero pixel values, and f denotes the input image. This expression is not a matrix equation but a shorthand notation for the situation shown in Figure 4.2. When μ points along the gray line, the rotation operator Rq [{f }] = q {f } q rotates a pixel value {f } through π radians in the plane normal to the luminance axis or gray line (see [2.3]). Like the classical Prewitt ﬁlter, these masks detect horizontal edges; a transposed version of both masks can be used to detect vertical edges. Where the mask overlies an image region of gradually varying color, the pixel values along the top and bottom rows of the mask will be similar. Since the top row is multiplied by unity (and therefore remains unchanged) and the bottom row is rotated through π about the gray line, the pixel values in the top row and bottom row are positioned opposite each other at an equal distance from the gray line after rotation. When added (using quaternion or vector addition) and scaled, the resultant pixel value lies on or near to the gray line, in a plane normal to the gray line positioned between the planes containing the original pixel values. In contrast, when the mask overlies a horizontal color edge, the result after rotation, vector addition and scaling does not lie near to the gray line and exhibits color. When the mask overlies a horizontal luminance edge where the pixels above and below the edge have similar or identical hues, the chrominance components of the pixels will almost cancel. The cancellation is less perfect where there is a large difference in luminance. Figure 4.3, provided in full color in the color section within this book, shows the result of applying the Prewitt ﬁlter to an image. The ﬁltered image is a full color image, despite appearing as a grayscale image. Where the original color image had gradually varying areas of color, the ﬁltered image has been reduced to gray (achromatic) pixels. Where the original image had a color edge (i.e. a sudden change of color), the ﬁltered image exhibits color and the “edge” pixels have a vector “strength”. The Sobel and Kirsch ﬁlters, ([4.5] and [4.6], respectively), are likewise generalized with the expressions: √ 1 2 1 ⎣ 0 0 8 √ q 2q ⎡

⎤⎡ ⎤⎡ 1 1 0 ⎦⎣ f ⎦⎣ 0 q q

√

2 0 √ 2q

⎤ 1 0 ⎦ q

[4.11]

Signal and Image Processing

73

and ⎡ √ √3q 1 ⎣ 30 √3q 3q

√

5 0 √ 3q

√ ⎤⎡ ⎤⎡ √ √5 √3q ⎦ ⎣ f ⎦ ⎣ 3q 5 √ √ 3q 3q

√

5 0 √ 3q

√ ⎤ √5 ⎦ √5 3q

[4.12]

where q is as in the Prewitt ﬁlter.

Figure 4.3. Filtered output (right) of color image (left) using Prewitt edge detector (see color section within this book for a full color version of this ﬁgure)

Finally, it is worth noting that the color of the edges in the output image is related to the colors involved in the edge transitions of the input image. That is, the colors of regions on either side of the edge determine the color of the ﬁnal output edge value. Also, reversing the transition, e.g., red to blue versus blue to red, will complement the output edge color. 4.1.4. Quaternion image spectrum Using the left-sided QFT of [3.41], with transform axis chosen to be √ μ = (i + j + k)/ 3 (which points along the gray line in color space), we obtain the quaternion image spectrum. This quaternion spectrum is an array of quaternion values, with one real and three imaginary parts. We customarily treat this spectrum as a quaternion image – that is, we regard each spectral coefﬁcient as a quaternion pixel. It is usual for all four parts of any given spectral coefﬁcient to be non-zero. A similar statement applies to the complex spectra of grayscale images, where the coefﬁcients usually have non-zero real and imaginary parts, even when the image pixels are purely real. Complex spectra are interpreted and visualized by creating images based on writing the spectrum in polar form: one image for modulus and the other for phase.

74

Quaternion Fourier Transforms for Signal and Image Processing

In a similar manner, a meaningful visualization of the quaternion spectrum can be based on quaternion polar form (see [1.47]). However, quaternions in polar form have a third quantity besides modulus and phase: i.e. the axis. This √ has no counterpart in the complex case because the imaginary operator I = −1 in the complex √ case is constant. In the quaternion polar form, this imaginary operator μ = −1 may be any pure unit quaternion. The modulus is visualized as a grayscale image, appropriately scaled using, for instance, a logarithmic scaling such as

l(q) =

log (1 + |q|) log (1 + K)

where K is the largest modulus in the image. The phase can be visualized by using a color representation based on the IHS color space. We create an RGB value from an IHS value, where H is the phase, I is the mid-range and S is unity. This means that a phase of zero is represented by a red hue, a phase of π/2 is represented by a green hue and a phase of π is represented by a green-cyan hue. This leaves us with the axis to represent. Because μ is a pure quaternion, it has three components. This suggests a representation using an RGB image. The standard RBG color space is restricted to one octant of 3-space, as each color component must be non-negative. The values of μ, however, may be anywhere on a unit sphere centered at the origin. By scaling this sphere by 12 and placing it at the center of a unit edge-length color cube, we can represent all possible values of μ. To further enhance the image, the μ value is scaled until it extends to the surface of the color cube, so that the resulting image is fully saturated. This is done by scaling μ so that the largest component has unit value. Figure 4.4, provided in full color in the color section within this book, shows the quaternion spectrum of an image, visualized as described in the previous text, with the constant term in the center of the image, as is conventional. Positive frequencies are, therefore, in the upper right quadrant.

Figure 4.4. Quaternion spectrum of an image. Left to right: original image, spectral modulus, phase and axis (see color section within this book for a full color version of this ﬁgure)

Signal and Image Processing

75

Given a quaternion spectrum image, F [u, v], its pixel-by-pixel symplectic decomposition with respect to Cμ is given as F [u, v] = F1 [u, v] + F2 [u, v]μ2 ×M where each part F1 and F2 are in CN . The orthogonal vector μ2 can be placed μ anywhere in the chromatic plane, but in the pictured example, it is parallel to the red line. If one takes the inverse transform of each component: f1 = F −1 {F1 } and f2 = F −1 {F2 }, then it happens that the simplex part f1 contains the luminance information and the perplex part f2 contains the chrominance information. This is not a coincidence, as μ is parallel to luminance and everything perpendicular to μ is in the chrominance plane. The simplex part will be purely imaginary, luminance requiring only one real number. The perplex part will be complex, chrominance requiring two real numbers. Figure 4.5, provided in full color in the color section within this book, shows the previous image decomposed in this manner. Both parts are still full color images, not grayscale images. The original image can be recovered by adding the simplex and perplex parts.

Figure 4.5. Symplectic decomposition of color image into simplex and perplex parts of spectra. Left to right: original image, simplex part, perplex part and direct sum of both parts (see color section within this book for a full color version of this ﬁgure)

The arrangement of coefﬁcients in the quaternion spectrum follows the same layout as that in a complex Fourier transform spectrum; the spectrum is thus divided into quadrants as shown in Figure 4.6. Any given pair of horizontal and vertical frequencies is represented in the spatial frequency domain by two quaternion values, the pair being taken from diagonal quadrants. Exceptions are the Nyquist frequency pair, and the DC term, and any pair of frequencies including a DC or Nyquist term. (The coefﬁcients in the two quadrants correspond to combinations of positive and negative horizontal and vertical spatial frequencies.) If the image is to be a pure quaternion image, which will be the case if color images are represented by the color components in the three components of the vector part of a quaternion pixel, the scalar parts of the two frequency components must cancel out. Certain symmetry conditions on the image (analogous to conjugate symmetries in complex images) will

76

Quaternion Fourier Transforms for Signal and Image Processing

ensure this; in general, an image with no symmetries whatever (e.g. a random image) must have a spectral domain representation that satisﬁes this constraint; otherwise the inverse QFT of the spectral coefﬁcients would not reconstruct the original image (the scalar part would be non-zero).

Figure 4.6. Spectral coefﬁcient layout

Consider the inverse transform that reconstructs an image from a pair of spectral coefﬁcients and the basis functions: M −1 N−1

f [n, m] =

eμ2π ( M

mu

+ nu N )

F [V, U ]

v=0 u=0

and extract two terms from the summation, which make up the component of the image at speciﬁc horizontal and vertical spatial frequencies. Let the arbitrary frequencies, denoted by their index values, be V and U (where 0 < V < M/2 and 0 < U < N − 1). From the symmetry of the coefﬁcients, those at M − V and N − U also represent the same frequencies; so for a given choice of V and U , we have, after some simpliﬁcation: f [n, m] = eμ2π( M

mu

+ nu N )

F [V, U ] + eμ2π(

−mu −nu M + N

) F [M − V, N − U ]

= eμ2π(+α+β) F [V, U ] + eμ2π(−α−β) F [M − V, N − U ] Any coefﬁcient of the Fourier spectrum may be separated into simplex and perplex parts with respect to μ: F [V, U ] = Fs [V, U ] + Fp [V, U ]μ2 , where Fs , Fp ∈ Cμ .

Signal and Image Processing

77

The simplex and perplex parts of the coefﬁcient multiply the exponential basis function of the Fourier transform to produce a cosinusoidal image component with horizontal and vertical frequencies determined by the indices of the coefﬁcient, and with amplitude, phase and orientation in color space determined by the value of the coefﬁcient. f [n, m] = eμ2π(+α+β) {Fs [V, U ] + Fp [V, U ]μ2 } + eμ2π(−α−β) {Fs [M − V, N − U ] + Fp [M − V, N − U ]μ2 } =

eμ2π(+α+β) Fs [V, U ] + eμ2π(−α−β) Fs [M − V, N − U ] + eμ2π(+α+β) Fp [V, U ]μ2 + eμ2π(−α−β) Fp [M − V, N − U ]μ2

Figure 4.7. Spectral coefﬁcient images for the constant term and ﬁrst three harmonics (horizontally and vertically) of color image. The DC term is the center tile on bottom row with positive horizontal frequencies to the right and negative horizontal frequencies on the left. The bottom row corresponds to the DC vertical term with each previous row increasing in vertical frequency – the left and right halves of the last row are mirrors of each other as they contain the same spectral coefﬁcient pairs (see color section within this book for a full color version of this ﬁgure)

Since the ﬁrst two terms in the last equality are isomorphic to the complex numbers (being in Cμ ), we see that the coefﬁcients F [V, U ] and F [M − V, N − U ] behave analogously to the classical complex case. That is, these coefﬁcients scale the magnitude and alter the phase of the Fourier exponentials that combine to form a cosinusoidal function. The same is true for the remaining terms in the last equality, after factoring out the μ2 term from both. This is shown in Figure 4.7, which is provided in full color in the color section within this book. When these two cosinusoidal functions are added, they describe an orbit in the color space centered at

78

Quaternion Fourier Transforms for Signal and Image Processing

mid-gray. A demonstration of the color orbit interpretation of the spectral components can be made possible through the use of a color-cube scatter plot. If each pixel value of an image is plotted in color space, the resulting scatter plot shows the distribution of the image’s color contents in color space. The ﬁrst row of Figure 4.8, provided in full color in the color section within this book, shows an image and its color-cube scatter plot. If the original image is quaternion Fourier transformed, a single harmonic (in one direction) extracted and inverse transformed, the resulting image will show a “rainbow grating” pattern due to harmonic oscillations in luminance and chrominance. The scatter plot of this rainbow grating image should draw an elliptic orbit about the center of the cube. The second through sixth rows of Figure 4.8 show that, indeed, each harmonic describes an elliptical orbit in the color space, which results in a rainbow grating at that harmonic. The rightmost column shows the resulting image built up by progressively adding each new harmonic (horizontal and vertical) to the previous harmonics.

Figure 4.8. Harmonic reconstruction of color image. The legend at the top shows a color-cube scatter plot for the original image. The lower grid shows plots of the color orbits for the ﬁrst ﬁve horizontal and vertical harmonics taken from the QFT spectrum of the image, along with their corresponding images. The ﬁnal, rightmost column shows the image resulting from the sum of the harmonics, starting from the constant (zero frequency) term to the corresponding harmonic (see color section within this book for a full color version of this ﬁgure)

Signal and Image Processing

79

Figure 4.9, provided in full color in the color section within this book, shows results of ﬁltering the image directly in the quaternion spectral domain. The ﬁrst column shows the original image and its spectral modulus. The second column shows the results of low-pass ﬁltering the image by zeroing the 12th and higher spectral harmonics. As expected, the resulting image is blurred consistent with this operation. The third column repeats the process in bandpass fashion, clearing the constant and ﬁrst two harmonics as well as the 26th and higher harmonics. As expected, only moderate-sized regions with uniform luminance or chrominance can be detected in the ﬁnal image. Finally, the last column shows high-pass ﬁltering by zeroing the ﬁrst ﬁve harmonics. Again, as expected, the image contains content only where there is rapid luminance or chrominance variation.

Figure 4.9. Direct ﬁltering in frequency domain of color image. Top row, left to right: original image, low-pass, band-pass and high-pass ﬁltered. The bottom row shows the corresponding spectral magnitude (see color section within this book for a full color version of this ﬁgure)

4.2. Generalized correlation Correlation is concerned with detecting and locating similarity in signals or images. Given a known feature such as a pulse of a particular shape and duration, or, in an image, an object seen from a known viewpoint, correlation can locate the feature, if present, within a longer, noisy, signal; or within a noisy or degraded image. For example, it may be desired to locate a target2 object within a defocused or blurred image; or in the context of color image processing, an image taken under 2 The use of the word “target” here is not meant to imply a military application, although of course, there are military applications for object location within images.

80

Quaternion Fourier Transforms for Signal and Image Processing

lighting conditions such that the colors of the target object are different from a reference image. More precisely, correlation will reveal the most likely location, or multiple locations, of the feature. A second use for correlation is to compute the temporal offset between two similar signals or the spatial offset between two images in order to bring one into registration with the other by applying a shift to one signal or image equal to the offset between them. As with any other fundamental technique, correlation can be used as a step within a more complex algorithm. A closely related concept is phase correlation [KUG 75] in which the normalized cross-power spectrum is computed, from which the inverse Fourier transform yields a correlation peak at an offset representing the best correlation between the two signals or images. Phase correlation works particularly well with images, where sharp edges are frequently present (these have frequency content that can be exploited well by the phase correlation technique). Correlation was generalized to quaternion signals and images in 2003 [MOX 03], and we consider it in detail in section 4.2.2. Phase correlation was successfully generalized to color images using QFTs in 2001 [SAN 01], and we consider it in detail in section 4.2.3. In the next section, we discuss the classical cases of correlation and phase correlation, before introducing some quaternion equivalents. The classical case is normally presented in the literature for real-valued signals, but in fact it can easily be generalized to complex-valued signals, which is a more useful starting point for the quaternion case in this book. The principal difference concerns the relationship between time (or spatial) domain reversal (of a signal or image), and the corresponding effect on the frequency-domain representation. Time reversal of a signal corresponds, in general, to interchange of the positive and negative frequencies in the spectrum of the signal. However, the Fourier spectrum of a real signal has conjugate symmetry, that is the positive and negative frequency coefﬁcients are conjugates of each other, and therefore the effect of interchanging positive and negative frequencies can be obtained more easily (without moving data) by conjugating the Fourier coefﬁcients, and it is the latter operation that is mentioned in most texts. Correlation can be deﬁned (and computed numerically) in both the time and frequency domains, but phase correlation has to be computed in the frequency domain, because it requires the phase of the Fourier coefﬁcients3. There is an equivalence between the two domains. In each domain, there is a formula – in the time domain, this requires an integral (or summation in the discrete case), whereas in 3 There is a normalization involved using the modulus of the products of the Fourier coefﬁcients, and this operation is not linear. It does not therefore have a straightforward equivalent in the time domain.

Signal and Image Processing

81

the frequency domain it is a point-wise product between Fourier coefﬁcients. The difference is computationally signiﬁcant (the Fourier domain approach is much faster for practical signals and images). However, it is necessary to know the frequency domain formula, and this depends on the QFT(s) used. In the quaternion cases presented after the next section, it will be shown that the formula is not a trivial generalization of the classical case. Formulas of this type are not easily discovered, as was noted by J.L. Synge in 1972: It is typical of quaternion formulae that, though they be difﬁcult to ﬁnd, once found they are immediately veriﬁable. [SYN 72] One of the ways that these formulas are veriﬁable is to implement them, of course, and check that the implementation behaves as expected, but a very useful fact to remember is that the formula must reduce to the classical case when the quaternion variables within it are constrained to a subspace that is isomorphic to the complex numbers, for example, when all the quaternions within the formula have the same axis (direction of the vector part). 4.2.1. Classical correlation and phase correlation Cross-correlation4 can be deﬁned directly in the time domain (without Fourier transforms) using the following formula [WEI 14a]: c(t) =

∞ −∞

f (τ )g(t + τ ) dτ

[4.13]

The conjugate in the formula is vital if the signals are complex; otherwise, the product of the imaginary parts of two similar signals could cancel with the products of the real parts. Conjugating one of the signals negates one of the imaginary parts and thus cancels out the −1 arising from the product of the two imaginary roots of −1. Of course, if the signals are real, the conjugate has no effect. Note that the formula is equivalent to convolution (compare [4.2]) with one signal time-reversed and conjugated [WEI 14a]. Note that one of the signals must be shifted through all possible shifts relative to the other signal. In the discrete case (with signals of ﬁnite length), the formula becomes: M −1

c[n] =

f [m]g[n + m]

[4.14]

m=0

4 If the mean, or zero-frequency, component of the signals is subtracted out (i.e. f and g in the formula have zero means), then the result is known as the cross-covariance.

82

Quaternion Fourier Transforms for Signal and Image Processing

assuming the signals are of the same length, M samples. For any reasonable length of signal or size of image, a much faster calculation is possible using Fourier transforms. In order to do this, it is necessary to know the operational formula that expresses the equivalent of [4.14] in the frequency domain. In the classical case, this is: c[n] = F −1 {F {f } F {g}}

[4.15]

Note that the product of the two Fourier spectra is computed element-wise (sometimes known as point-wise), that is, each coefﬁcient in the spectrum is multiplied by the corresponding coefﬁcient in the other spectrum (and one of the two is also conjugated). A subtlety about the use of discrete Fourier transforms is that they assume that the signal is periodically repeated over all time (in image processing this is often referred to as tiling – the image is assumed to repeat over an inﬁnite plane both horizontally and vertically). A consequence is that the formula in the time domain must be modiﬁed if it is to yield exactly the same result as the formula based on discrete Fourier transforms, by implementing cyclic wrap around of the shifted signal. This is easily done by modulo arithmetic on the sample index, so that the value m − n in [4.14] is computed modulo M , and hence always has a result falling in the range between 0 and M − 1. A more common approach with real signals is to use zero-padding. Figure 4.10 illustrates classical cross-correlation with complex signals by showing the cross-correlation of two noisy signals, each containing a Gaussian-modulated complex exponential pulse. The cross-correlation was computed using [4.15]. The correlation result is complex, with a waveform matching the two signals that were correlated, and a position along the time axis corresponding to the displacement in time between the two original pulses. Phase correlation is a simple modiﬁcation of the above ideas, requiring only one extra step. Originally developed for grayscale image processing [KUG 75], the technique works best on signals with edges, rather than smoothly varying signals. (Images tend to have sharp transitions between amplitude levels – at least in a properly focused image, due to the edges of objects.) The technique was shown to work for complex images with pixels representing the chrominance components of a color image [SAN 98b, section 12.2.1]. In the complex case, given two signals or images f and g, the formula required to compute phase correlation is: F −1

F G |F G|

[4.16]

Signal and Image Processing

0.5 0.4 0.3 0.2

ℑ

0.1 0 −0.1 −0.2 −0.3 −0.4 0

−0.6 −0.4

500

−0.2

1000

0

1500

0.2 2000

0.4

ℜ

Time

15 10 5

ℑ

0 −5 −10 −15 0

−20 500

−10 1000

0 1500

10 2000

20

ℜ

Time

Figure 4.10. Complex cross-correlation of two Gaussian-modulated complex exponentials. Upper: original signals. Lower: the complex cross-correlation result. In both plots, the real and imaginary parts are projected onto the lower plane and rear plane of the bounding axes

83

84

Quaternion Fourier Transforms for Signal and Image Processing

where F = F {f } and G = F {g}. The normalized element-wise or point-wise product of F and G (whose inverse Fourier transform yields the phase correlation) is called the cross-power spectrum. Care must be taken with the normalization because a divide-by-zero is possible if any frequency coefﬁcient has zero amplitude (this is unlikely with a real signal, but it can easily happen with simulated signals) 5. Usually, further processing is applied to the result to extract the location of the phase correlation peak or peaks. The simplest method is to take the modulus of the phase correlation result and look for the maximum sample or pixel value by thresholding. Figure 4.11 shows the result of applying phase correlation to a pair of signals similar to those illustrated in Figure 4.10. An important difference is that the amplitude of the complex exponential has been discretized in order to introduce “edges” into the signal. The phase correlation technique is less suited to smoothly varying signals, but well suited to signals with abrupt amplitude transitions. Notice the sharp peak that results from suppression of modulus information in the cross-power spectrum (essentially the result is based on the phase only, hence the name of the technique). Note that the location of the peak produced by phase correlation can sometimes be determined to subsample or subpixel accuracy even though the image is sampled on a discrete grid. (As a simple example, consider what happens when the offset between the two signals is half a sample period: the correlation peak will have equal amplitude at two adjacent samples and therefore the location of the peak can be inferred as midway between the two samples). Cross-correlation and phase correlation, as we have seen from the examples just presented, are deﬁned for both real and complex signals or images. To extend the ideas to vector-valued signals or images (with three or four components per sample), it is obviously possible to compute independent classical correlations on each component (marginal processing). It should be clear without much thought that this is unlikely to work well, although in many cases it will yield a workable result. An obvious difﬁculty is that three or four separate marginal results will be obtained, and these have to be combined in some way, meaning that an additional algorithmic step is needed, such as a majority vote, or some numerical combination of the marginal results. In this book, we argue that a hypercomplex approach is a much more sensible way to proceed, and correlation is no exception. 5 Obviously, there can also be numerical problems with frequency coefﬁcients that have very small values, and it is wise to apply a threshold criterion to suppress such coefﬁcients.

Signal and Image Processing

2 1.5 1

ℑ

0.5 0 −0.5 −1 −1.5 −2 0

−2 500

−1 1000

0 1500

1 2000

2

ℜ

Time

0.7

0.6

Amplitude

0.5

0.4

0.3

0.2

0.1

0

0

200

400

600

800

1000 Time

1200

1400

1600

1800

2000

Figure 4.11. Phase-correlation example. Upper: Gaussian modulated complex exponential pulses with discretized amplitude. Lower: modulus of the phase correlation, showing a sharp peak corresponding to the displacement between the two pulses

85

86

Quaternion Fourier Transforms for Signal and Image Processing

4.2.2. Quaternion correlation To generalize correlation to quaternion-valued signals and images, we use [4.13], and replace the complex conjugate by the quaternion conjugate. That this is the correct approach should be evident from the following: the formula must reduce to the classical complex case when all the quaternion values involved come from the same subﬁeld (i.e. they all have the same axis 6 in their vector parts). This is easily seen to be the case, since if this is so, we may replace the axis value by the complex root of −1 and replace each quaternion by an isomorphic complex number, and the formula then reduces to [4.13] because the quaternion conjugate under these conditions reduces to the complex conjugate. We can then immediately compute cross-correlation using [4.14], replacing all the complex operations by their quaternion counterparts. As already discussed, this is a slow method, particularly for images where the run-time depends on the square of the image size, but it provides a very useful reference case for veriﬁcation of a frequency domain implementation. To work in the frequency domain, we need the quaternion equivalent of [4.15], which is unfortunately not obtained simply by changing all the complex operations to their quaternion counterparts. Furthermore, the frequency domain formula depends on our choice of QFT. Unlike the complex case, where there is only one formulation of the Fourier transform, in the quaternion case, because of non-commutative multiplication, there are many possibilities. Section 3.1.5 explained that a frequency-domain correlation formula may be best derived using a mix of left- and right-sided transforms. This was ﬁrst done for the quaternion case in 2000 [ELL 00b] using 2D one-sided QFTs. What follows therefore is an example of the derivation of a frequency-domain correlation formula. A classical complex derivation is given in [WEI 14b], and we follow a similar process here, but in the reverse order, because the example given in [WEI 14b] would require us to choose the Fourier transforms as the ﬁrst step. Instead, we choose the inverse transform as the ﬁrst step and manipulate the integrals until we can see which forward transforms are needed. At various points in this derivation we could make a different choice, and thereby arrive at a different end result. Practically, once we have one formula that is correct and implemented, there is little to gain by deriving others, since the speed of computation is unlikely to differ, and all the formulas will compute the same result. We start with the quaternion equivalent of [4.13]. This derivation could be done with the discrete form of correlation, as in [4.14], but we choose to use the integral 6 See [1.36] for the deﬁnition of the axis of a quaternion.

Signal and Image Processing

87

form in order to make use of the results given in section 3.1.1. Hence, we have the following deﬁnition of quaternion cross-correlation: ∞

c(t) =

−∞

f (τ )g(t + τ ) dτ

[4.17]

where c, f and g are now quaternion-valued signals. Take the right-sided QFT of both sides, as in [3.5], but for simplicity, we will omit the scale factor, and we will also leave implicit the transform axis (+μ unless stated otherwise): F R {c(t)} = r(ω) =

∞

∞

−∞

−∞

f (τ )g(t + τ ) dτ eμωt dt

Since the exponential is constant with respect to the inner integration, and dτ is real (and therefore commutes with the exponential), we may re-arrange to obtain: r(ω) =

∞

∞

−∞

−∞

f (τ )g(t + τ )eμωt dτ dt

Changing the order of integration and moving f (τ ) outside the new inner integral: r(ω) =

∞ −∞

f (τ )

∞ −∞

g(t + τ )eμωt dt dτ

we recognize the inner integral as a right-sided transform and applying the shift theorem from [3.16], we can replace it with F R {g(t)}eμωτ : r(ω) =

∞ −∞

f (τ )F R {g(t)}eμωτ dτ

Now the problem is one of ordering. We need to place the exponential that we have just inserted adjacent tof (τ ). This cannot be done without a decomposition because the exponential does not commute with F R {g(t)}. However, by decomposing F R {g(t)} into components parallel and perpendicular to μ using [3.21], it becomes possible to change the order, in one case requiring a conjugate. r(ω)

= =

∞ −∞ ∞ −∞

R f (τ ) F R {g(t)} + F⊥ {g(t)} eμωτ dτ

f (τ )F R {g(t)}eμωτ dτ +

∞ −∞

R f (τ )F⊥ {g(t)}eμωτ dτ [4.18]

88

Quaternion Fourier Transforms for Signal and Image Processing

The parallel component commutes with the exponential, whereas the perpendicular component anti-commutes: when the order is changed, the exponential must be conjugated:

r(ω) =

∞ −∞

f (τ )eμωτ dτ F R {g(t)} +

∞ −∞

R f (τ )e−μωτ dτ F⊥ {g(t)}

The remaining integrals can now be recognized as right-sided transforms, in one case with a negated transform axis, hence: R R r(ω) = F R f (t) F R {g(t)} + F−μ f (t) F⊥ {g(t)}

[4.19]

and c(t) = F −R {r(ω)}. Note that we have arrived at a result using only right-sided transforms, all in the same direction (forward), but in one case with a negated axis. Only one inverse transform is required as the ﬁnal step. This means that if there is asymmetry in the scale factors between the forward and inverse transforms, this formula requires no modiﬁcation, as it would if we mixed forward and inverse transforms. It remains to explain the meaning of the quaternion values in the correlation result, since it is a quaternion-valued function. It was shown in 2003 [MOX 03] that the quaternion values represent an averaged rotation that in some sense is a best ﬁt rotation to map the quaternion values in one signal into the values in the other signal for each offset between the signals. This value has a more precise meaning when there is a clear correlation “peak”, the position of which identiﬁes the offset that best matches one signal to the other. 4.2.3. Quaternion phase correlation Phase correlation works well on images because of the presence of sharp edges (in most images). Therefore, we present here an example of color image phase correlation. Figure 4.12, provided in full color in the color section within this book, shows two images with identical backgrounds. Each image contains a car, of identical model, but different colors (one is red and the other is yellow). The two cars are placed at different positions in each image (approximately displaced, but also slightly rotated relative to the camera). Figure 4.13, provided in full color in the color section within this book, shows results from phase correlation based on the same operational formula as was given in the previous section, apart from normalization of each value in the cross-power spectrum by dividing by its

Signal and Image Processing

89

modulus (exactly as in the complex case and subject to the same provisos about very small values noted immediately after [4.16]): p = F −R

r(ω) |r(ω)|

[4.20]

where r(ω) is as given in [4.19]. The maximum amplitude peak in the phase correlation surface is found, and its coordinates used as an offset to shift (cyclically) the second car image. When added to the ﬁrst car image (and normalized by dividing by 2), the result in Figure 4.13 is obtained. Note that the two cars have been superimposed correctly, although a small rotation is apparent (the two cars do not align precisely). This result is achieved despite the complex background of the image (which admittedly is invariant between the two images). As can be seen from the phase correlation surface, correct detection of the offset between the two cars is only just achieved. There are 20 points with a magnitude greater than 75% of the amplitude of the largest peak, showing that the detection of the offset is close to being incorrect. It was shown in 1998 that complex phase correlation working on chrominance data could also achieve a phase correlation result for this pair of images [THO 98], and the result is in agreement with that obtained above.

Figure 4.12. Two color images with identical backgrounds but a foreground object (the car) in different positions with different colors. (Images by and reproduced with the permission of Dr. A.L. Thornton [THO 98]) (see color section within this book for a full color version of this ﬁgure)

Tests with grayscale versions of the two images showed that classical phase correlation was unable to establish the offset between the two cars. More than 6,000 points in the phase correlation surface had magnitudes greater than the point at the correct offset of (4,169) found in the quaternion/color case.

90

Quaternion Fourier Transforms for Signal and Image Processing

Figure 4.13. Phase correlation results using the color images in Figure 4.12. Upper: the phase correlation surface. Lower: the two images superimposed after shifting one of them by the coordinate offset found from the maximum peak in the phase correlation surface, showing the two cars correctly superimposed (see color section within this book for a full color version of this ﬁgure)

In a practical application, a better result would be obtained by removing the background from a “reference” image of the car so that the extraneous detail in the reference image would not correlate with details in the “target” image.

Signal and Image Processing

91

Quaternion phase correlation has been shown to work when correlating a color image with a grayscale image [SAN 01], indicating that there could be applications in image registration (e.g. to register an aerial photograph with a map). 4.3. Instantaneous phase and amplitude of complex signals The notions of instantaneous phase and instantaneous amplitude have been well-known concepts in the signal and image areas for a long time [GAB 46]. A very intuitive parallel, which was made in the early years of time-frequency analysis [VIL 48] by physicists, is to link instantaneous phase to the phase velocity and instantaneous amplitude to the group velocity. Instantaneous phase and amplitude can be considered as the simplest versions of a time-frequency analysis for single frequency signals [FLA 98]. In this section, we show how it is possible to deﬁne an instantaneous phase and amplitude for a complex-valued signal using similar techniques as those used for realvalued signals, but by means of quaternion rather than complex Fourier transforms. The work presented here can be seen as an alternate approach to the bivariate signal analysis developed by Lilly et al. [LIL 10], where the authors used vector and matrix algebra together with complex Fourier transforms to analyze bivariate (or equivalently complex-valued) signals. Here, the proposed approach is to stick to scalar (univariate) complex-valued signals and use the QFTs to analyze them. The presented results are then straightforward generalizations of existing results for real-valued signals, which become a special case of the theory developed in the following. In [LIL 10], the authors analyzed complex-valued signals z(t) as bivariate signals that can be represented as trajectories in the plane, as displayed in Figure 4.14. Here, we will rather consider complex signals as trajectories in 3D space, where the third dimension is time t. This is illustrated in Figure 4.15. This choice for a 3D trajectory representation will become obvious when deﬁning the instantaneous phase and frequency concepts in sections 4.3.4 and 4.3.5 as they will be associated with the 3D behavior of the complex signal z(t). 4.3.1. Important properties of 1D QFT of a complex signal z(t) Recall that here we consider complex signals z(t) = zr (t) + izi (t), with real part zr (t) and imaginary part zi (t), which belong to L1 (R; C) ∩ L2 (R; C). Such complex signals have a right-sided QFT with axis j, denoted as ZjR (ω), given by: ZjR (ω) = FjR {z}(ω) =

+∞ −∞

z(t)e−jωt dt

[4.21]

92

Quaternion Fourier Transforms for Signal and Image Processing

Figure 4.14. Representations of a complex signal z(t) as a trajectory in a 2D plane with constant time step from its initial point (black dot) at t = 0 to its ﬁnal point (black circle) at t = 1

Note that we consider a right-sided QFT here, and that similar results could be obtained by choosing a left-sided QFT. Also, in this section, we take the axis of the transform to be μ = j, a special case of the general deﬁnition of the right-sided QFT given in [3.5]. Now, given the quaternion-valued Fourier transform with respect to axis j, denoted Zj (ω) ∈ L1 (R; H) ∩ L2 (R; H), then its inverse transform is: z(t) =

+∞ −∞

Zj (ω)ejωt dω

[4.22]

where z(t) is in general complex-valued. Note that in some cases to be discussed later on, the signal so constructed z(t) may be quaternion-valued, depending on the values taken by Zj (ω). Now, we present some properties of the QFT that will be useful in the deﬁnition of instantaneous attributes. The key property that allows a quaternion extension to what is known for a real signal and its complex-valued analytic signal comes from the symmetry properties of the QFT.

Signal and Image Processing

93

Figure 4.15. 2D and 3D representations of a complex signal z(t). The 2D trajectory in the plane (black signal) is the same as in Figure 4.14. The 3D trajectory (gray signal) is obtained by unfolding z(t) along a third dimension, which corresponds to t. In the 3D representation, each point has coordinate [ (z(t)), (z(t)), t]. The initial point (black dot) is displayed at time t = 0.2 and the ﬁnal point (black circle) is displayed at time t = 1.2 (the 0.2 shift is used for display convenience)

4.3.1.1. Symmetry Given a real-valued signal x(t), it is well known that its (complex) Fourier transform fulﬁlls Hermitian symmetry, such that X(−ω) = X ∗ (ω), where ∗ denotes the complex conjugate. Now, it is also well known that the (complex) Fourier transform of a complex-valued signal z(t) has no symmetry. This lack of symmetry suggests that an alternate Fourier transform may be used that possesses some symmetry properties. This is where the QFT can be the alternative to a classical Fourier transform. The QFT with axis j of a complex signal z(t) ∈ Ci , denoted as Zj (ω), is i-Hermitian symmetric, i.e.: Zj (ω) = Zj (−ω)

i

[4.23]

94

Quaternion Fourier Transforms for Signal and Image Processing

Figure 4.16 illustrates this symmetry property. (Zj (ω))

i (Zj (ω))

✻

✻

✲ ω

0

✲ ω

0

j (Zj (ω))

k (Zj (ω))

✻

✻

✲ ω

0

0

✲ ω

Figure 4.16. Symmetry of the four components of the QFT Zj (ω) of a complex i

signal z(t). The QFT satisﬁes the i-Hermitian symmetry: Zj (ω) = Zj (−ω) , which means that (Zj (ω)) and i (Zj (ω)) are even while j (Zj (ω)) and k (Zj (ω)) are odd

An additional symmetry property is also of interest when considering complex signals z(t). The QFT of the conjugate of z is given by: FjR {z }(ω) = FjR {z}(ω)

j

[4.24]

Signal and Image Processing

95

Note that good symmetries do not exist for the (classical) complex Fourier transform of a complex signal and that the QFT allows us to deﬁne a Fourier transform with symmetries for complex signals. 4.3.1.2. Orthocomplex modulation A key property of the 1D QFT is its behavior with respect to modulation. It is well known that the classical (complex) Fourier transform of a modulated signal consists of a frequency-shifted version of the unmodulated version of the signal. That is, for a given real or complex signal z(t) with FT Z(ω), the FT of z(t)eiω0 t is Z(ω − ω0 ) for a given constant frequency ω0 . For the QFT, things are different when one considers a complex-valued signal z(t). This can be seen by inserting w(t) = z(t)eiω0 t in place of z(t) in [4.21] which gives: Fj {w}(ω) = Wj (ω) =

+∞ −∞

w(t)e−jωt dt =

+∞ −∞

z(t)e−iω0 t e−jωt dt

[4.25]

Now, recall from [1.38] that: eiα e−jβ = e(iα−jβ)

[4.26]

for any α, β ∈ R, then the property of modulation does not extend trivially to the 1D QFT case. In fact, the QFT of w(t) reads: Wj (ω) =

1 (Zj (ω − ω0 ) + Zj (ω + ω0 ) − iZj (ω − ω0 )j + iZj (ω + ω0 )j) [4.27] 2

which cannot be considered as a modulation in the sense that several copies of the QFT of z(t) are shifted and modiﬁed. In order to recover the same frequency shift behavior as in the well-known complex case, it is necessary to use another axis for the exponential that multiplies the signal. In fact, the well-known phase shift appears when the axis of the exponential is the same as the transform axis. So now, consider the signal y(t) = z(t)e−jω0 , then its QFT (with axis j), denoted as Yj (ω) = Fj {y}(ω), is given by: Yj (ω) =

+∞ −∞

z(t)e−j(ω−ω0 )t dt = Zj (ω − ω0 )

[4.28]

Note also that a special case of orthocomplex modulation is that the QFT with axis j of exp(jω0 ) is simply δ(ω − ω0 ). This orthocomplex modulation property extends to QFTs with any axis μ, provided that the exponential axis is also μ.

96

Quaternion Fourier Transforms for Signal and Image Processing

Note that the orthocomplex modulated signal y(t) = z(t)e−jω0 obtained by orthocomplex modulation of the complex signal z(t), is quaternion-valued, i.e. y(t) ∈ H. Thus, when taking a base-band complex signal and orthocomplex modulating it, we get a quaternion-valued signal. As will also be shown in section 4.3.3, a signal with a right-sided QFT is quaternion-valued in the time domain. This will be exploited to propose the quaternion extension of a complex signal in section 4.3.3. 4.3.1.3. Convolution Here, we are interested in the convolution of two complex-valued signals. The very interesting property that exists for such a convolution reads as follows. Consider two complex signals z(t) and w(t) with QFTs Zj (ω) and Wj (ω), respectively, then the following stands: Fj {z ∗ w} (ω) = Zj (ω) where

Wj (ω)

[4.29]

denotes the bicomplex product of [1.65].

Note that the bicomplex product (deﬁned in 1.4.2) is commutative so that, as one would hope, the order in the convolution product does not affect the result, i.e.: Fj {z ∗ w} (ω) = Fj {w ∗ z} (ω)

[4.30]

for any two complex-valued signals z(t) and w(t). This property is rather important as it allows us to deal with convolution products of complex signals and their QFTs just as we are used to doing with real signals and their complex FTs. This simply means that the presented results are generalizations of existing facts and that the QFT offers a thorough generalization of known results to deal with complex-valued signals. 4.3.2. Hilbert transform and right-sided quaternion spectrum As we are interested in deﬁning instantaneous amplitude and phase for complex signals, a simple way to do it is to mimic the Gabor [GAB 46] and Ville [VIL 48] approach that adds a signal in quadrature to the original signal. In their work, Gabor and Ville considered real-valued signals and added a quadrature signal in the imaginary part of a complex-valued signal associated with the original. This complex signal is named the “analytic signal” associated with a real-valued signal. The quadrature signal was obtained by the Hilbert transform of the original signal. Here, our original signal is complex-valued. So, in order to proceed in the manner of Gabor and Ville, we now need to add a complex signal in quadrature to the ﬁrst

Signal and Image Processing

97

one and to make this addition so that the ﬁnal signal is 4D. This is where quaternionvalued signals enter the scene. Before building this quaternion signal, we now take a look at how to build the quadrature signal. For this, we would like to simply build the Hilbert transform of a given complex signal z(t). It is constructed as follows. First, we are interested in the QFT of the real-valued signal x(t) = 1/πt. It is trivial to notice that, due to the behavior of the QFT with respect to real-valued signals (see Chapter 3), we get: Fj {x}(ω) = −j sgn(ω)

[4.31]

This allows us to deﬁne, in a straightforward way, the Hilbert transform of a complex signal z(t) ∈ C, denoted zH (t), as follows: zH (t) = p.v. z ∗

1 πt

(t)

[4.32]

Now, as 1/πt is a real-valued signal, if we look at the QFT of zH (t), it is the QFT of the convolution product between a real and a complex signal. Thus, the previously introduced property [4.29] for the transformation of a convolutive product applies in a degenerate form where one of the two signals has no imaginary part. In this case, as can be checked directly from [4.29], the bicomplex product boils down to a classical product and reads: ZjH (ω) = −j sgn(ω)

Zj (ω) = −j sgn(ω)Zj (ω)

[4.33]

A proof can be found in [LE 14]. As we have deﬁned the quadrature signal zH (t) of a given complex signal z(t), we need to check which type of signal corresponds to a spectrum with only positive frequencies. Consider a spectrum Zj (ω) obtained from the QFT of a complex signal z(t) = zr (t) + izi (t). Then, just as shown in Figure 4.16, the spectrum Zj (ω) has positive and negative frequencies. In order to build from Zj (ω) a spectrum with only positive frequencies, we simply need to consider (1 + sgn(ω)) Zj (ω) = Wj (ω). It is obvious that Wj (ω) = 0, ∀ω < 0, and by linearity of the inverse QFT, we get: Fj−1 {Wj } (t) = Fj−1 {(1 + sgn(.)) Zj } (t) = Fj−1 {Zj } (t) + Fj−1 {sgn(.)Zj } (t) = Fj−1 {Zj } (t) + jFj−1 {−j sgn(.)Zj } (t)

[4.34]

98

Quaternion Fourier Transforms for Signal and Image Processing

where the last equality is valid due to the fact that we are considering a right-sided QFT. From the last equality in [4.34] and equation [4.33], we can see that the signal w(t) = Fj−1 {Wj } (t) is simply the following combination: w(t) = z(t) + jzH (t)

[4.35]

which is the sum of the original signal z(t) and its quadrature version that is put in two extra dimensions spanned by the basis vectors j and k. As a result, we can see that a one-sided spectrum QFT signal w(t) is indeed a quaternion-valued signal, that is w(t) ∈ H. This is the direct counterpart of what is known for the classical FT, where a one-sided FT corresponds to a complex-valued signal. As an example, we mention here what happens when the signal z(t) is band-limited and then orthocomplex modulated. This type of signal was introduced previously and referred to as an orthocomplex modulated signal. So, consider a signal z(t) with a band-limited QFT, i.e. Zj (ω) = 0 for |ω| > ωm for a given ωm . Then, for a frequency ω0 > ωm , the orthocomplex modulated signal z(t)ejωO t has a right-sided QFT, and, as a result, should be quaternion-valued. This can be directly checked from its literal expression as: z(t)ejω0 t

= (zr (t) + izi (t)) (cos (ω0 t) + j sin (ω0 t)) = zr (t) cos (ω0 t) + izi (t) cos (ω0 t) +jzr (t) sin (ω0 t) + kzi (t)sin (ω0 t)

[4.36]

Illustration of this orthocomplex modulated signal z(t)ejω0 is given in Figure 4.17, where the QFT spectrum Zj (ω − ω0 ) is drawn. We now move to the problem of associating a quaternion-valued signal with any complex signal, making use of the behavior of the QFT just presented. Even (Zj (ω)) i (Zj (ω))

Odd j (Zj (ω)) k (Zj (ω))

Table 4.1. Symmetry properties of the components of the QFT Zj (ω) of a complex signal z(t)

4.3.3. The quaternion signal associated with a complex signal Due to the results presented in the previous section 4.3.2, we can now turn to the deﬁnition of a quaternion-valued signal associated with any complex signal z(t). We have just seen that a quaternion-valued signal corresponds to a one-sided spectrum. It

Signal and Image Processing

99

seems now rather clear how to associate a given complex signal z(t) with its quaternion companion. However, one could imagine several ways to make such an association. The one presented here is similar to the association between a real signal and its associated analytic signal: they share the same spectral content. This is explained in detail in [PIC 97], from which we just recall the main idea and follow the lines. (Zj (ω − ω0 ))

i (Zj (ω

✻

✻

ω0

0

j (Zj (ω

✲ ω

ω0

0

k (Zj (ω

− ω0 ))

✲ ω

− ω0 ))

✻

✻

0

− ω0 ))

ω0

✲ ω

0

ω0

✲ ω

Figure 4.17. Four components of the quaternion-valued Fourier transform Zj (ω − ω0 ) of the orthocomplex modulated signal z(t)ejω0 t and where the complex signal z(t) has a QFT verifying |Zj (ω)| = 0 for |ω| > ω0 . The spectrum of z(t) is displayed with the standard symmetries holding for a QFT (with axis j) of complex signals

100

Quaternion Fourier Transforms for Signal and Image Processing

Given a signal z(t), there is an inﬁnite number of choices of corresponding quaternion-valued signals qz (t). In order to uniquely determine this signal, the most natural choice is to ask that the frequency content of the signal should be the same as that of z(t). Now, remembering the symmetry properties of the QFT depicted in Figure 4.16, it is obvious that only half of the spectrum (positive or negative frequencies) of z(t) is necessary to reconstruct the entire spectrum due to the symmetry properties of the QFT (with axis j) of z(t). As a consequence, we can simply choose to keep the right-sided spectrum that contains all the information about z(t). As we just showed in section 4.3.2, the signal constructed as such will be a quaternion-valued signal (obtained, for example, by inverse QFT of the right-sided spectrum), and this quaternion-valued signal, denoted as qz (t) previously, will share exactly the same frequency properties as the “original” signal z(t). In fact, this way of building a quaternion-valued signal from a complex signal allows us to make a one-to-one correspondence between the two signals (through the content of their QFTs). To summarize, if we are given a complex signal z(t), then we can associate with it a quaternion signal qz (t) that we will call the quaternion extension of z(t) and that is simply constructed as: qz (t) = z(t) + jHj {z} (t)

[4.37]

As can be deduced from what we exposed, the quaternion extension of a complex signal can be seen as the equivalent of the analytic signal associated with a real signal as ﬁrst introduced by Gabor [GAB 46]. The property of analyticity is also shared by qz (t) in a certain way. In fact, the quaternion extension can be called pair-wise analytic in the following sense. Consider qz (t) ∈ H, a quaternion extension signal corresponding to a complex signal z(t), with the following Cj -pair form: qz (t)

= ( (qz (t)) + j = (q1 (t), q2 (t))

j

(qz (t)) ,

i

(qz (t)) + j

k

(qz (t)))

[4.38]

Then, the quaternion extension qz (t) is said to be pair-wise analytic as the two complex-valued signals q1 (t) and q2 (t) are analytic signals associated, respectively, with the real and imaginary parts of z(t). Analyticity is understood in the classical sense here, with the slight difference that the complex signals considered are Cj -valued. This is a minor issue as Cj is isomorphic to the classical C ﬁeld. The demonstration of the pair-wise analyticity of q(t) is trivial. It can be guessed easily by understanding how the QFT acts on a complex signal z(t). The real and imaginary parts of the signal are transformed separately and the H-valued spectrum obtained by a QFT can be seen as a pair of Fourier transforms applied to the real and imaginary parts of z(t). The and j parts of the QFT contain information from the real part

Signal and Image Processing

101

of z(t) and the i and k parts contain information from the imaginary part of z(t). Then, by linearity of the Hilbert transformation, we can see that the process of constructing the quaternion extension consists of processing the information contained in the real and imaginary parts of z(t) in parallel. In other words, the quaternion extension of a complex signal z(t) consists of associating with it an analytic pair of signals with components that are the analytic signals associated with the real and imaginary parts of z(t). These analytic signals are then the components of the Cj -pair form of the quaternion extension qz (t). 4.3.4. Instantaneous amplitude and phase We have just constructed a quaternion signal qz (t) that is deﬁned in a unique way from a complex signal z(t), following the tracks of the construction of the analytic signal associated with a real signal. The natural question to ask now is: Can we obtain an instantaneous amplitude and phase of a complex signal z(t) from its quaternion extension qz (t), and if so, what do the amplitude and phase represent? In the analytic signal case, provided that the spectral content of the signal is “simple”,7 the analytic signal is understood as the simplest time-frequency representation. Now it is no longer the best representation to use if the spectral content of the signal is richer. In such cases, we should refer to more sophisticated time-frequency analyses (e.g. see the Wigner-Ville representation in section 4.3.7). However, the links between the analytic signal and time-frequency representations (such as the Wigner-Ville distribution) are well known [VIL 48]. They are investigated for the quaternion case in sections 4.3.8 and 4.3.9. From what was explained in section 4.3.3, one might think that the most natural way to handle the quaternion extension qz (t) is to consider it as a pair of complex signals. It is tempting to deﬁne a pair of instantaneous amplitudes as well as a pair of instantaneous phases. This was done by Lilly et al. [LIL 10] and extended to signals with more than two components in [LIL 11]. However, here we would like to consider the signal qz (t) as a single entity rather than analyzing it component-wise. This is because we would like to take advantage of the geometric nature of quaternions. To do so, it is possible to use the richness of quaternion representations/forms to express the quaternion extension signal in terms of a modulus and a phase. Clearly, there are two polar forms that could be used: the polar and the CD polar form. We will use the latter here because it includes a complex-valued envelope which is of particular interest in the analysis of a complex signal. 7 Here, “simple” means that the spectral content is made up of one bandpass contribution.

102

Quaternion Fourier Transforms for Signal and Image Processing

Given a complex signal z(t) and its quaternion extension qz (t), there is a canonical pair [ρz (t), φz (t)], with ρz (t) ∈ C and φz (t) ∈ C uniquely associated with z(t) by the polar CD form of qz (t). Recall that the quaternion extension qz (t) possesses the same spectral information as z(t) and that it is made up of the signal z(t) and its orthogonal Hilbert transform with respect to the j axis (see [4.37]). Then, using the polar CD form introduced in 1.4.1.5, the quaternion extension qz (t) is: qz (t) = ρz (t)eφz (t)j

[4.39]

The principle of canonical pairs (modulus and phase) associated with a given signal was originally used by Picinbono [PIC 97] to demonstrate the one-to-one relation between a real signal and its instantaneous amplitude and phase deﬁned through the analytic signal [VIL 48]. There is a one-to-one correspondence between a complex signal z(t) and the canonical pair associated with z(t) through its quaternion extension qw (t) (see details in [LE 14]). Now, looking at the polar CD form of the quaternion extension qz (t) in [4.39], it is noticeable that it looks very much like an orthocomplex modulation. Recall that the quaternion extension has, by construction, no negative frequencies. As a consequence, it can be shown (see [LE 14] for details) that φz (t) is not complex but real-valued8. Thus, the polar CD form of qz (t) should be understood as made up of a “low frequency” part ρz (t) (the instantaneous amplitude/modulus) subject to an orthocomplex modulation with the exponential eφz (t)j where φz (t) is the instantaneous phase of z(t). Once again, note that this is a very thorough extension of the concepts well known for real signals when considering their complex extension (analytic signal). The presented results incorporate the analytic case as a special case, which happens when z(t) is real-valued. Some illustrative examples of the concept of instantaneous amplitude and phase for different types of complex signals z(t) are presented in section 4.3.6. 4.3.5. The instantaneous frequency of a complex signal Of great interest in non-stationary signal processing is the concept of instantaneous frequency, which provides information on the time-varying spectral content of a signal. However, it must be remembered that this can only lead to a proper interpretation when the spectral content of the signal is not too complicated. 8 Note that for quaternion-valued signals with a one-sided QFT, the argument/phase in their polar CD form is always real. However, there are some quaternion-valued signals having a two-sided QFT spectrum, but with no speciﬁc symmetry relations. For such signals, the argument/phase of their polar CD form can be complex-valued.

Signal and Image Processing

103

When such cases arise, the notion of instantaneous frequency of a complex signal z(t), denoted by Ωz (t), is directly linked to the instantaneous phase deﬁned previously, as: Ωz (t) =

dφz (t) dt

[4.40]

Once again, remember that the considered phase is real-valued, leading to a realvalued instantaneous frequency. This deﬁnition is again a generalization of what is known for real-valued signals [VIL 48]. Now, recalling that the signal we are looking at, i.e. z(t), is complex-valued and bearing in mind the 2D planar or 3D trajectory that characterizes z(t), it is obvious that we need more than just an instantaneous frequency to analyze such a signal. Due to the 2D nature of its samples, it is possible to think about z(t) as a trajectory evolving in the 2D plane as mentioned earlier. However, here we prefer to consider z(t) as a 3D curve, as illustrated in Figure 4.15. As we have obtained an instantaneous frequency Ω(t) for z(t) and as it is scalar (real) valued, there should be more information necessary to describe the instantaneous behavior of z(t). In fact, in addition to the value of the instantaneous frequency Ωz (t), we need to know in which plane in 3D space is the signal oscillating at each time t. This notion can be understood as the instantaneous osculating plane of z(t). This osculating plane can be obtained from the instantaneous amplitude introduced in section 4.3.4, as explained in [LE 14]. This can be inferred from the classical case of the analytic signal associated with a real signal x(t). The instantaneous amplitude, also known as the envelope of the signal, is obtained from the modulus of the analytic signal. Now, as shown in Figure 4.18, the envelope ρx (t) together with its ﬁrst-order derivative ρx (t) deﬁne the plane in which the signal x(t) oscillates. In the case where x(t) is real-valued, it is trivial as the plane does not change with time (the plane is deﬁned by the time axis and the signal magnitude axis). Now, in the case where the signal z(t) is complex-valued, and remembering that we consider that z(t) is a 3D trajectory (as displayed in Figure 4.15), the envelope is complex-valued, meaning that it can also be interpreted as a 3D trajectory, just like z(t). As a consequence, the derivative of the envelope, ρz (t) is the set of vectors tangent to ρz (t), just like the velocity vector associated with a trajectory (see [DO 76] for example). Now, the osculating plane is simply deﬁned by ρz (t) × ρz (t) where × is the cross product. In the Frenet-Serret formulation [DO 76], the vector deﬁned by this cross product is the binormal vector. It is orthogonal to the osculating plane. This osculating plane concept for complex signals will be illustrated in section 4.3.6. Now, we have shown that in the case of a complex signal z(t), it is possible to deﬁne the concept of an instantaneous frequency Ωz (t). However, additional

104

Quaternion Fourier Transforms for Signal and Image Processing

information is necessary to describe the dynamical behavior of z(t), and this information is provided by the complex envelope ρz (t). First, just as in the real case, it provides a complex envelope to the signal z(t), i.e. an instantaneous amplitude, and in addition to this, it provides the orientation in 3D space of the plane in which z(t) oscillates. The information about this plane and the instantaneous frequency provide complete information about the local behavior of z(t) in 3D. The introduced concepts are illustrated using two simple examples in the following section.

Figure 4.18. In the classical (real) case when x(t) ∈ R , the instantaneous amplitude ρx (t) (dashed black line) of a modulated real signal x(t) (gray) is given by the envelope of its analytic signal. Assuming that ρx (t) is at least C 1 , its derivative ρx (t) (represented at t = t0 ) together with ρx (t) deﬁnes the osculating plane of the signal x(t)

4.3.6. Examples In order to show the information carried by the instantaneous phase and amplitude of the quaternion extension associated with a complex signal, we present two toy examples that should highlight the concepts of instantaneous phase and envelope that were previously introduced. The presented examples are orthocomplex modulated complex signals with different kinds of behavior for their instantaneous frequency.

Signal and Image Processing

105

4.3.6.1. Example 1: orthocomplex modulated signal with constant frequency In this example, we consider a complex signal z(t), which is displayed in blue in Figure 4.19 (provided in full color in the color section within this book) together with its real and imaginary parts (projections). This signal z(t) with band-limited QFT spectrum is also displayed (real and imaginary parts) in Figure 4.20. Note that it is not analytic in the classical sense [VIL 48], meaning that its real and imaginary parts are not in quadrature at any time instant. In Figure 4.19, a 3D representation of the signal z(t) is given together with its projections in the and planes. Figure 4.20 presents the projections of the signal, i.e. its real and imaginary parts (z(t)) and ((z(t)).

Figure 4.19. Orthocomplex modulated complex signal with constant frequency z(t). 3D representation of the complex signal z(t) (blue solid line) and its complex instantaneous amplitude ρz (t) (red solid line) obtained from the polar CD form of its quaternion extension qz (t). The projections in and planes (black solid lines) correspond to the signals presented in Figure 4.20 (see color section within this book for a full color version of this ﬁgure)

From the original complex signal z(t), it is possible to estimate its complex envelope ρz (t) as explained in section 4.3.4. This is done by computing the QFT of z(t), and then by canceling out the negative frequencies and computing the inverse QFT. The result is a quaternion-valued signal qz (t). The complex envelope ρz (t) of this equivalent quaternion signal qz (t), obtained from its polar CD form section 1.4.1.5) is displayed in Figures 4.19 and 4.20. One can see from these ﬁgures (3D and projections) that ρz (t) is a complex-valued envelope. This is thus a

106

Quaternion Fourier Transforms for Signal and Image Processing

complex-valued signal that corresponds to a base-band version of z(t). Now, as explained previously, as qz (t) is quaternion-valued, its polar CD form is: qz (t) = ρz (t)eΦz (t)j

[4.41]

from which we already mentioned the envelope ρz (t). In this example, the signal z(t) was constructed such that it contains only one frequency. This can be seen from Figures 4.19 and 4.20 as the oscillations of z(t) are not varying with time, showing a single frequency content behavior. This means that its instantaneous frequency is constant. This could be seen as well by plotting the estimated Φz (t). As it is a rather trivial case, the display is omitted here and an exposition of instantaneous frequency is postponed to example 2.

Figure 4.20. Orthocomplex modulated complex signal. Real part (top) and imaginary part (bottom) of the complex signal z(t) (gray solid lines) and its complex instantaneous amplitude ρz (t) (dashed black lines) obtained from the CD polar of its quaternion extension qz (t). The envelopes computed using the “classical” analytic signal from (z(t)) and (z(t)) are displayed in black solid lines

Signal and Image Processing

107

It must be noticed from the complex envelope displayed in Figures 4.19 and 4.20 that it is different from envelopes that could be obtained from the classical analytic signals computed separately on the real and imaginary parts of z(t). This can be seen in Figure 4.20 where the envelopes obtained using the classical analytic signal computed on the real and imaginary parts of z(t) separately are displayed. As both the envelopes computed using the classical analytic signal can only be positive valued (as they are moduli of complex-valued signals), there is no chance that they give an envelope that is a base-band version of the original complex signal z(t). It must be emphasized that this is a major advantage of using the QFT to study complex-valued signals, as the complex envelope ρz (t) has a direct interpretation, while any process computing the real and imaginary parts of z(t) may not carry as much insight as the QFT approach does. Another approach to extend the concept of instantaneous attributes to complex signals has been proposed in [LIL 10]. It is a parametric approach called the modulated elliptical signal (MES), but it does not provide a complex-valued envelope for a given complex signal z(t). A comparison of the MES approach with the instantaneous amplitude and frequency presented here is available in [LE 14]. 4.3.6.2. Example 2: orthocomplex modulated signal with linearly varying frequency In this example, we consider a complex signal with an instantaneous frequency linearly varying with time. The original z(t) signal is displayed in Figure 4.21 in 3D and in Figure 4.22 as real and imaginary parts. The frequency content of z(t) is a single frequency that linearly changes with time. This can be guessed from Figures 4.21 and 4.22 where it can be seen that the frequency increases with time from t = 0 on, followed by a small decrease just before time t = 1. The complex envelope ρz (t) obtained from the polar CD form of the quaternion signal associated with z(t) is displayed in 3D in Figure 4.21. In Figure 4.22, the real and imaginary parts of z(t) are presented, together with the complex envelope. The envelopes computed using the classical analytic signal computed separately on the real and imaginary parts of z(t) are superimposed. As we can see in Figure 4.22, the “classical” envelope cannot take negative values, which means that it cannot lead to the type of result obtained with ρz (t). However, for positive values, the “classical” envelope and ρz (t) take similar values. Finally, to illustrate the instantaneous frequency concept, Figure 4.23 presents the estimated instantaneous frequency obtained by numerical differentiation of the phase Φz (t) given by the polar CD form of the quaternion signal qz (t). It must be noted that in theory, the phase Φz (t) could be complex-valued (see the deﬁnition in [1.60] and [1.61]). However, it can only take real values when the real and imaginary parts of the signal share the same frequency, which is the case here and is for most of the time in real-world applications. For example, it is well known that polarized signals can

108

Quaternion Fourier Transforms for Signal and Image Processing

only be called so if both components of the electric ﬁeld (classically denoted by Ex (t) and Ey (t)) have frequencies in common [BRO 98]. As can be seen in Figure 4.23, the estimated frequency contains some errors. These errors occur at time instants where the signal z(t) vanishes, leading to aberrant values in the phase estimation. However, the frequency content of the signal is recovered as can be seen from Figure 4.23 where the increasing frequency behavior together with the decrease at the end of the time window can be observed.

Figure 4.21. Orthocomplex modulated complex signal with linearly varying frequency z(t). 3D representation of the complex signal z(t) (blue solid line) and its complex instantaneous amplitude ρz (t) (red solid line) obtained from the polar CD form of its quaternion extension qz (t). The projections in and planes (black solid lines) correspond to the signals presented in Figure 4.22 (see color section within this book for a full color version of this ﬁgure)

This demonstrates the ability of the polar CD form of the quaternion extension of a complex signal z(t) to access its instantaneous frequency. This is a direct extension of what is known for real signals and the ability of the polar form of the analytic signal to provide the instantaneous phase and frequency. Together, the instantaneous amplitude (envelope), the instantaneous phase and the instantaneous frequency obtained from the polar CD form of qz (t), the quaternion extension of the complex signal z(t), are the faithful generalizations of the well-known instantaneous attributes obtained from the analytic signal of a real-valued signal. This shows the position of the QFT with respect to complex-valued signals: it is the counterpart of the (classical) FT for real signals.

Signal and Image Processing

109

This is due to the inherent symmetries of the QFT and its intrinsic four dimensions. The QFT should thus be considered on a systematic basis for the study of complex-valued signals as it has the required symmetries to extract any spectral information from bivariate signals.

Figure 4.22. Orthocomplex modulated complex signal. Real part (top) and imaginary part (bottom) of the complex signal z(t) (gray solid lines) and its complex instantaneous amplitude ρz (t) (dashed black lines) obtained from the CD polar form of its quaternion extension qz (t). The envelopes computed using the “classical” analytic signal from (z(t)) and (z(t)) are displayed as black solid lines

4.3.7. The quaternion Wigner-Ville distribution of a complex signal We end this section with some results that provide the basic foundations for a generalization of time-frequency representations inspired by the QFT9. 9 A simple time-frequency analysis of non-stationary complex signals can also be carried out using a short-time QFT that can be easily computed using sliding time-windowed QFTs.

110

Quaternion Fourier Transforms for Signal and Image Processing

Time-frequency representations are well known in signal processing for the study of non-stationary signals [FLA 98]. Here, we show how the previously introduced concepts open the door to a new way to study non-stationary complex signals. Some parts of the presented material can be found in [ELL 12].

Figure 4.23. Orthocomplex modulated complex signal with linearly varying frequency. Estimated instantaneous frequency obtained from the numerical derivative of the instantaneous phase computed as the phase of the polar CD form of qz (t), the quaternion signal associated with the complex signal z(t). The estimation errors are due to the signal vanishing at certain times

This section shows the link between the just deﬁned quaternion extension and instantaneous frequency of a complex signal and its quaternionic Wigner-Ville distribution (QWVD). First, we introduce the deﬁnition of the quaternion Wigner distribution (QWD) of a complex signal and then the QWVD of its quaternion extension. Then, we introduce some of the remarkable properties of the QWVD. Finally, we show how the time marginal and mean frequency formulas [FLA 98] extend to the quaternionic case. 4.3.7.1. Deﬁnitions Just as in the case of real signals, we will make the distinction between Wigner and the Wigner-Ville distribution. Recall that classically the Wigner distribution (WD) is built for real-valued signals, while the Wigner-Ville is built using the complex analytic signal associated with a real signal. The latter version was originally proposed by Ville [VIL 48]. Motivations for using Wigner-Ville distributions are given, for example, in [BOA 88].

Signal and Image Processing

111

First, we start with the Wigner distribution deﬁnitions. Given two complex-valued signals z(t) and w(t), their quaternionic (cross-) Wigner distribution with respect to the j-axis, denoted by Wz,w (t, ω), is: +∞

Wz,w (t, ω) =

−∞

z(t + τ /2)w (t − τ /2)e−jωτ dτ

[4.42]

while the quaternionic (auto-) Wigner distribution of z(t), denoted by Wz (t, ω), is given by: Wz (t, ω) =

+∞ −∞

z(t + τ /2)z (t − τ /2)e−jωτ dτ

[4.43]

Now, given two complex-valued signals z(t) and w(t) and their respective quaternion extensions qz (t) and qw (t), their quaternionic (cross-) Wigner-Ville distribution with respect to the j-axis, denoted by W Vz,w (t, ω), is: W Vz,w (t, ω) =

+∞ −∞

qz (t + τ /2)qw (t − τ /2)e−jωτ dτ

[4.44]

while the quaternionic (auto-) Wigner-Ville distribution of z(t), denoted by W Vz (t, ω), is given by: W Vz (t, ω) =

+∞ −∞

qz (t + τ /2)qz (t − τ /2)e−jωτ dτ

[4.45]

Note that in the Wigner-Ville deﬁnitions, the transform consists of a QFT with axis j of a quaternion-valued function with respect to the τ variable. Just as would happen with the complex FT of a complex signal, and because the quaternion extensions qz (t) and qw (t) have one-sided QFT spectra, the quaternionic auto- and cross-Wigner-Ville distributions are one-sided (i.e. equal to zero for negative frequencies). In the sequel, we will consider Wigner and Wigner-Ville distributions, depending on the purpose we wish to illustrate. Mainly, we will stick to the Wigner case when dealing with standard properties, and then move to the Wigner-Ville case when considering the link with the instantaneous amplitude, phase and frequency of a complex signal.

112

Quaternion Fourier Transforms for Signal and Image Processing

4.3.7.2. Basic properties Here, we list several properties of the QWD. Most of the proofs can be obtained by direct calculation, and can be found in [ELL 00a]. First, the QWD of the conjugate signal z (t) is simply: Wz (t, ω) = Wz (t, ω)

[4.46]

and the (cross-) QWD fulﬁlls the following symmetry relation when reversing the order of z(t) and w(t): Ww,z (t, ω) = Wz,w (t, ω)

k

[4.47]

The behavior of the QWD with respect to time shifting is as follows: WTχ z,Tχ w (t, ω) = Wz,w (t − χ, ω)

[4.48]

where Tχ z(t) = z(t − χ) is a translation operator. 4.3.7.3. Inversion It is possible to extend the well-known inversion formulas for the WD [FLA 98] to the quaternion version of the WD. Given two complex signals z(t) and w(t), the following equality holds: z(t)w (0) =

+∞ −∞

Wz,w

t , ω ejωt dω 2

[4.49]

provided that w (0) = 0. Now, in the case of the (auto-) QWD of a complex signal z(t), this becomes: z(t)z (0) =

+∞ −∞

Wz (t, ω) ejωt dω

[4.50]

This means that z(t) can be recovered up to a phase factor z (0) by computing the inverse QFT of its QWD. Now, we previously introduced the orthocomplex modulation which shares some behavior with the classical modulation when considering complex FTs. As was highlighted in section 4.3.1.2, the orthocomplex modulation is a central concept in the study of band-limited complex signals with right-sided quaternionic spectrum. It

Signal and Image Processing

113

is thus interesting to know how this peculiar modulation behaves through the QWD, which is given, for two complex signals z(t) and w(t), by: WOω0 z,Oω0 w (t, ω) = Wz,w (t, ω − ω0 )

[4.51]

where [Oω0 z] (t) = z(t)ejω0 t and [Oω0 w] (t) = w(t)ejω0 t deﬁnes the orthocomplex modulation operator. In addition to the presented properties of the QWD, we now give the two characteristics that make the connection between the instantaneous phase and amplitude introduced in section 4.3.4 and some characteristics of the QWVD, namely the time marginal and the mean frequency formula. 4.3.8. Time marginal We investigate here what information is available at each time instant t in the QWD of a complex signal z(t) by integrating Wz (t, ω) over the frequency range. The result is summarized by the following expression: +∞ −∞

Wz (t, ω) dω = z(t)z (t) = |z(t)|

2

[4.52]

This is simply the squared instantaneous magnitude of the signal z(t), i.e. its instantaneous energy. Doing the same with the QWVD would lead to the squared modulus of the quaternion extension of z(t), which carries the same information in terms of energy. It would read as follows: +∞ −∞

W Vz (t, ω) dω = q(t)qz (t) = |qz (t)|

2

[4.53]

In the case of a real signal, it is well known that proceeding in the same way with the analytic signal associated with the original real signal leads to a similar result. Thus, the QWD (and QWVD) simply generalize this existing result to the case of an arbitrary complex signal. 4.3.9. The mean frequency formula A more interesting result, which is known for the classical Wigner-Ville distribution, is that its mean value is the instantaneous frequency [FLA 98]. It states that the normalized mean value of the Wigner-Ville distribution is in fact the instantaneous frequency. We show here how this extends to the QWVD. The

114

Quaternion Fourier Transforms for Signal and Image Processing

generalization leads to a more complicated expression for the mean value of the quaternion-valued QWVD; however, it is a generalization of the classical case and is given in [4.66]. As an illustration of calculations with a QFT, we give the lines to follow to obtain this property of the QWVD. First, recall that if a complex signal z(t) has a QFT denoted by Zj (ω), then the following property holds: Fj

d z (ω) = Zj (ω)jω dt

[4.54]

This extends to the QFT of the quaternion extension qz (t) directly by: Fj

d qz (ω) = Zj+ (ω)jω dt

[4.55]

where we now have that Zj+ (ω) = 0 for ω < 0. We also have the classical initial value theorem that holds for the QFT of qz (t), i.e.: +∞ −∞

Zj+ (ω) dω = qz (t)|t=0 = qz (0)

[4.56]

Now, from the deﬁnition of the QWVD, we can write that: qz t +

τ τ qz t − 2 2

Fj

−→

W Vz (t, ω)

[4.57]

Now, remembering the above-mentioned property of the QFT with respect to derivation and that in the QWVD deﬁnition, integration is with respect to τ , we get: d τ τ qz (t + )qz (t − ) dτ 2 2

Fj

−→

W Vz (t, ω)jω

[4.58]

Now, expanding the derivative, denoting qz (t) = dqz (t)/dt and taking its value for τ = 0, we get: +∞

W Vz (t, τ )jω dω =

−∞

1 τ τ τ τ qz t + qz t − − qz t + qz t − 2 2 2 2 2

τ =0

[4.59]

Signal and Image Processing

115

Note that we left the integration over the range ]−∞, +∞[ even though we know that the QWVD is right-sided. This does not affect the calculation anyway. Now, rightmultiplying both sides of the equality by j leads to the following equality: +∞ −∞

ωW Vz (t, ω) dω = −

1 q (t)qz (t) − qz (t)qz (t) j 2 z

Dividing both sides of the equality by result in [4.53], we get: +∞ ωW Vz (t, ω) dω −∞ +∞ W Vz (t, ω) dω −∞

=

1 2

= −

+∞ −∞

[4.60]

W Vz (t, ω) dω and remembering the

qz (t)qz (t) − qz (t)qz (t) j qz (t)qz (t) 1 q (t)qz−1 (t) − q −1 z (t)qz (t) j 2 z

[4.61]

= − V qz (t)qz−1 (t) j It sufﬁces now to remember that the quaternion extension qz (t) can be expressed with its polar CD form (as explained in [4.39]), so that qz (t) = ρz (t)eφz (t)j . This means that its time derivative is given by: qz (t) = ρz (t)eφz (t)j + ρz (t)φz (t)jeφz (t)j = {ρz (t) + ρz (t)φz (t)j} eφz (t)j

[4.62]

Now, it follows that: +∞ ω W Vz (t, ω) dω −∞ +∞ W Vz (t, ω) dω −∞

= −V

qz (t) qz (t)

j

ρ (t) + ρz (t)φz (t)j = −V z ρz (t)

[4.63] j

It is important to remember that ρz (t) is complex-valued and that φz (t) can also be complex-valued (even if we already mentioned that it is most likely to be real-valued in real-world applications). Writing ρz (t) = |ρz (t)| exp ϕρz (t) i , we ﬁrst evaluate the quantity qz (t)qz−1 (t) as: qz (t)qz−1 (t) =

|ρz (t)| + iϕρz (t) + exp ϕρz (t) i φz (t)j |ρz (t)|

[4.64]

which is a full quaternion term as |ρz (t)| / |ρz (t)| ∈ R, and iϕρz (t) is a pure i part and exp ϕρz (t) i j is made of a j and a k part. Now, the ﬁnal step consists of

116

Quaternion Fourier Transforms for Signal and Image Processing

taking the vector part of this term and multiplying it on the right by j, which leads to the following: −V

ρz (t) + ρz (t)φz (t)j ρz (t)

j = exp ϕρz (t) i φz (t) − kϕρz (t)

[4.65]

which ﬁnally gives the following expression for the “mean frequency formula”: +∞ ωW Vz (t, ω) dω −∞ +∞ W Vz (t, ω) dω −∞

= exp ϕρz (t) i φz (t) − kϕρz (t)

[4.66]

As said previously, this expression is more complicated than its “classical” version that involves the standard Wigner-Ville distribution of a real signal. Here, however, the instantaneous frequency of the complex signal z(t) (deﬁned as Ωz (t) = φz (t) in section 4.3.5) can be recovered from the normalized mean of the QWVD. First, note that in [4.66], the right-hand side of the equation has two terms that do not mix together in the 4D quaternion space. The term kϕρz (t) is a pure k term that provides information on the complex envelope of z(t). In fact, it is the instantaneous frequency of the complex envelope/amplitude ρz (t) of z(t). Next, the term eϕρz (t) i φz (t) is made of the instantaneous frequency of z(t), multiplied by a pure phase factor containing again information on the complex envelope. Thus, by simply taking the modulus of this term, we recover the instantaneous frequency φ (t) of the complex signal z(t). This is true in the cases where this instantaneous amplitude is real, which is, as already mentioned, the case in most real-world applications. The mean frequency formula thus provides information on the instantaneous frequency of the processed signal, just as in the classical case, and with some extra terms due to the additional information available in a complex signal (with respect to information available in a real-valued signal). As a conclusion, we can see that the QWVD generalizes some well-known results from “classical” WVD, and may be a very powerful tool to analyze complex signals that arise in many applications. There is still some work to be done to completely discover the rest of the properties of the QWVD and of time-frequency representations based on quaternionic Fourier analysis.

Bibliography

[ALF 07] A LFSMANN D., G ÖCKLER H., S ANGWINE S.J. et al., “Hypercomplex algebras in digital signal processing: beneﬁts and drawbacks”, Proceedings of 15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, pp. 1322–1326, 3–7 September 2007. [ALT 86] A LTMANN S.L., Rotations, Quaternions and Double Groups, Dover, 1986. [BAE 02] BAEZ J.C., “The octonions”, Bulletin of the American Mathematical Society, vol. 39, no. 2, pp. 145–205, 2002. [BOA 88] B OASHASH B., “Note on the use of the Wigner distribution for time-frequency signal analysis”, IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 36, no. 9, pp. 1518–1521, 1988. [BRA 74] B RANETS V.N. S HMYGLEVSKIY I.P., Application of Quaternions to Rigid Body Rotation Problems, NASA Technical Translation no. NASA TT F-15414, National Aeronautics and Space Administration, Washington, DC, July 1974, Translation of Primeneniye Kvaternionov v Zadachakh Orientatsii Tverdogo Tela, Nauka Press, Moscow, 1973. Available at http://archive.org/details/nasatechdoc19740020927. [BRA 86] B RACEWELL R.N., The Hartley Transform, Oxford University Press, New York, 1986. [BRA 00] B RACEWELL R.N., The Fourier Transform and its Applications, 3rd edition, McGraw–Hill, Boston, 2000. [BRA 13] B RACKX F., H ITZER E., S ANGWINE S.J., “History of Quaternion and Clifford Fourier transforms and wavelets”, in H ITZER E., S ANGWINE S.J. (eds), Quaternion and Clifford Fourier Transforms and Wavelets, Birkhäuser/Springer, pp. xi–xxvii, 2013. [BRI 13] B RISTOW P.A., KORMANYOS C., H OLIN H., et al., “Boost C++ Libraries, Math Toolkit”, Software Library, 2006–2010, 2012, 2013. Available at http://www.boost.org/. [BRO 98] B ROSSEAU C., Fundamentals of Polarized Light: A Statistical Optics Approach, Wiley-Blackwell, 1998.

118

Quaternion Fourier Transforms for Signal and Image Processing

[BÜL 99] B ÜLOW T., Hypercomplex spectral signal representations for the processing and analysis of images, PhD Thesis, University of Kiel, Germany, August 1999. [BÜL 01] B ÜLOW T., S OMMER G., “Hypercomplex signals a novel extension of the analytic signal to the multidimensional case”, IEEE Transactions on Signal Processing, vol. 49, no. 11, pp. 2844–2852, November 2001. [CHE 00] C HEN B., K AUFMAN A., “3D volume rotation using shear transformation”, Graphical Models, vol. 62, pp. 308–322, 2000. [CON 03] C ONWAY J.H., S MITH D.A., On Quaternions and Octonions, A.K. Peters, Ltd., Natick, MA, 2003. [COX 46] C OXETER H.S.M., “Quaternions and reﬂections”, The American Mathematical Monthly, vol. 53, no. 3, pp. 136–146, 1946. [COX 74] C OXETER H. S.M., Regular Polytopes, Dover Publications Inc., 1974. [DEL 88] D ELSUC M.A., “Spectral representation of 2D NMR spectra by hypercomplex numbers”, Journal of Magnetic Resonance, vol. 77, no. 1, pp. 119–124, March 1988. [DIC 14] D ICKSON L., Linear Algebras, Cambridge University Press, London, 1914. [DO 76] D O C ARMO M.P., Differential Geometry of Curves and Surfaces, Pearson, 1976. [ELL 00a] E LL T.A., S ANGWINE S.J., “Decomposition of 2D hypercomplex Fourier transforms into pairs of complex Fourier transforms”, in G ABBOUJ M., K UOSMANEN P., (eds.), Proceedings of EUSIPCO 2000, 10th European Signal Processing Conference, Tampere, Finland, vol. 2, pp. 1061–1064, 5–8 September 2000. [ELL 00b] E LL T.A., S ANGWINE S.J., “Hypercomplex Wiener-Khintchine theorem with application to color image correlation”, IEEE International Conference on Image Processing (ICIP 2000), Vancouver, Canada, vol. 2, pp. 792–795, 11–14 September 2000. [ELL 07a] E LL T.A., “Hypercomplex color afﬁne ﬁlters”, IEEE International Conference on Image Processing (ICIP 2007), San Antonio, TX, vol. 5, pp. 249–252, 16–19 September 2007. [ELL 07b] E LL T.A., On systems of linear quaternion functions, arXiv:math/0702084v1, February 2007. [ELL 07c] E LL T.A., S ANGWINE S.J., “Hypercomplex Fourier transforms of color images”, IEEE Transactions on Signal Processing, vol. 16, no. 1, pp. 22–35, January 2007. [ELL 07d] E LL T.A., S ANGWINE S.J., “Quaternion involutions and anti-involutions”, Computers and Mathematics with Applications, vol. 53, no. 1, pp. 137–143, 2007. [ELL 12] E LL T.A., L E B IHAN N., “Quaternion Wigner-Ville distribution”, Applied Geometric Algebras in Computer Science and Engineering (AGACSE), La Rochelle, France, 2012. [ERN 87] E RNST R.R., B ODENHAUSEN G., W OKAUN A., Principles of Nuclear Magnetic Resonance in One and Two Dimensions, Oxford University Press, Oxford, 1987. [FLA 98] F LANDRIN P., Time-Frequency/Time-Scale Analysis, Academic Press, 1998.

Bibliography

119

[FRI 05] F RIGO M., J OHNSON S.G., “The design and implementation of FFTW3”, Proceedings of the IEEE, vol. 93, no. 2, pp. 216–231, 2005. Special issue on “Program generation, optimization, and platform adaptation”. [GAB 46] G ABOR D., “Theory of communication”, Journal of the IEE, vol. 93, no. 26, pp. 429–457, 1946. [GIL 08] G ILMORE R., Lie Groups, Physics and Geometry, Cambridge University Press, 2008. [GOL 96] G OLUB G.H., VAN L OAN C.F., Matrix Computations, The Johns Hopkins University Press, 3rd edition, Baltimore and London, 1996. [HAN 06] H ANSON A.J., Visualizing Quaternions, Morgan Kaufmann, 2006. [HIT 07] H ITZER E., “Quaternion Fourier transform on quaternion ﬁelds and generalizations”, Advances in Applied Clifford Algebras, vol. 17, no. 3, pp. 497–517, 2007. [HIT 13] H ITZER E., S ANGWINE S.J., Eds., Quaternion and Clifford Fourier Transforms and Wavelets, Birkhäuser/Springer, Basel, Switzerland, 2013. [ICK 70] I CKES B.P., “A new method for performing digital control system attitude computations using quaternions”, AIAA Journal, vol. 8, no. 1, pp. 13–17, January 1970. [JOL 05] J OLY C.J., A Manual of Quaternions, Macmillan, 1905. [KAN 89] K ANTOR I., S OLODOVNIKOV A., Hypercomplex Numbers: Introduction to Algebras, Springer-Verlag, 1989.

An Elementary

[KUG 75] K UGLIN C.D., H INES D.C., “The phase correlation image alignment method”, Proceedings of IEEE Conference on Cybernetics and Society, San Francisco, CA, IEEE, pp. 163–165, 23–25 September 1975. [KUI 02] K UIPERS J., Quaternions and Rotation Sequences: A Primer with Applications to Orbits, Aerospace and Virtual Reality, Princeton University Press, 2002. [LE 14] L E B IHAN N., S ANGWINE S.J., E LL T.A., “Instantaneous frequency and amplitude of orthocomplex modulated signals based on quaternion Fourier transform”, Signal Processing, vol. 94, pp. 308–318, January 2014. [LIL 10] L ILLY J.M., O LHEDE S.C., “Bivariate instantaneous frequency and bandwidth”, IEEE Transactions on Signal Processing, vol. 58, no. 2, pp. 591–603, 2010. [LIL 11] L ILLY J.M., “Modulated oscillations in three dimensions”, IEEE Transactions on Signal Processing, vol. 59, no. 12, pp. 5930–5943, 2011. [LOU 01] L OUNESTO P., Clifford Algebras and Spinors, Cambridge University Press, Cambridge, 2001. [MAC 06] M ACFARLANE A., Vector Analysis and Quaternions, John Wiley & Sons, New York, 1906. [MCI 98] M C I LHAGGA W., “Colour vision”, in S ANGWINE S.J., H ORNE R.E.N. (eds), The Colour Image Processing Handbook, Chapter 2, Chapman and Hall, pp. 7–25, 1998. [MES 83] M ESERVE B.E., Fundamental Concepts of Geometry, Corrected reprint of 1959 edition, Dover, New York, 1983.

120

Quaternion Fourier Transforms for Signal and Image Processing

[MOX 03] M OXEY C.E., S ANGWINE S.J., E LL T.A., “Hypercomplex correlation techniques for vector images”, IEEE Transactions on Signal Processing, vol. 51, no. 7, pp. 1941–1953, July 2003. [PAL 98] PALUS H., “Representations of colour images in different colour spaces”, in S ANGWINE S.J., H ORNE R.E.N. (eds), The Colour Image Processing Handbook, Chapter 4, Chapman and Hall, pp. 67–90, 1998. [PIC 97] P ICINBONO B., “On instantaneous amplitude and phase of signals”, IEEE Transactions on Signal Processing, vol. 45, no. 3, pp. 552–560, 1997. [POR 95] P ORTEOUS I.R., Clifford Algebras and the Classical Groups, Cambridge University Press, 1995. [PRI 91] P RICE G.B., An Introduction to Multicomplex Spaces and Functions, Marcel Dekker Inc., New York, 1991. [ROB 68] ROBINSON P.D., Fourier and Laplace Transforms, Routledge and Kegan Paul, London, 1968. [ROC 04] ROCHON D., T REMBLAY S., “Bicomplex quantum mechanics: I. the generalized Schrödinger equation”, Advances in Applied Clifford Algebras, vol. 14, no. 2, pp. 231–248, 2004. [ROW 14] ROWLAND T., W EISSTEIN E.W., “Matrix exponential”, MathWorld. Available at Wolfram Web Resource, http://mathworld.wolfram.com/ MatrixExponential.html, January 2014. [SAN 96] S ANGWINE S.J., “Fourier transforms of colour images using quaternion, or hypercomplex, numbers”, Electronics Letters, vol. 32, no. 21, pp. 1979–1980, 10 October 1996. [SAN 98a] S ANGWINE S.J., H ORNE R.E.N., (eds.), The Colour Image Processing Handbook, Chapman and Hall, London, 1998. [SAN 98b] S ANGWINE S.J., T HORNTON A.L., “Frequency domain methods”, in S ANGWINE S.J., H ORNE R.E.N. (eds), The Colour Image Processing Handbook, Chapter 12, Chapman and Hall, pp. 228–241, 1998. [SAN 98c] S ANGWINE S.J., “Colour image edge detector based on quaternion convolution”, Electronics Letters, vol. 34, no. 10, pp. 969–971, 14 May 1998. [SAN 00] S ANGWINE S.J., E LL T.A., “Colour image ﬁlters based on hypercomplex convolution”, IEE Proceedings – Vision, Image and Signal Processing, vol. 147, no. 2, pp. 89–93, April 2000. [SAN 01] S ANGWINE S.J., E LL T.A., M OXEY C.E., “Vector phase correlation”, Electronics Letters, vol. 37, no. 25, pp. 1513–5, 6 December 2001. [SAN 04] S ANGWINE S.J., E LL T.A., “Gray-centered RGB color space”, Second European Conference on Color in Graphics, Imaging and Vision (CGIV 2004), Technology Center AGIT, Aachen, Germany, pp. 183–186, 5–8 April 2004. [SAN 10] S ANGWINE S.J., L E B IHAN N., “Quaternion polar representation with a complex modulus and complex argument inspired by the Cayley-Dickson form”, Advances in Applied Clifford Algebras, vol. 20, no. 1, pp. 111–120, March 2010.

Bibliography

121

[SAN 12] S ANGWINE S.J., E LL T.A., “Complex and hypercomplex discrete Fourier transforms based on matrix exponential form of Euler’s formula”, Applied Mathematics and Computation, vol. 219, no. 2, pp. 644–655, October 2012. [SAN 13a] S ANGWINE S.J., “Perspectives on color image processing by linear vector methods using projective geometric transformations”, in H AWKES P.W. (ed.), Advances in Imaging and Electron Physics, Academic Press: Elsevier Inc., vol. 175, pp. 283–307, 2013. [SAN 13b] S ANGWINE S.J., L E B IHAN N., “Quaternion Toolbox for Matlab®, Version 2 with support for Octonions”, Software Library , 2013. Available at http://qtfm.sourceforge.net/. [SHO 85] S HOEMAKE K., “Animating rotation with quaternion curves”, SIGGRAPH, vol. 19– 3, pp. 245–254, 1985. [SNE 61] S NEDDON I.N., Fourier Series, Routledge and Kegan Paul, London, 1961. [SOM 01] S OMMER G., (ed.), Geometric Computing with Clifford Algebras: Theoretical Foundations and Applications in Computer Vision and Robotics, Springer-Verlag, London, UK, 2001. [STI 10] S TILLWELL J., Mathematics and Its History, Undergraduate Texts in Mathematics, 3rd edition, Springer Science+Business Media, LLC, 2010. [SYN 72] S YNGE J.L., Quaternions, Lorentz Transformations, and the Conway-DiracEddington Matrices, Communications of the Dublin Institute for Advanced Studies, Series A no. 21, Dublin Institute for Advanced Studies, Dublin, 1972. [THO 98] T HORNTON A.L., Colour object recognition using a complex colour representation and the frequency domain, PhD Thesis, Department of Engineering, The University of Reading, 1998. [TRU 53] T RUMPLER R.J., W EAVER H.F., Statistical Astronomy, University of California Press, Berkeley and Los Angeles, CA, 1953. [VIL 48] V ILLE J., “Théorie et Applications de la Notion de Signal Analytique”, Cables et Transmission, vol. 2A, pp. 61–74, 1948. [WAR 97] WARD J.P., Quaternions and Cayley Numbers: Algebra and Applications, Kluwer, Dordrecht, vol. 403, 1997. [WEI 14a] W EISSTEIN E.W., “Cross-correlation”, MathWorld. Available at Wolfram Web Resource, http://mathworld.wolfram.com/Cross-Correlation.html, January 2014. [WEI 14b] W EISSTEIN E.W., “Cross-correlation theorem”, MathWorld. Available at Wolfram Web Resource, http://mathworld.wolfram.com/Cross-CorrelationTheorem.html, January 2014. [YEH 08] Y EH M.-H., “Relationships among various 2-D quaternion Fourier transforms”, Signal Processing Letters, IEEE, vol. 15, pp. 669–672, 2008. [ZHA 97] Z HANG F.Z., “Quaternions and matrices of quaternions”, Linear Algebra and its Applications, vol. 251, pp. 21–57, January 1997.

Index

A absolutely integrable, 40 algebra Clifford, 65 division, 1, 3, 5 amplitude instantaneous, 101 analytic signal, 96, 100 anti-involution, 4, 23 arc, spherical, 27 associativity, 3 axis, 11 B Baker-Campbell-Hausdorff formula, 9 basis, 2 change of, 61 bicomplex product commutativity, 16 binormal vector, 103 bivariate signal analysis, 91 Boost, 57 C Cayley-Dickson form, 12 polar, 14 change of basis, 61 change of coordinates, 2D, 52 Cj -pair, 15 color image, 70

decomposition, 75 color space RGB, 70 commutativity non-, 3 complex exponential, 82 complex signal instantaneous amplitude, 91 instantaneous phase, 91 properties of 1D QFT, 91 complex subﬁeld, 86 conjugate, 4 time/frequency reversal, 45 conjugation quaternion, 23 convolution, 67 classical, 2D, 67 complex, 36, 96 mask, 68 quaternion, 70 theorem, 47, 68 coordinates change of, 52 homogeneous, 28 correlation classical, 81 color and grayscale image, 91 complex, 36 frequency domain, 86 peak interpretation, 88

124

Quaternion Fourier Transforms for Signal and Image Processing

phase, 80–82, 88 sub-sample precision, 84 quaternion, 79, 86 theorem, 47 cross product, 3 cross-covariance, 81 cross-power spectrum, 80 D de Moivre, 18 decomposition symplectic, 45 decompositions, 42 2D, 54 derivatives, 42, 53 dilation, 3D, 24 dirac delta function, 42 dirichlet conditions, 39 division algebra 1, 3, 5 E edge detector, 68 envelope, 103 complex, 105, 107 Euler angles, 12 formula, 8, 11, 64 exponential complex, 82 matrix, 63 quaternion, 7, 9 F FFTW library, 61 Fourier convolution theorem, 68 transform coefﬁcients, 65 complex, 35, 37 computation by decomposition, 60 direct DFT coding, 58 direct FFT coding, 60 discrete-time, 65 dual-axis form, 48 factored form, 48 fast, 60 FFT, 60

matrix exponential, 63 pairs, 1D, 40 pairs, 2D, 56 pairs, complex, 37 quaternion, 1D, 38 quaternion, 2D, 48 tiling, 82 veriﬁcation, 62 transforms Clifford coding, 57 computation of, 57 Frenet-Serret formulas, 103 frequency instantaneous, 91 mean Wigner-Ville distribution, 113 negative, 66, 80, 97 Nyquist, 66, 75 frequency domain, 81 G Gabor, D., 96 geometry, 21 Euclidean, 21 spherical, 26 great circle, 27 group general linear, 32, 33 group velocity, 91 H Hamilton, Sir W. R., 1 Hilbert transform, 96, 102 homogeneous coordinates, 28 I image color, 70 grayscale, 67 registration, 91 impulse, 42, 53 instantaneous amplitude, 91, 101 frequency, 91 of a complex signal, 102

Index

phase, 91, 101 integrable absolutely, 40 involution, 4, 5, 15 isometry, 33 K, L, M Kirsch, 68 linearity, 40, 52 logarithm quaternion, 7 marginal processing, 84 MATLAB®, 57, 64 matrix exponential, 63 representations, 17 mean frequency Wigner-Ville distribution, 113 modulated elliptical signal, 107 modulation orthocomplex, 95, 112 modulo arithmetic, 82 modulus, 4 N, O norm 1 product of, 3 quaternion, 3 Nyquist frequency, 66, 75 one-sided spectrum, 97 orthocomplex modulation, 95, 112 orthogonal split, 2D, 45 osculating plane, 103 P phase Fourier coefﬁcients, 80 instantaneous, 101 phase correlation, 80, 82 classical, 81 phase velocity, 91 plane, osculating, 103 point, weighted, 28 polar form, 101 Cayley-Dickson, 102 Prewitt, 68

product scalar, 3 vector, 3 projective space, 29 (3D), 28 Q QTFM, 57 quadrature, 105 signal, 96 quaternion axis, 11 basis, 2, 6 Cartesian form, 1 Cayley-Dickson form, 12 polar, 14 Cj -pair, 15 conjugate, 4 convolution, 70 de Moivre formula, 18 exponential, 7, 9 Fourier transform, see Fourier functions linear, 31 ijk, 1 inverse, 5 involution, 5, 15 properties, 5 logarithm, 7, 10 modulus, 4 norm, 3, 4 perplex part, 13 polar form, 8, 11 Euler angles, 12 power, 18 product, 3 bicomplex, 15 matrix-vector form, 32 pure, 2 ratio of, 5 representations, 11 matrix, 17 scalar part, 1 simplex part, 13 subﬁelds, 18 sum, 2

125

126

Quaternion Fourier Transforms for Signal and Image Processing

symplectic form, 13 systems of equations, 31 vector part, 1 R reﬂection 3D, 22 4D, 25 registration image, 91 reversal time, 80 time/frequency, 45, 46 implementation, 65 rotation, 71 3D, 22 4D, 25 composition, 23 S scalar part, 1 using conjugate, 4 scalar product, 3 scaling, 40, 52 shear, 24 shift frequency, 95 shifting, 40, 53 in correlation, 81 signal analytic, 100 complex instantaneous frequency, 102 non-stationary, 109 quaternion extension, 101 envelope, 103 quadrature, 96, 98 quaternion associated with complex, 98 signal analysis bivariate, 91 Sobel, 68 spectrum cross-power, 80 one-sided, 97 quaternion

layout, 75 visualization, 74 quaternion image, 73 spherical triangle, 28 split even-odd, 42 orthogonal 2D plane, 45 subﬁeld complex, 86 superposition, 33 swap rule, 13 symmetry, 42, 75, 80, 94, 100 complex signal, 93 QFT, 109 QFT of complex signal properties of, 98 symplectic decomposition, 45 symplectic form, 13 Synge, J.L., 81 T tangent vector, 103 tiling, 82 time domain, 80 time reversal, 80 time-frequency analysis, 91, 101 time/frequency reversal, 65 transform Hilbert, 96, 102 transformation afﬁne, 33 dilation, 24 Euclidean, 33 linear, 33 projective, 33 reﬂection 3D, 22 4D, 25 rotation, 71 3D, 22 4D, 25 shear, 24 similarity, 33 translation vector, 28 transmutation, 32 triangle, spherical, 28

Index

V

W, Z

vector translation, 28 vector part using conjugate, 4 versor, 6, 10 Ville, J., 96 visual system human, 70 visualization, quaternion spectrum, 74

Wigner-Ville distribution mean frequency, 113 quaternion, 109, 112 zero padding, 82

127

Figure 4.3. Filtered output (right) of color image (left) using Prewitt edge detector

Figure 4.4. Quaternion spectrum of an image. Left to right: original image, spectral modulus, phase and axis

Figure 4.5. Symplectic decomposition of color image into simplex and perplex parts of spectra. Left to right: original image, simplex part, perplex part and direct sum of both parts

Figure 4.7. Spectral coefﬁcient images for the constant term and ﬁrst three harmonics (horizontally and vertically) of color image. The DC term is the centre tile on bottom row with positive horizontal frequencies to the right and negative horizontal frequencies on the left. The bottom row corresponds to the DC vertical term with each previous row increasing in vertical frequency. (The left and right halves of the last row are mirrors of each other as they contain the same spectral coefﬁcient pairs.)

Figure 4.8. Harmonic reconstruction of color image. The legend at top shows a color-cube scatter plot for the original image. The lower grid shows plots of the color orbits for the ﬁrst ﬁve horizontal and vertical harmonics taken from the QFT spectrum of the image, along with their corresponding images. The ﬁnal, right-most column shows the image resulting from the sum of the harmonics, starting from the constant (zero frequency) term to the corresponding harmonic

Figure 4.9. Direct ﬁltering in frequency domain of color image. Top row, left to right: original image, low-pass, band-pass and high-pass ﬁltered. The bottom row shows the corresponding spectral magnitude

Figure 4.12. Two color images with identical backgrounds but a foreground object (the car) in different positions with different colors. (Images by and reproduced with the permission of Dr A.L. Thornton [THO 98])

Figure 4.13. Phase correlation results using the color images in Figure 4.12. Upper: the phase correlation surface. Lower: the two images superimposed after shifting one of them by the coordinate offset found from the maximum peak in the phase correlation surface, showing the two cars correctly superimposed

Figure 4.19. Orthocomplex modulated complex signal with constant frequency z(t). 3D representation of the complex signal z(t) (blue solid line) and its complex instantaneous amplitude ρz (t) (red solid line) obtained from the polar CD form of its quaternion extension qz (t). The projections in and planes (black solid lines) correspond to the signals presented in Figure 4.20

Figure 4.21. Orthocomplex modulated complex signal with linearly varying frequency z(t). 3D representation of the complex signal z(t) (blue solid line) and its complex instantaneous amplitude ρz (t) (red solid line) obtained from the polar CD form of its quaternion extension qz (t). The projections in and planes (black solid lines) correspond to the signals presented in Figure 4.22