U.S. patent application number 12/482237 was filed with the patent office on 2010-07-29 for signal processing with fast s-transforms.
This patent application is currently assigned to UTI LIMITED PARTNERSHIP. Invention is credited to Robert Brown, Richard Frayne, M. Louis Lauzon.
Application Number | 20100191792 12/482237 |
Document ID | / |
Family ID | 41417186 |
Filed Date | 2010-07-29 |
United States Patent
Application |
20100191792 |
Kind Code |
A1 |
Brown; Robert ; et
al. |
July 29, 2010 |
Signal Processing with Fast S-Transforms
Abstract
The ability to examine the frequency content of a signal is
critical in a variety of fields, and many techniques have been
proposed to fill this need, including the Fourier and wavelet
family of transforms. One of these, the S-transform, is a Fourier
based transform that provides simultaneous time and frequency
information similar to the wavelet transform but uses sinusoidal
basis functions to produce true frequency and globally referenced
phase measurements. It has been shown to be useful in several
medical imaging applications but its use is limited due to high
computational requirements of the original, continuous form. The
described embodiments include a general framework for describing
linear time-frequency transforms, using the Fourier, wavelet and
S-transforms as examples. As an illustration of the utility of this
formalism, a fast discrete S-transform algorithm is developed that
has the same computational complexity as the fast Fourier
transform.
Inventors: |
Brown; Robert; (Ille des
Soeurs, CA) ; Lauzon; M. Louis; (Calgary, CA)
; Frayne; Richard; (Calgary, CA) |
Correspondence
Address: |
FULBRIGHT & JAWORSKI L.L.P.
600 CONGRESS AVE., SUITE 2400
AUSTIN
TX
78701
US
|
Assignee: |
UTI LIMITED PARTNERSHIP
|
Family ID: |
41417186 |
Appl. No.: |
12/482237 |
Filed: |
June 10, 2009 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
61060315 |
Jun 10, 2008 |
|
|
|
Current U.S.
Class: |
708/404 ;
708/400 |
Current CPC
Class: |
G06F 17/148
20130101 |
Class at
Publication: |
708/404 ;
708/400 |
International
Class: |
G06F 17/14 20060101
G06F017/14 |
Claims
1. A method comprising performing a general time-frequency
transform on a dataset associated with a signal according to a fast
S-transform having a computational complexity of O(N log N), where
N is the total number of samples in the acquired-dataset.
2. The method of claim 1, further comprising performing the general
time-frequency transform on the dataset associated with the signal
according to a fast S-transform having a resource usage of
O(N).
3. The method of claim 1, where the dataset is associated with a
multi-dimensional signal.
4. The method of claim 1, further comprising: receiving a signal;
sampling the signal to generate an signal dataset; transforming the
signal dataset representation with a fast Fourier transform to
obtain the frequency-domain representation; dividing the
frequency-domain into sub-regions; multiplying each sub-region by a
window function; inverse transforming each sub-region individually
with an inverse fast Fourier transform; arranging the inverse
transformed sub-regions to form a representation of the frequency
content and changes in frequency content along each axis of the
signal dataset; and generating an output comprising the
representation of the time-frequency content.
5. The method of claim 4, further comprising dividing the frequency
domain into a plurality of overlapping sub-regions.
6. The method of claim 1, further comprising: receiving a signal;
sampling the signal to generate an signal dataset; transforming the
signal dataset with a fast Fourier transform to obtain a
frequency-domain representation; dividing the frequency-domain into
a high-frequency sub-region and a low-frequency sub-region;
generating a second dataset from an inverse fast Fourier transform
of the high-frequency sub-region; multiplying the second dataset by
a complex sinusoid of a predetermined frequency and with the same
number of dimensions as the signal dataset to form a third dataset;
selecting a first portion of the third dataset; multiplying the
selected portion of the third dataset by a window function; and
computing a sum of the results of multiplying the selected portion
of the third dataset by the window function to obtain the frequency
measurement at a predetermined sample location in the signal
dataset.
7. The method of claim 6, further comprising repeating steps with a
second selected portion of the third data set to obtain the
frequency measurements at one or more additional frequencies the
signal dataset.
8. The method of claim 7, further comprising: selecting the low
frequency sub-region from as the frequency-domain data; generating
a second dataset from an inverse fast Fourier transform of the
low-frequency sub-region; multiplying the second dataset by a
complex sinusoid of a predetermined frequency and with the same
number of dimensions as the signal dataset to form a fourth
dataset; selecting a first portion of the fourth dataset;
multiplying the selected portion of the fourth dataset by a window
function; computing a sum of the results of multiplying the
selected portion of the fourth dataset by the window function to
obtain the frequency measurement at a predetermined sample location
in the signal dataset; repeating steps with a second selected
portion of the fourth data set to obtain the frequency measurements
at one or more additional frequencies the signal dataset; and
arranging the resultant data to form a representation of the
frequency content and changes in frequency content along each axis
of the acquired dataset.
9. A system comprising a signal processor suitably configured to
perform a general time-frequency transform on a dataset associated
with a signal according to a fast S-transform having a
computational complexity of O(N log N) and resource usage of O(N),
where N is the total number of samples in the acquired-dataset.
10. The system of claim 9, where the signal processor is further
configured to performing the general time-frequency transform on
the dataset associated with the signal according to a fast
S-transform having a resource usage of O(N).
11. The system of claim 9, where the dataset is associated with a
multi-dimensional signal.
12. The system of claim 9, further comprising: a signal receiver
coupled to the signal processor, the signal receiver configured: to
receive a signal; and sample the signal to generate an signal
dataset.
13. The system of claim 12, the signal processor further configured
to: transform the signal dataset representation with a fast Fourier
transform to obtain the frequency-domain representation; divide the
frequency-domain into sub-regions; multiply each sub-region by a
window function; inverse transform each sub-region individually
with an inverse fast Fourier transform; arrange the inverse
transformed sub-regions to form a representation of the frequency
content and changes in frequency content along each axis of the
signal dataset; and generate an output comprising the
representation of the frequency content.
14. The system of claim 13, the signal processor further configured
to divide the frequency domain into a plurality of overlapping
sub-regions.
15. The system of claim 12, the signal processor further configured
to: transform the signal dataset with a fast Fourier transform to
obtain a frequency-domain representation; divide the
frequency-domain into a high-frequency sub-region and a
low-frequency sub-region; generate a second dataset from an inverse
fast Fourier transform of the high-frequency sub-region; multiply
the second dataset by a complex sinusoid of a predetermined
frequency and with the same number of dimensions as the signal
dataset to form a third dataset; select a first portion of the
third dataset; multiply the selected portion of the third dataset
by a window function; and compute a sum of the results of
multiplying the selected portion of the third dataset by the window
function to obtain the frequency measurement at a predetermined
sample location in the signal dataset.
16. The system of claim 15, the signal processor further configured
to repeat the steps of claim 10 with a second selected portion of
the third data set to obtain the frequency measurements at one or
more additional frequencies the signal dataset.
17. The system of claim 16, the signal processor further configured
to: select the low frequency sub-region from as the
frequency-domain data; generate a second dataset from an inverse
fast Fourier transform of the low-frequency sub-region; multiply
the second dataset by a complex sinusoid of a predetermined
frequency and with the same number of dimensions as the signal
dataset to form a fourth dataset; select a first portion of the
fourth dataset; multiply the selected portion of the fourth dataset
by a window function; and compute a sum of the results of
multiplying the selected portion of the fourth dataset by the
window function to obtain the frequency measurement at a
predetermined sample location in the signal dataset; repeat steps
with a second selected portion of the fourth data set to obtain the
frequency measurements at one or more additional frequencies the
signal dataset; and arrange the resultant data to form a
representation of the frequency content and changes in frequency
content along each axis of the acquired dataset.
Description
CROSS-REFERENCE TO RELATED APPLICATIONS
[0001] This application claims priority to U.S. Provisional
Application No. 61/060,315 filed Jun. 10, 2008, the entire text of
which disclosure is specifically incorporated by reference herein
without disclaimer.
BACKGROUND OF THE INVENTION
[0002] 1. Field of the Invention
[0003] The present invention relates generally to the fields of
signal processing. More particularly, it concerns methods and
systems for signal processing with fast S-transforms.
[0004] 2. Description of the Related Art
[0005] The ability to detect the frequency content of a signal is a
desirable capability in diverse applications from electrical
engineering to neuroscience. Many transforms for accomplishing this
task by transforming signals of one or more dimensions into
frequency or frequency-analogue spaces are known to those of skill
in the art.
[0006] While the original Fourier transform (FT) is an extremely
important signal and image analysis tool, it assumes that a signal
is stationary, i.e., that the frequency content is constant at all
times in a signal, or at all locations in an image. Since most
interesting signals are non-stationary, a series of techniques have
been developed to characterize signals with dynamic frequency
content. These are the foundation of the field of time-frequency
analysis.
[0007] A simple approach to the problem of non-stationary signal
analysis is the short-time Fourier transform (STFT). In this
technique changes in frequency over time are captured using a
window function to provide time localization. However, the choice
of window function represents a compromise. Narrower windows
provide better time resolution but poorer frequency resolution,
while wider windows provide the converse. Ideally, the window width
is chosen to produce the best representation of particular features
of interest in the signal, but this requires a priori
knowledge.
[0008] The wavelet transform (WT), which has been applied to a wide
variety of signal processing problems, improves on the STFT by
introducing the concept of progressive resolution. The WT provides
the equivalent of finer time resolution at high frequencies and
finer frequency resolution at low frequencies. However, the WT does
not measure frequency but only an analogue, called scale.
Additionally, the WT provides either no phase information, or phase
measurements which are all relative to different local references.
This is in contrast to the conventional concept of phase, as
provided by the FT, where all phase measurements are relative to a
global reference.
[0009] The S-transform (ST) exhibits the globally referenced phase
and true frequency measurements of the FT and STFT, as well as the
progressive resolution of the WT. This combination of desirable
features have shown promise in a wide variety of applications,
including, for example, detecting abnormalities in the heart,
identifying genetic abnormalities in brain tumors, analyzing
electroencephalograms, transmitting medical images, characterizing
the behavior of liquid crystals, detecting disturbances in
electrical power distribution networks monitoring high altitude
wind patterns, and detecting gravitational waves. However, while
very efficient methods have been developed for calculating the FT
and WT, the computational demands of the ST have limited its
utility.
SUMMARY OF THE INVENTION
[0010] From the foregoing discussion, it should be apparent that a
need exists for systems and methods for signal processing with fast
S-transforms.
[0011] A method for signal processing with fast S-transforms is
presented. In one embodiment, the method includes performing a
general time-frequency transform on a dataset associated with a
signal according to a fast S-transform having a computational
complexity of O(N log N), where N is the total number of samples in
the acquired-dataset. In a further embodiment, the method may
include performing the general time-frequency transform on the
dataset associated with the signal according to a fast S-transform
having a resource usage of O(N). In still a further embodiment of
the method, the dataset is associated with a multi-dimensional
signal.
[0012] In one embodiment, the method includes receiving a signal,
sampling the signal to generate an signal dataset, transforming the
signal dataset representation with a fast Fourier transform to
obtain the frequency-domain representation, dividing the
frequency-domain into sub-regions, multiplying each sub-region by a
window function, inverse transforming each sub-region individually
with an inverse fast Fourier transform, arranging the inverse
transformed sub-regions to form a representation of the frequency
content and changes in frequency content along each axis of the
signal dataset, and generating an output comprising the
representation of the time-frequency content. The method may
further include dividing the frequency domain into a plurality of
overlapping sub-regions.
[0013] In another embodiment, the method may include receiving a
signal, sampling the signal to generate an signal dataset,
transforming the signal dataset with a fast Fourier transform to
obtain a frequency-domain representation, dividing the
frequency-domain into a high-frequency sub-region and a
low-frequency sub-region, generating a second dataset from an
inverse fast Fourier transform of the high-frequency sub-region,
multiplying the second dataset by a complex sinusoid of a
predetermined frequency and with the same number of dimensions as
the signal dataset to form a third dataset, selecting a first
portion of the third dataset, multiplying the selected portion of
the third dataset by a window function, and computing a sum of the
results of multiplying the selected portion of the third dataset by
the window function to obtain the frequency measurement at a
predetermined sample location in the signal dataset. The method may
further include repeating these steps with a second selected
portion of the third data set to obtain the frequency measurements
at one or more additional frequencies the signal dataset.
[0014] In one embodiment, the method may include selecting the low
frequency sub-region from as the frequency-domain data, generating
a second dataset from an inverse fast Fourier transform of the
low-frequency sub-region, multiplying the second dataset by a
complex sinusoid of a predetermined frequency and with the same
number of dimensions as the signal dataset to form a fourth
dataset, selecting a first portion of the fourth dataset,
multiplying the selected portion of the fourth dataset by a window
function, and computing a sum of the results of multiplying the
selected portion of the fourth dataset by the window function to
obtain the frequency measurement at a predetermined sample location
in the signal dataset, repeating steps with a second selected
portion of the fourth data set to obtain the frequency measurements
at one or more additional frequencies the signal dataset, and
arranging the resultant data to form a representation of the
frequency content and changes in frequency content along each axis
of the acquired dataset.
[0015] A system for signal processing with fast S-transforms is
also presented. In one embodiment, the system includes a signal
processor suitably configured to perform a general time-frequency
transform on a dataset associated with a signal according to a fast
S-transform having a computational complexity of O(N log N) and
resource usage of O(N), where N is the total number of samples in
the acquired-dataset. In a further embodiment, signal processor is
further configured to performing the general time-frequency
transform on the dataset associated with the signal according to a
fast S-transform having a resource usage of O(N). Additionally, the
dataset is associated with a multi-dimensional signal.
[0016] In a further embodiment, the system may include a receiver
coupled to the signal processor. The receiver may receive a signal,
and sample the signal to generate an signal dataset.
[0017] In one embodiment, the signal processor may be configured to
transform the signal dataset representation with a fast Fourier
transform to obtain the frequency-domain representation, divide the
frequency-domain into sub-regions, multiply each sub-region by a
window function, inverse transform each sub-region individually
with an inverse fast Fourier transform, arrange the inverse
transformed sub-regions to form a representation of the frequency
content and changes in frequency content along each axis of the
signal dataset, and generate an output comprising the
representation of the frequency content. Additionally, the signal
processor may divide the frequency domain into a plurality of
overlapping sub-regions.
[0018] In another embodiment, the signal processor may be
configured to transform the signal dataset with a fast Fourier
transform to obtain a frequency-domain representation, divide the
frequency-domain into a high-frequency sub-region and a
low-frequency sub-region, generate a second dataset from an inverse
fast Fourier transform of the high-frequency sub-region, multiply
the second dataset by a complex sinusoid of a predetermined
frequency and with the same number of dimensions as the signal
dataset to form a third dataset, select a first portion of the
third dataset, multiply the selected portion of the third dataset
by a window function, and compute a sum of the results of
multiplying the selected portion of the third dataset by the window
function to obtain the frequency measurement at a predetermined
sample location in the signal dataset. In such an embodiment, the
signal processor may repeat these steps with a second selected
portion of the third data set to obtain the frequency measurements
at one or more additional frequencies the signal dataset.
[0019] In still a further embodiment, the signal processor may
select the low frequency sub-region from as the frequency-domain
data, generate a second dataset from an inverse fast Fourier
transform of the low-frequency sub-region, multiply the second
dataset by a complex sinusoid of a predetermined frequency and with
the same number of dimensions as the signal dataset to form a
fourth dataset, select a first portion of the fourth dataset,
multiply the selected portion of the fourth dataset by a window
function, and compute a sum of the results of multiplying the
selected portion of the fourth dataset by the window function to
obtain the frequency measurement at a predetermined sample location
in the signal dataset, repeat steps with a second selected portion
of the fourth data set to obtain the frequency measurements at one
or more additional frequencies the signal dataset, and arrange the
resultant data to form a representation of the frequency content
and changes in frequency content along each axis of the acquired
dataset.
[0020] The term "coupled" is defined as connected, although not
necessarily directly, and not necessarily mechanically.
[0021] The terms "a" and "an" are defined as one or more unless
this disclosure explicitly requires otherwise.
[0022] The term "substantially" and its variations are defined as
being largely but not necessarily wholly what is specified as
understood by one of ordinary skill in the art, and in one
non-limiting embodiment "substantially" refers to ranges within
10%, preferably within 5%, more preferably within 1%, and most
preferably within 0.5% of what is specified.
[0023] The terms "comprise" (and any form of comprise, such as
"comprises" and "comprising"), "have" (and any form of have, such
as "has" and "having"), "include" (and any form of include, such as
"includes" and "including") and "contain" (and any form of contain,
such as "contains" and "containing") are open-ended linking verbs.
As a result, a method or device that "comprises," "has," "includes"
or "contains" one or more steps or elements possesses those one or
more steps or elements, but is not limited to possessing only those
one or more elements. Likewise, a step of a method or an element of
a device that "comprises," "has," "includes" or "contains" one or
more features possesses those one or more features, but is not
limited to possessing only those one or more features. Furthermore,
a device or structure that is configured in a certain way is
configured in at least that way, but may also be configured in ways
that are not listed.
[0024] Other features and associated advantages will become
apparent with reference to the following detailed description of
specific embodiments in connection with the accompanying
drawings.
BRIEF DESCRIPTION OF THE FIGURES
[0025] The following drawings form part of the present
specification and are included to further demonstrate certain
aspects of the present invention. The invention may be better
understood by reference to one or more of these drawings in
combination with the detailed description of specific embodiments
presented herein.
[0026] FIG. 1 is a schematic block diagram illustrating one
embodiment of a system for signal processing with fast
S-transforms;
[0027] FIG. 2 is a schematic block diagram illustrating one
embodiment of a computing device suitably configured for signal
processing with fast S-transforms;
[0028] FIG. 3 is a schematic block diagram illustrating another
embodiment of a signal processing device configured for signal
processing with fast S-transforms;
[0029] FIGS. 4A-F show an example of the construction of the
.alpha.-domain for a simple signal, and the associated
S-domain;
[0030] FIGS. 5A-D illustrate one embodiment of a band-passed region
of the signal shown in FIG. 4A;
[0031] FIGS. 6A-C illustrate one embodiment sampled and shifted
signals in the .alpha.-domain;
[0032] FIG. 7 illustrates one embodiment of benchmarks comparing
the computation speed of the FFT (dotted black line), DST (solid
black line) and fast GFT (gray line) for one-dimensional signals of
various lengths;
[0033] FIGS. 8A-C illustrate one embodiment of an application of
the fast S-transform to an ECG signal;
[0034] FIGS. 9A-D illustrate one embodiment of a .nu.-.tau.
sampling scheme in accordance with the present embodiments;
[0035] FIGS. 10A-B illustrate embodiments of sampling schemes for
the discrete S-transform of a complex 8 sample signal;
[0036] FIGS. 11A-C illustrate another embodiment of a fast
S-transform of a sample signal in accordance with the present
embodiments;
[0037] FIGS. 12A-D illustrate embodiments of S-transform kernels
and windows in accordance with the present embodiments;
[0038] FIGS. 13A-H illustrated embodiments of discrete S-transform
window-kernel product in the .nu.-domain for .nu.=4 in accordance
with the present embodiments;
[0039] FIGS. 14A-D illustrate a comparison of the ST and FST in
accordance with the present embodiments;
[0040] FIG. 15 illustrates one embodiment of a spatial frequency
spectra calculated according to the present embodiments;
[0041] FIG. 16 illustrates one embodiment of a spatial frequency
spectra indicating a difference in the texture of CT images of
trabecular bone from normal controls and osteoporosis patients
calculated according to the present embodiments; and
[0042] FIG. 17 illustrates one embodiment of a Fast GFT spectrum of
an audio recording of a volunteer speaking the English vowel sounds
a, e, i, o and u calculated according to the present
embodiments.
DETAILED DESCRIPTION
[0043] Various features and advantageous details are explained more
fully with reference to the nonlimiting embodiments that are
illustrated in the accompanying drawings and detailed in the
following description. Descriptions of well known starting
materials, processing techniques, components, and equipment are
omitted so as not to unnecessarily obscure the invention in detail.
It should be understood, however, that the detailed description and
the specific examples, while indicating embodiments of the
invention, are given by way of illustration only, and not by way of
limitation. Various substitutions, modifications, additions, and/or
rearrangements within the spirit and/or scope of the underlying
inventive concept will become apparent to those skilled in the art
from this disclosure.
[0044] The present embodiments describe a general Fourier Family
transform (GFT), which may describe linear transforms based on the
complex sinusoid, including the FT, STFT and ST. Additionally, a
fast algorithm for calculating the GFT is described that is much
more computationally efficient than existing methods for
calculating the ST, making the computation of the ST more
practical. This greatly improved computational efficiency will
improve existing techniques as well as allowing more advanced
applications than have previously been investigated. Certain
benefits of the described embodiments include the ability to
analyze longer signals and larger images, as well as transform data
with three or more dimensions, e.g., volumetric data from
tomographic medical imaging techniques. As described below, the
present embodiments may be extended to higher dimensional data.
[0045] The terminology used in the literature may inconsistent
between the S-transform, wavelet and Fourier transforms. As used
herein, the term "continuous" refers to transforms expressed as
integrals. The term "discrete" refers to formulations using
summation. The term "one-to-one" refers to discrete transforms that
are minimally sampled: i.e., the transform of an N-point signal
produces an N-point spectrum. This may reduce the memory
requirements of the algorithm and also implies that each unique
signal is transformed into a unique spectrum, and vice versa. The
term "fast" transform means an algorithm for calculating a
discrete, one-to-one transform that eliminates remaining
redundancies to achieve a reduction in computational
complexity.
[0046] The present embodiments may involve several transforms, each
of which having several variants. The continuous Fourier transform
is commonly used in analytic time-frequency analysis. The discrete
Fourier transform (DFT) may include one-to-one transform.
Alternatively, a fast Fourier transform (FFT) may be implemented
which offers O(N log N) computational complexity as compared to
O(N.sup.2) complexity for the DFT. The short-time Fourier transform
(STFT) may have continuous and discrete formulations. The latter
may be either one-to-one or not, and fast algorithms may exist to
calculate the STFT efficiently. Though there is a true continuous
wavelet transform, the term is often used to describe a discrete
approximation. In the present embodiments, the wavelet transform
refers to a discrete wavelet transform (DWT), which may not be
one-to-one. What is commonly called the discrete wavelet transform
is specifically referred to as the one-to-one discrete wavelet
transform in the present embodiments.
[0047] Finally, the S-transform may have a number of variations.
For example, the S-transform may include a continuous
time-frequency transform that uses Gaussian windows to provide time
localization. Alternatively, the discrete (but not one-to-one) ST
(DST), and the "fast" ST, which is a method to calculate the
discrete ST from the frequency domain and, though more efficient
than the time-domain version, is still quite slow when compared to
other time-frequency transforms. As used herein, the DST calculated
in the frequency domain is called the frequency-domain ST, or
frequency-domain DST.
[0048] Additional transforms include a discrete orthonomal
S-transform (DOST). The DOST refers to a discrete, though still not
one-to-one, transform that uses sinc shaped windows (boxcar
functions in the frequency domain). Since Gaussian windows are not
supported, the DOST is not really an S-transform.
[0049] The GFT of the present embodiments may include a general
transform that encompasses the FT, STFT, ST, GST and DOST. The fast
GFT algorithm described herein may be a fast, one-to-one, discrete
transform. Due to the generality of the GFT, this algorithm may
represent an efficient means of calculating the STFT, ST and GST,
as well as a faster method of computing the DOST. Besides unifying
these transforms in one framework, the fast GFT algorithm greatly
improves on the computational efficiency of the ST, GST and DOST. A
signal that would take days to transform using the frequency-domain
DST algorithm may be transformed using the fast GFT in less than a
second.
[0050] The present embodiments also describe a general framework
for describing linear time-frequency transforms that allows the
Fourier, wavelet and S-transforms to be compared directly. The fast
discrete S-transform algorithm described herein may greatly
diminished computational and storage demands. This transform is
much more computationally efficient than original continuous
approximations of the ST, and so, allows the ST to be used in acute
clinical situations as well as allowing more advanced applications
than have been investigated to date, including analyzing longer
signals and larger images, as well as transforming data with three
or more dimensions, e.g., volumes obtained by magnetic resonance
(MR) or computed tomography (CT) imaging.
System Architecture
[0051] FIG. 1 illustrates one embodiment of a system 100 configured
for signal processing with fast S-transforms. In one embodiment the
system 100 includes a signal receiver 102 and a signal processor
104. The system 100 may additionally include a signal transmitter
106. In a further embodiment, the system 100 may include one or
more user interface devices.
[0052] According to one embodiment, the signal receiver 102 may
include an input port and/or coupler (not shown). The input port
may receive an incoming signal from, e.g., an antenna, scanner,
sensor, or the like. In a further embodiment, the signal receiver
may include an analog to digital (A/D) converter configured to
sample the incoming signal. The A/D converter may generate a signal
dataset comprising values representing sampled signal levels.
[0053] In a further embodiment, the signal receiver 102 may be
coupled to the signal processor 104. The signal receiver 102 may
provide the signal dataset to the signal processor 104. The signal
processor 104 may be configured to compute a fast ST of the signal
dataset as described below. In one embodiment, the signal processor
104 may calculate the fast ST according to a frequency-domain
algorithm. Alternatively, the signal processor may compute the fast
ST according to a time-domain algorithm. Both the frequency-domain
algorithm and the time-domain algorithm are described in greater
detail below.
[0054] In a further embodiment, the signal processor 104 may
perform one or more operations on a resultant fast ST dataset. The
signal processor 104 may then convert the fast ST dataset to a
second signal dataset by using an inverse fast ST as described in
greater detail below. In one embodiment, the signal processor 104
may provide the second signal dataset to the signal transmitter
106. The signal transmitter 106 may either transmit an analogue or
a digital signal corresponding to the second signal dataset. In one
embodiment, the signal transmitter 106 may include a digital to
analogue (D/A) converter and one or more output ports and/or
couplers (not shown).
[0055] In a further embodiment, the signal processor 104 may be
coupled to a user interface device 108. A user interface device 108
may include, but is not limited to, a display device, a keyboard, a
mouse device, or the like. The signal processor 104 may communicate
information with the user interface device 108. For example, the
signal processor 104 may send information to the user interface
device 108 configured to generate a graphical user interface, an
image display, or the like. In a further embodiment, the signal
processor 104 may be configured to receive one or more user inputs,
such as control commands, from the user interface device 108.
[0056] FIG. 2 illustrates a computer system 200 adapted according
to certain embodiments of the signal processor 104. The central
processing unit (CPU) 202 may be coupled to the system bus 204. The
CPU 202 may be a general purpose CPU or microprocessor. The present
embodiments are not restricted by the architecture of the CPU 202
as long as the CPU 202 supports the modules and operations as
described herein. The CPU 202 may execute the various logical
instructions according to the present embodiments. For example, the
CPU 202 may execute machine-level instructions according to the
exemplary operations described below with reference to the
frequency-domain fast ST or the time-domain fast ST.
[0057] The computer system 200 also may include Random Access
Memory (RAM) 208, which may be SRAM, DRAM, SDRAM, or the like. The
computer system 200 may utilize RAM 208 to store the various data
structures used by a software application configured to perform the
fast ST. The computer system 200 may also include Read Only Memory
(ROM) 206 which may be PROM, EPROM, EEPROM, or the like. The ROM
may store configuration information for booting the computer system
200. The RAM 208 and the ROM 206 hold user and system 100 data.
[0058] The computer system 200 may also include an input/output
(I/O) adapter 210, a communications adapter 214, a user interface
adapter 216, and a display adapter 222. The I/O adapter 210 and/or
user the interface adapter 216 may, in certain embodiments, enable
a user to interact with the user interface device 108 in order to
input information for, e.g., setting parameters of the fast ST.
[0059] The I/O adapter 210 may connect to one or more storage
devices 212, such as one or more of a hard drive, a Compact Disk
(CD) drive, a floppy disk drive, a tape drive, to the computer
system 200. The storage device 212 may be a tangible computer
readable medium comprising computer readable instructions that,
when executed by the CPU 202, cause the computer system 200 to be
suitably programmed to perform signal processing with the fast ST
according to the present embodiments. The communications adapter
214 may be adapted to couple the computer system 200 to the network
106, which may be one or more of a LAN and/or WAN, and/or the
Internet. The user interface adapter 216 couples user input
devices, such as a keyboard 220 and a pointing device 218, to the
computer system 200. The display adapter 222 may be driven by the
CPU 202 to control the display on the display device 224.
[0060] The present embodiments are not limited to the architecture
of system 200. Rather the computer system 200 is provided as an
example of one type of computing device that may be adapted to
perform the functions of a signal processor 104. For example, any
suitable processor-based device may be utilized including without
limitation, including personal data assistants (PDAs), computer
game consoles, and multi-processor servers. Moreover, the present
embodiments may be implemented on application specific integrated
circuits (ASIC) or very large scale integrated (VLSI) circuits. In
fact, persons of ordinary skill in the art may utilize any number
of suitable structures capable of executing logical operations
according to the described embodiments.
[0061] For example, FIG. 3 is a schematic block diagram
illustrating one embodiment of logical modules of a single
processor 104. In one embodiment, the signal processor 104 may
include Digital Signal Processing (DSP) device that may be
configured to perform the fast ST according to the present
embodiments. For example, one embodiment of a DSP device that may
be used in accordance with the present embodiments is a Super
Harvard Architecture Single-Chip Computer (SHARC) type DSP. In one
embodiment, the DSP 104 may include an .alpha.-domain module 302
for performing the .alpha.-domain ST, a time-domain module 304 for
performing the time-domain ST, and a multi-dimensional ST module
306 for extending the fast ST to multi-dimensional data. Further
embodiments of the function of each module are described below with
reference to the .alpha.-domain ST, the time-domain ST, and the
multi-dimensional ST respectively.
.alpha.-Domain S-Transform
The Fourier Family of Transforms
[0062] The Fourier transform (FT) describes a signal in terms of a
series of complex sinusoids, with varying amplitude and phase. The
sinusoidal basis functions of the FT are purely periodic and
infinite in extent, and the FT converts entirely between the time
and frequency domains, with no direct temporal information
remaining after the transform. As noted in the introduction, this
process is ideal for stationary signals, where the frequency
content does not change over time. However, in many applications,
non-stationary signals are common, and frequency changes over time
or space are often not only present, but are in fact the very
feature the inventors are interested in examining. To allow
examination of non-stationary signals, a number of solutions have
been proposed, including the STFT and, more recently, the ST.
[0063] The ST of a time signal, s(t), may be expressed as:
S ( .tau. , v ) = .intg. - .infin. + .infin. s ( t ) v 2 .pi. - (
.tau. - t ) 2 v 2 2 - 2 .pi. vt t ( 1 ) ##EQU00001##
where s(t) is the signal to be transformed and .tau. and .nu. are
the transform time and frequency coordinates. This equation has the
same form as that of the FT, but adds a normalized area Gaussian
window for time localization. It is important to note that the
window also depends on .nu., causing the width of the window to
decrease with increasing frequency. This automatically adjusts the
ST window to provide a progressive trade off between time and
frequency resolution for each frequency, with finer frequency
resolution at low frequencies and better time resolution at high
frequencies. Like the FT and STFT, the ST produces a complex
spectrum that includes both frequency and globally referenced phase
information.
[0064] The ST may be similar in form to the Morlet wavelet. By way
of comparison, the ST basis functions are not a proper wavelet as
they may not satisfy the condition of zero mean for admissible
wavelets. Additionally, the basis functions of the ST, and more
generally of the GFT, may not be shifted and scaled versions of a
mother function, and are not self-similar. Due to this, the GFT may
produce a globally referenced phase spectrum where all phase
measurements are relative to the origin. Though complex wavelets
exist, these typically produce locally referenced phase
measurements where the reference point changes as the wavelet is
shifted. The globally-referenced phase of the GFT allows the
Fourier spectrum to be directly recovered from any GFT, and
facilitates the use of phase analysis techniques such as
cross-spectral analysis.
[0065] In one embodiment, Magnetic Resonance (MR) imaging may
directly measure the complex Fourier spectrum of the object being
imaged. Many MR techniques, such as phase contrast MR, encode
velocity information in the phase. Furthermore, many artifacts in
MR imaging, such as those caused by motion and magnetic field
non-uniformities, manifest themselves as changes in the phase.
Since phase is critically important in this modality, many
filtering, reconstruction and analysis techniques are based on the
FT. These techniques could potentially be improved by using the ST,
since it offers the capability to perform local analysis of
non-stationary signals.
[0066] In one embodiment, the ST may be generalized to include
arbitrary windows as follows:
S ( .tau. , v ) = .intg. - .infin. + .infin. s ( t ) w ( .tau. - t
, v ) - 2 .pi. vt t ( 2 ) ##EQU00002##
where w(.tau.-t,.nu.) is a general window function that
satisfies:
.intg. - .infin. + .infin. w ( .tau. - t , v ) .tau. = 1 for all v
. ( 3 ) ##EQU00003##
[0067] The original ST may be a special case of this general
S-transform, where w(.tau.-t,.nu.) is a Gaussian function.
The General Fourier Family Transform
[0068] The ST may be effectively further generalized into a general
Fourier family transform (GFT), which may encompass all Fourier
family transforms as follows:
S ( .tau. , v ) = .intg. - .infin. + .infin. s ( t ) w ( t - .tau.
, .sigma. ) - 2 .pi. vt t . ( 4 ) ##EQU00004##
If .sigma.=.nu. (i.e., the windows are adaptive), the GFT is a
general ST or, if .sigma.=C, it is a short time Fourier transform
(STFT) with a constant window width C. The Fourier transform itself
can be understood in this context by using a constant window,
w(.tau.-.tau.,.sigma.)=1.
[0069] In the GFT formalism the complex sinusoid,
e.sup.-i2.pi..nu.t, may be labeled as the transform kernel, and
w(.sigma.-.tau.,.sigma.) as the transform window. The product of
the two forms the transform basis functions. In this context, the
FT, STFT, ST and GST are examples of the GFT, differing by
properties of their respective windows.
Discrete Fourier Family Transforms
[0070] Calculation of the discrete ST, or of the discrete GST,
presents a considerable computational challenge. Calculated from
the time domain according to Eq. (1), the DST of an N-point signal
generates an N.times.N-point spectrum. Each of those N.sup.2 points
requires a numerical approximation of the integral in Eq. (1),
which is a summation over the N-point signal multiplied by the
basis function. Therefore, the total computational complexity of
this operation is O(N.sup.3) and the resulting spectrum requires
O(N.sup.2) storage space. Alternatively, the DST may be calculated
somewhat more efficiently by operating on the Fourier spectrum and
taking advantage of the FFT algorithm. This alternative requires
O(N.sup.2 log N) operations and has the same memory requirements,
O(N.sup.2). In comparison, the DFT (and some discrete wavelet
transforms) can be calculated in place, with the transform
requiring no more space than the original signal. The
one-dimensional one-to-one DFT has a computational complexity of
O(N.sup.2). The FFT improves on this, with a computational
complexity of O(N log N).
[0071] Due to the high computational complexity of the ST,
transforms of longer signals and even small images has proven to be
challenging to compute and store. As an example, the DST of a
512.times.512 pixel image, which is a typical image size
encountered in medical imaging applications, may currently take
more than a day to compute and require around 1 TB of storage
(double precision complex floating point), compared to the FFT,
which currently takes milliseconds and requires 4 MB of
storage.
[0072] Though parallel and distributed computing techniques for
calculating the ST on multiprocessor machines and computing
clusters have been developed, computation of the DST for large
signals or images is prohibitive, severely limiting potential
applications of the ST. A more efficient fast DST algorithm may
both allow (i) implementation of ST-based techniques that have
shown promise, and (ii) application of the ST in many other areas
that benefit from signal processing.
Characterization of the Generalized Fourier Family Transform
[0073] Since any operation on a signal implies a concomitant
operation on the frequency spectrum of the signal, it is often
useful to consider such modifications simultaneously from both
perspectives. For example, the convolution of two signals is
equivalent to multiplication of the Fourier transforms of those
signals: convolution in the "time domain" is equivalent to
multiplication in the "frequency domain." In one embodiment, it may
be possible to make use of two more domains: the "(generalized)
S-domain," S(.sigma., .nu.), defined by Eq. (4), and the
".alpha.-domain," defined as:
.alpha. ( v ' , v ) = .intg. - .infin. + .infin. S ( .tau. , v ) -
2 .pi. v ' .tau. .tau. ( 5 ) ##EQU00005##
where the .nu.'-axis is produced by taking the Fourier transform
along the .tau.-axis of the S-domain.
[0074] The .alpha.-domain is also related to the frequency domain
by:
.alpha.(.nu.',.nu.)=F(.nu.'+.nu.)W(.nu.',.sigma.) (6)
as shown in the Appendix, where F(.nu.') and W(.nu.',.sigma.) are
the FTs of the signal and the window function along the
.sigma.-axis, respectively. The S-domain is recoverable from the
.alpha.-domain via:
S ( .tau. , v ) = .intg. - .infin. + .infin. .alpha. ( v ' , v ) 2
.pi. v ' .tau. v ' ( 7 ) ##EQU00006##
[0075] In one embodiment, a frequency domain DST algorithm may be
calculated via the .alpha.-domain: the Fourier transform of the
signal may be shifted and multiplied by the Fourier transform of
the window, then inverse Fourier transformed to build up the
S-domain. FIG. 4 shows an example of the construction of the
.alpha.-domain for a simple signal, and the associated S-domain. In
the described embodiment, the signal is a simple boxcar signal as
shown in FIG. 4A that has a Fourier transform which is a sinc
depicted in FIG. 4B. As defined in Eq. (6), the .alpha.-domain may
include shifted copies of the FT shown in FIG. 4C multiplied by the
FT of the window functions illustrated in FIG. 4D. The product
shown in FIG. 4E can be Fourier transformed along the .nu.'-axis to
produce the S-domain as shown in FIG. 4F. This example uses
Gaussian windows that scale with frequency, and so the result in
FIG. 4F is a discrete S-transform.
[0076] Since the DST produces an N.times.N result from an N-point
signal, it may contain redundant information. In the S-domain that
redundancy is illuminated by the Heisenberg uncertainty principle:
the more finely energy is localized in the time domain, the more
uncertain is its location in the frequency domain, and vice versa.
This basic property is not reflected by the N.times.N S-domain
sampling scheme since it produces uniform time (.tau.) and
frequency (.nu.) resolution at all points. However, it is not
entirely clear how to perform efficient sampling in the S-domain,
particularly for arbitrary windows. Though the DOST attempts this
using a dyadic sampling scheme, it still produces an S-domain with
more than N points. Furthermore, the DOST is restricted to windows
that are a boxcar function in the frequency domain, or
equivalently, a sinc function in the time domain, which has several
undesirable properties.
[0077] The redundancy of the N.times.N S-domain is more clearly
revealed in the corresponding N.times.N .alpha.-domain (FIG. 4).
Since each horizontal line (along .nu.') at a particular .nu. value
is a shifted version of the FT of the signal multiplied by
W(.nu.',.sigma.), the full frequency domain, and thus the signal,
is recoverable from any one of these horizontal lines, provided
W(.nu.',.sigma.) does not have any zero values. Similarly, each
vertical line at a particular .nu.'-value is also a shifted version
of F(.nu.) windowed by the known W(.nu.',.sigma.). In one
embodiment, the full Fourier spectrum may be obtained from any
.tau.-line of the S-domain by Fourier transforming the line,
dividing by W(.nu.',.sigma.), and shifting accordingly.
Additionally, by the Fourier Projection Slice Theorem, a projection
through the S-domain along the .tau.-direction produces the .nu.'=0
line of the .alpha.-domain, which is also a complete Fourier
spectrum.
The Fast General Fourier Family Transform
[0078] The .alpha.-domain representation and the Heisenberg
uncertainty principle clearly suggest two optimizations that can be
made to the N.times.N sampling scheme: (i) band pass filtering
along .nu.' and (ii) octave sampling along .nu.. The .alpha.-domain
is weighted by the Fourier transform of the window functions. The
purpose of the window is to provide localization in the time
domain, but via the uncertainty principle it also controls the
degree of localization in the frequency domain. To provide
localization in the frequency domain, the window may act as a
band-pass filter, i.e., frequencies at a distance from the window
centre will be attenuated. This is illustrated in FIG. 4D for the
Gaussian windows of the ST. A degree of optimization can be
achieved by discarding the portions of each .nu.'=line that are not
band-passed by the window. In transforms where .sigma. (the width
of the window) varies, the extent of the remaining .nu.'-segment
will vary to match the band-pass characteristics of the window
(FIG. 5B).
[0079] Further optimization is obtained by sampling along the
.nu.-axis in an octave fashion (FIG. 5C). Taken together, these two
optimizations produce a dyadic sampling scheme (FIG. 5C). The
result can be Fourier transformed to produce an N-point S-domain
(FIG. 5D).
[0080] By reversing the shift applied to the frequency domain as
part of the construction of the .alpha.-domain, the relationship of
this sampling scheme to the Fourier transform of the signal may
become clear (FIG. 6). The frequency domain may be partitioned and
windowed (FIG. 6C), then inverse Fourier transformed in parts to
produce the S-domain. This calculation is described in the pseudo
code for Algorithm 1, which is a fast GFT, synthesizing the
S-domain from the frequency domain.
Algorithm 1
[0081] The in-place fast general Fourier family transform is
described below. By way of illustration, A[*] is an array of N
points. The subscript on A may indicate its domain: t is the time
domain, .nu. is the frequency domain and .tau. represents the
S-domain samples along the .tau.-axis. In this example a dyadic
frequency sampling scheme may be used, but this could be adapted,
for example in the STFT.
TABLE-US-00001 LET A.sub.t[0:N] be an array of N samples
representing a signal 1. Calculate the FFT of the signal:
A.sub.v[0:N] = FFT.sub.t.fwdarw.v(A.sub.t[0:N]) FOR f start in 0 ,
1 , 2 , 4 , , N 8 , N 4 DO : ##EQU00007## IF f.sub.start < 2:
LET f.sub.width = 1, f.sub.centre = f.sub.start + 1 ELSE: LET f
width = f start , f centre = f start + f width 2 + 1 ##EQU00008##
LET f.sub.end = f.sub.start + f.sub.width LET f.sub.start' =
-f.sub.start, f.sub.end' = -f.sub.end and f.sub.centre' =
-f.sub.centre 2. Calculate the window function and FFT: W[0:N] =
FFT.sub.t.fwdarw.v (w[0:N]) 3. Shift the window function to be
centered on the frequency band: W[0:N] = shift (W[0:N],
f.sub.centre) 4. Apply the window to the frequency band:
A.sub.v[f.sub.start:f.sub.end] = A.sub.v[f.sub.start:f.sub.end]
W[f.sub.start:f.sub.end] 5. Calculate the in-place inverse FFT of
the frequency band: A.sub..tau.[f.sub.start:f.sub.end] =
FFT.sup.-1.sub.v.fwdarw..tau.(A.sub.v[f.sub.start:f.sub.end]) Note:
The result of this step, A.sub..tau.[f.sub.start:f.sub.end], maps
to S[0:.tau..sub.max, v.sub.centre] in the S-domain. 6. Repeat
steps 3 to 5 substituting f.sub.start', f.sub.end' and
f.sub.centre' to calculate the negative frequencies. END FOR
Computational Complexity of Algorithm 1
[0082] For a signal of N points, Algorithm 1 consists of an N point
FFT (step 1), a shift (step 3), an M point multiply (step 4) and an
M point inverse FFT (step 5). Steps 3 through 5 are performed
twice, once for the positive frequencies and once for the negative
frequencies (step 6). Although Algorithm 1 shows the window
functions calculated inside the loop (step 2) for brevity, they may
be pre-calculated.
[0083] Of the steps in the loop, the one with the highest order is
step 5, an FFT of M points, which can be accomplished in O(M log
M). For a dyadic sampling scheme, as in Algorithm 1, M takes values
in
0 , 1 , 2 , 4 , , N 8 , N 4 , ##EQU00009##
but all one-to-one GFT sampling schemes have
i = 1 K M i = N . ##EQU00010##
If the inventors let M.sub.i=f.sub.iN, where f.sub.i is the
fraction of points, i.e. f.sub.i.epsilon.(0,1), then
i = 1 K M i = N i = 1 K f i N = N i = 1 K f i = N , .thrfore. i = 1
K f i = 1. ##EQU00011##
[0084] Using this relation the inventors can write the order of
step 5 as
i = 1 K M i log ( M i ) . ##EQU00012##
Substituting M.sub.i=f.sub.iN the inventors have:
i = 1 K f i N log ( f i N ) = N i = 1 K f i [ log ( f i ) + log ( N
) ] = N log ( N ) i = 1 K f i + N i = 1 K f i log ( f i ) .
##EQU00013##
From above,
i = 1 K f i = 1 , ##EQU00014##
so this simplifies to
N log ( N ) + N i = 1 K f i log ( f i ) . ##EQU00015##
Since
[0085] f i .di-elect cons. ( 0 , 1 ) , i = 1 K f i log ( f i )
.ltoreq. 0. ##EQU00016##
In fact, it is zero only if K=1, which implies f.sub.i=1, the
trivial case corresponding to the FFT. For K>1 the inventors
have
i = 1 K f i log ( f i ) < 0 so i = 1 K M i log ( M i ) < N
log ( N ) . ##EQU00017##
[0086] Since the complexity of step 5 is less than N log N, the
complexity of the loop is dominated by the N point FFT in step 1,
which gives the order of the algorithm: O(N log N).
[0087] As shown above, Algorithm 1 has a computational complexity
of O(N log N) and transforms an N-point signal into an N-point
S-domain. This performance is a great improvement over (i) the time
domain DST, which is O(N.sup.3) and produces an N.sup.2-point
S-domain, (ii) the frequency domain DST, which is O(N.sup.2 log N)
and also produces an N.sup.2-point S-domain, and (iii) the DOST,
which is O(N.sup.2) and produces a CN-point S-domain, where
1<C.ltoreq.N. The computational complexity fast GFT is of the
same order as the Fast Fourier Transform, which is also O(N log N)
and produces an N-point result. Since the GFT formalism describes
the STFT, GST and DOST, the fast GFT algorithm provides a
computationally and memory efficient method for calculating all
three transforms.
[0088] The fast GFT algorithm presented here is the minimally
sampled, one-to-one transform. This configuration is efficient to
compute and store, but some applications may benefit from a less
optimal implementation. If the frequency domain windows are allowed
to overlap, introducing some redundancy, a smoother window shape
may be achieved. This redundancy comes at the cost of some
computation time and additional storage requirements. For example,
if each window is allowed to overlap its neighbors by 50%, each FFT
in Algorithm 1 (step 5) will require transforming twice as many
points, increasing the computation time and producing an S-domain
with 2N points. However, for constant overlap factors, such as in
this example, the computational complexity of the algorithm will
remain at O(N log N).
[0089] Actual computation times for the frequency domain DST, the
fast GFT (Algorithm 1), and an optimized FFT library are compared
in FIG. 7. Due to memory constraints, the largest signal that could
be calculated with the DST algorithm had N=2.sup.13 points. While
the DST of this signal took over 30 s, the fast GFT and the FFT
each took only a few milliseconds: 5 ms and 2 ms, respectively.
[0090] The efficiency of the fast GFT allows computation of the
transform for larger signals and images. For example, FIG. 8A shows
the first 8 s of one embodiment of an electrocardiogram (ECG)
recording. In this case, ST-T refers to ST- and T-waves, which are
features of an ECG and are not related to the S-transform. This
sample is from the V4 lead of record e0103, a 62-year-old male
subject complaining of mixed angina. The complete recording
contains over one million points. FIG. 8B is a portion of the
spectrum produced by a transform of the first 2.sup.20 points from
this signal, using the fast GFT with adaptive Gaussian windows. In
this configuration the GFT is a fast ST. The transform took 0.77 s
to calculate. The classic frequency domain S-transform algorithm,
if it could be calculated entirely in memory, would take over three
days. However, since the N.times.N spectrum requires 16 TB of
storage space (64-bit complex data), swapping to an array of hard
drives would greatly increase the computation time beyond even this
figure. FIG. 8C is the same spectrum with bilinear
interpolation.
The Inverse Fast General S-Transform
[0091] The forward continuous GFT is calculated via the
.alpha.-domain by shifting and windowing the Fourier transform of
the signal, according to Eq. (6), then inverse Fourier transforming
each horizontal .nu.'=line to produce the S-domain, according to
Eq. (7). The inverse continuous GFT reverses these steps, using the
FT of each horizontal .tau.-line in the S-domain to recover the
.alpha.-domain, according to Eq. (5). As discussed in Section
above, the signal can then be obtained from any horizontal line in
the .alpha.-domain by dividing by W(.nu.',.sigma.), shifting, and
taking the inverse FT. Alternately, the complete frequency domain
can be reconstructed directly from the S-domain by summing over
.tau..
[0092] The fast GFT avoids redundancy by using only the central
band-passed portion of the .alpha.-domain along the .nu.'-axis and
sampling on a dyadic scale along the .nu.-axis. This scheme samples
along horizontal line segments in the .alpha.-domain (FIG. 5C). The
inverse FFT of each line segment produces a horizontal line of the
S-domain (FIG. 5D). Since the .alpha.-domain is sampled using a
total of N-points, the corresponding S-domain will also be sampled
with a total of N-points.
[0093] Since the S-domain is sampled using only N points, the
inverse fast GFT may require all N points to reconstruct the
signal. The FFT of each horizontal line of the S-domain may recover
a line segment in the .alpha.-domain. When each segment is divided
by W(.nu.',.sigma.) it may represent an unweighted portion of the
frequency domain. As illustrated in FIG. 6, when shifted by .nu.,
these segments combine to produce a full frequency domain, which
can be inverse fast Fourier transformed to recover the signal.
Algorithm 2 calculates this inverse fast GFT in-place, making it
optimally memory efficient.
[0094] The following pseudo code describes one embodiment of a
process for computing the inverse in-place fast general Fourier
family transform, using a dyadic sampling scheme for this
example:
Algorithm 2
TABLE-US-00002 [0095] LET A.sub..tau.[0:N]be an array of N samples
of the S-domain as produced by FOR f start in 0 , 1 , 2 , 4 , , N 8
, N 4 DO : ##EQU00018## IF f.sub.start < 2: LET f.sub.width = 1,
f.sub.centre = f.sub.start + 1 ELSE: LET f width = f start , f
centre = f start + f width 2 + 1 ##EQU00019## LET f.sub.end =
f.sub.start + f.sub.width LET f.sub.start' = -f.sub.start,
f.sub.end' = -f.sub.end and f.sub.centre' = -f.sub.centre 1.
Calculate the window function and FFT: W[0:N] = FFT.sub.t.fwdarw.v
(w[0:N]) 2. Calculate the in-place FFT of the S-domain band:
A.sub.v[f.sub.start:f.sub.end] =
FFT.sub..tau..fwdarw.v(A.sub..tau.[f.sub.start:f.sub.end]) 3. Shift
the window function to be centered on the frequency band: W[0:N] =
shift (W[0:N], f.sub.centre) 4. Remove the window: A v [ f start :
f end ] = A v [ f start : f end ] W [ f start : f end ]
##EQU00020## 5. Repeat steps 2 to 4 substituting f.sub.start',
f.sub.end' and f.sub.centre' to calculate the negative frequencies.
END FOR 6. Recover the signal through the inverse FFT: A.sub.t[0:N]
= FFT.sup.-1.sub.v.fwdarw.t(A.sub.v[0:N])
Time-Domain S-Transform
General Transform
[0096] Though the different transforms, particularly the Fourier
and wavelet transforms, are often considered to be distinct
entities, they have many similarities. To aid comparison of the ST
with other transforms and help translate techniques developed for
one transform to be used with another, the inventors present
several common transforms in a unified context. Previous
investigators have noted the similarities between the FT, STFT and
wavelet transform and the utility of representing them in a common
context. To this end, generalized transforms that describe all
three have been constructed. However, to our knowledge, previous
generalized formalisms do not explicitly specify separate kernel
and window functions. Separating the two may better illustrate the
relationships between the transforms, particularly when the ST is
included. Eq. 2 describes a generalized integral transform of a
signal, s(t), multiplied by a kernel function, .phi.(t,.tau.,.nu.),
and a window, w(t,.tau.,.nu.).
T ( .tau. , v ) = .intg. - .infin. + .infin. s ( t ) w ( t , .tau.
, v ) .phi. ( t , .tau. , v ) t ( 9 ) ##EQU00021##
[0097] The result is a spectrum, T(.tau.,.nu.), which has both
time, .tau., and frequency (or analogue), .nu., axes. The kernel
determines the target domain of the transform and provides
localization in that domain (the .nu. axis). The window provides
localization in the original domain (the .tau. axis). The variable
.tau. determines the centre of the window and is equivalent to the
shift parameter of the wavelet transform.
[0098] In many common cases w and .phi. are symmetric, so the
correlation in Eq. (9) is also a convolution, making this product,
w.phi., the impulse response function or kernel of the system. To
remain as consistent as possible with the literature, the inventors
will use the term kernel to refer to .phi. alone, while the product
of w and .phi. will be referred to as a basis function.
[0099] Six transforms are presented for comparison in Table 1 using
the general transform formalism of Eq. (9). These can be (and
usually are) divided into two families by considering whether the
basis functions are shifted and scaled versions of each other or
not. The Fourier family, including the FT, STFT, ST and sinc
transforms (which is discussed later), use the familiar complex
sinusoidal kernel e.sup.i2.pi..nu.t, which does not depend on
.tau., and produce a globally referenced phase spectrum. The
wavelets, including the Morlet and Haar wavelets, use kernels that
are shifted and scaled versions of a mother wavelet. In the general
transform formalism this shifting is characterized as a dependence
of the kernel on .tau., and accounts for the locally referenced
phase spectrum of complex wavelets like the Morlet.
[0100] In the context of the general transform, the four common
members of the Fourier family considered here differ only in their
window functions. The Fourier transform is the simplest and uses a
constant for the window. This provides no time localization and
maximal frequency resolution (the Fourier transform of a constant
is a delta-function), transforming all the way from the time to the
frequency domain. The Fourier kernel is not shifted across the
signal (it does not depend on .tau.), ensuring that all phase
measurements are globally referenced. Since neither the kernel nor
the window depend on .tau., the variable .tau. can be eliminated
from the equation, leaving the familiar definition of the Fourier
transform. The STFT can be understood as a transform using a
Fourier kernel and a window, often a boxcar or a Gaussian, that
does not depend on .nu., i.e. the window is not frequency adapted.
The window provides a constant time (and thus frequency)
resolution.
[0101] In terms of the general transform, the S-transform is very
similar to the STFT. It also uses the Fourier kernel and a Gaussian
window, with one important difference: in this case the window
function does depend on .nu. and so the window is adapted depending
on the frequency under examination, providing coarse frequency (and
finer time) resolution at higher frequencies, much as the wavelet
transform does.
[0102] The wavelet transform family does not use the Fourier
kernel, but rather shifted and scaled versions of a mother wavelet.
Two common transforms are the Morlet and Haar wavelets (Table 1).
The Morlet mother wavelet consists of a complex sinusoid multiplied
by a Gaussian, for time localization. This mother wavelet is then
scaled and shifted across the signal during calculation of the
transform. In the context of the general transform, the inventors
can split this mother wavelet into a complex sinusoidal kernel and
a Gaussian window multiplied by a boxcar to make it finite, both of
which are similar (though not identical) to those of the
S-transform. In one embodiment, the ST may be derived as a phase
correction to the Morlet wavelet. Substitution of these functions
into the general transform produces an equation that is equivalent
to the definition of the Morlet wavelet transform, differing only
in the convention of having the kernel function depend on .sigma.
and .tau. directly, instead of defining a mother wavelet and
subsequently shifting and scaling it.
[0103] Finally, the Haar wavelet utilizes a kernel that appears to
depart entirely from the Fourier kernel. The Haar mother wavelet,
even when divided into window and kernel functions, is non-periodic
and discontinuous, unlike the other transforms considered
previously. For consistency with the other transforms, this
difference can be accommodated by convolving the defined Haar
mother wavelet with a comb function (making it periodic) then using
a boxcar function as the window to zero all but the desired
period.
The Fast S-Transform
[0104] Though calculating a discrete approximation of a continuous
transform is useful, as with the continuous wavelet and
S-transforms, a fully discrete approach makes optimal use of
knowledge of the sampling process applied to the signal to decrease
the computational and memory resources required. In this section a
discrete fast S-transform (FST) is developed by utilizing
properties that apply to all of the discrete versions of transforms
described by the general transform of Eq. (9).
[0105] A sampled signal has two important features--the sampling
period, .DELTA.t, and the number of samples, N. Multiplying these
two values gives the total signal length, W.sub.t. Sampling the
signal and limiting it to finite length imposes several limitations
on the transformed spectrum. The Fourier transform is the simplest
case. The discrete Fourier transform (DFT) of a signal is limited
to the same number of samples, N, as the original signal,
conserving the information content of the signal. The highest
frequency that can be represented, .nu..sub.max, is the Nyquist
frequency, which is half the signal sampling frequency,
1 2 .DELTA. t . ##EQU00022##
The sampling period of the frequency spectrum, .DELTA..nu., is the
reciprocal of the signal length, 1/W.sub.t.
[0106] The result of the continuous S-transform of a
one-dimensional signal is S(.tau.,.nu.), a spectrum with both time
and frequency axes. A fast ST may sample the .tau.-.nu. plane
sufficiently such that the transform can be inverted (i.e., without
loss of information) but also avoid unnecessary oversampling. If
the original signal contains N points, the inventors possess N
independent pieces of information. Since information cannot be
created by a transform (and may be conserved by an invertible
transform) an efficient ST will produce an N point spectrum, as
does the DFT. The original definition of the ST produces a spectrum
with N.sup.2 points. Therefore, for all discrete signals where
N>1, the continuous ST is oversampled by a factor of N.
[0107] In addition, the ST varies temporal and frequency resolution
with changing frequency: higher frequency/lower time resolution at
low frequencies and the converse at higher frequencies, but this is
not reflected in the uniform N.sup.2-point .tau.-.nu. plane
sampling scheme. According to the uncertainty principle, at higher
frequencies the FST should produce .tau.-.nu. samples that have a
greater extent along the frequency axis and lesser extent along the
time axis, as is seen in the discrete wavelet transform. The cost
of this oversampling is evident in the increased computational
complexity and memory requirements of the ST.
[0108] The dyadic sampling scheme used by discrete wavelet
transforms provides a progressive sampling scheme that matches
underlying resolution changes. In light of the similarities between
the WT and ST illustrated by the general transform formalism, a
dyadic sampling scheme can be used to construct a discrete ST. In
the case of the ST of a complex signal, a double dyadic scheme may
be used to cover both the positive and negative frequency spectra.
In this arrangement, time resolution is at a minimum for low
frequencies, both positive and negative, and increases with the
absolute value of frequency. For comparison, the .tau.-.sigma.
plane sampling schemes for the signal, Fourier transform,
continuous S-transform and discrete wavelet transform are shown in
FIG. 9. In FIG. 9A is a uniformly sampled signal with N=8, FIG. 9B
is the Fourier transform, FIG. 9C is the S-transform, and FIG. 9D
is a discrete Wavelet transform. FIG. 9A provides no frequency
localization, FIG. 9B provides no time localization and both FIG. 9
C and FIG. 9D provide time and frequency information
simultaneously. While FIGS. 9A, B, and D each have N samples, FIG.
9C has N.sup.2. In FIG. 9D, because .nu. corresponds to wavelet
scale and not frequency, it is positive. The double dyadic scheme
of the FST is illustrated in FIG. 10, and a particular algorithm
for calculating the FST with this type of sampling scheme is
presented in the pseudo code for computational Algorithm 3
described below. No N-sample octave arrangement with integral
sample locations exists for the complex FST of a power of two
signal.
Algorithm 3: The Discrete S-Transform with 4.times.
Oversampling
TABLE-US-00003 1. Calculate the Fourier transform of the signal , H
( n NT ) ##EQU00023## 2. Pre-calculate the required kernel
functions: .phi. + = ( kT , n NT ) = e - i 2 .pi. kn N and .phi. -
= ( kT , n NT ) = e i 2 .pi. kn N ##EQU00024## 3. Pre - calculate
the window functions : w ( kT , n NT ) ##EQU00025## FOR n in { N 2
, N 4 , N 8 , , 4 , 2 , 1 } DO : ##EQU00026## 4. Band pass filter
H(.kappa.): H ' ( .kappa. ) = H ( .kappa. ) where n 2 < .kappa.
.ltoreq. n ##EQU00027## and inverse FT to obtain h'(t). FOR every
point j in h'(t) DO: 5. Calculate the transform samples: S ( jT , 3
n 4 NT ) = k = 0 N - 1 h ' ( kT ) w ( kT - T , 3 n 4 NT ) .phi. + (
kT , 3 n 4 NT ) ##EQU00028## S ( jT , - 3 n 4 NT ) = k = 0 N - 1 h
' ( kT ) w ( kT - T , 3 n 4 NT ) .phi. - ( kT , 3 n 4 NT )
##EQU00029## END FOR END FOR
[0109] The sampling in FIG. 10A is very similar to the discrete
wavelet sampling scheme but adds negative frequencies, resulting in
oversampling by a factor of approximately 2. However, samples are
distributed consistently. FIG. 10B depicts an additional 2.times.
oversampling of the .tau.-axis. This allows the representation of
odd to even point transitions in the signal in the appropriate
frequency band, and also allows the use of windows that are not
perfect Gaussians. Since the oversampling in both FIG. 10A and FIG.
10B is a constant factor, it does not affect the computational
order of the FST.
[0110] The result from Algorithm 3 is presented for a sample signal
in FIG. 11, along with the corresponding spectrum from the
continuous ST. FIG. 11A consisting of a low frequency for the first
half of the signal, followed by a higher frequency for the second
half A high frequency chirp is present in the middle of the first
half of the signal and a delta in the second. FIG. 11B is the
magnitude of the positive frequency half of the continuous
S-transform of FIG. 11A and FIG. 11C is the magnitude of the
positive frequency half of the discrete S-transform of FIG.
11A.
[0111] In Algorithm 3, the highest two frequencies (positive and
negative) are calculated at the full resolution of the signal.
Using dyadic sampling, the next frequencies to be sampled are half
the magnitude of the first two. Since the frequency may be lower,
the signal may be downsampled (.DELTA.t can be increased), if the
Nyquist criterion is satisfied. This is accomplished in steps 5 and
6 of Algorithm 3. Since Algorithm 3 uses four points per wavelength
instead of two, evaluating the transform with the window centered
at every point results in a factor of two oversampling of the
.tau.-axis (FIG. 10B) compared to the discrete wavelet
transform.
[0112] The Gaussian windows of the ST are infinite in extent but
when calculating the transform of a finite length signal, the
Gaussians are necessarily multiplied by a boxcar window the width
of the signal itself. Therefore, the actual window for any finite
length signal is a Gaussian multiplied by a boxcar. This situation
is particularly apparent at lower frequencies, where the Gaussians
are wider (FIG. 12). FIGS. 12A-D illustrate embodiments of
S-transform kernels (gray) and windows at .tau.=0.5 (black) for
.nu.=1 (FIG. 12A), .nu.=2 (FIG. 12B), .nu.=4 (FIG. 12C) and .nu.=8
(FIG. 12D). Dotted lines indicate the imaginary component of the
kernel while solid lines are the real part. Black circles
superimposed on the real kernel indicate the sampling positions for
the FST. As the frequency increases the number of sampling points
required also increases. In the continuous ST, the Gaussian window
is scaled with frequency but the boxcar is not. This contrasts with
the Morlet wavelet, which, using the general formalism of Eq. (9),
can also be defined with a window that is a composite of a Gaussian
and a boxcar. However, in the Morlet wavelet the Gaussian and
boxcar are scaled together. Scaling of both parts of the window
function is a key refinement, as it both produces more consistent
windows and significantly decreases the computational complexity of
the FST.
[0113] According to one embodiment, the signal processor 104 may
only compute the sums in step 4 of Algorithm 3 for non-zero points
(those that are inside the boxcar window). The boxcar may be wide
enough to contain at least one entire wavelength, but the width may
not be a multiple of the wavelength. This effectively decreases
W.sub.t, with the full signal length required only for calculating
the DC component, and shorter portions used at higher frequencies.
Since the inventors are downsampling in step 6, .DELTA.t is
smallest at high frequencies and becomes progressively larger at
lower frequencies and so the summation operation in step 4 may be
over a constant number of points. The examples according to the
present embodiments may use 8 points. This reduces the complexity
of step 4, which is nested inside two FOR loops, from O(N) to
constant complexity: O(C). As an additional benefit, adjusting the
width of the boxcar window greatly reduces the number of kernels
and windows that may be pre-calculated in steps 2 and 3 of
Algorithm 3, since the kernels and windows remain essentially
constant while the signal's length and resolution are
manipulated.
[0114] Algorithm 3 produces a discrete fast S-transform that
samples the .tau.-.nu. plane in an octave pattern, with the number
of .tau. samples increasing with frequency while the number of .nu.
samples decreases. This sampling scheme produces 2N samples,
representing a two times oversampling, compared to the theoretical
minimum number of samples. This oversampling, as compared to the
DWT, may be used to achieve a globally referenced phase
measurement. The phase spectrum requires the use of complex
sinusoid basis functions, which may be represented by at least
three points (in this case, the inventors use four), and, for
complex signals, necessitates the calculation of the negative
frequency spectrum, in contrast to the wavelet. However, this
oversampling is a constant factor and does not affect the
computational complexity of the transform.
[0115] The constant term in step 3 of Algorithm 3 is used to widen
the Gaussian window slightly. This term is also used in the
continuous ST. For c>1 its effect is to smooth out the
S-spectrum along the .tau.-axis, at the expense of temporal
resolution. This produces a more intuitive spectrum appearance: for
example, a sinusoidal signal appears as a uniform amplitude line in
the S-spectrum rather than a line with slightly varying amplitude.
As long as c.gtoreq.1, its value does not affect the invertibility
of the transform, and c=1.5 in one example.
[0116] The computational complexity of the FST, calculated in the
appendix, is O(N log N)--the same as that of the Fourier transform.
The storage requirements are O(N), like the FFT and discrete
wavelet transforms. Since the FST uses the same Gaussian windows as
the ST and these windows are infinite in extent both in the time
and frequency domains, the FST basis functions are not orthogonal.
However, the Gaussian windows in both domains decay quickly so the
spectrum is a reasonable representation, both in frequency and
time, of a true localized spectrum, as it is with the continuous
ST. In addition, with a computational complexity equal to that of
the Fourier transform, the redundant information in the FST is not
overly burdensome.
The Inverse Fast S-Transform
[0117] A discrete version of the S-transform may be invertible by
the same procedure as the inverse continuous ST: summation of the
transform space over the .tau.-axis, producing the Fourier
spectrum, which can then be inverse discrete or fast Fourier
transformed to obtain the original signal. The inverse discrete and
fast Fourier transforms require a coefficient for each integer
value of the frequency index variable .nu. from
- v max 2 ##EQU00030##
to
+ v max 2 . ##EQU00031##
Since the FST uses an octave system along the .nu.-axis, the
missing coefficients may be calculated from known points and a
priori information.
[0118] Multiplying a signal by a window is equivalent to convolving
the Fourier transform of the signal with the Fourier transform of
the window. Thus, for a given .nu.-line in the FST, by examining
the Fourier transform of w(.tau.,t,.parallel..nu.|).phi.(t,.nu.)
the inventors can examine the extent and position of the windows in
the .nu.-domain. FIG. 13 shows this for .nu.=4. Recall from
Algorithm 3 that calculating the .nu.=4 line involves positioning
the window at eight locations along .tau., and the signal itself
has been downsampled to sixteen points (four per wavelength). The
effect of moving the window across the signal is to apply a phase
ramp to the Fourier transform of the window. In the ST, and for
calculated values of .nu. in the FST, when summed, the .nu.-space
windows interfere and collapse to a delta at a frequency of
.nu..
[0119] For example, FIG. 13 illustrates the discrete S-transform
window-kernel product in the .nu.-domain for .nu.=4 at .tau.=0 (A),
.tau.=1 (B), .tau.=2 (C), .tau.=3 (D), .tau.=4 (E), .tau.=5 (F),
.tau.=6 (G) and .tau.=7 (H). In this figure dotted lines may
indicate the imaginary part while solid lines may indicate the real
component and circles may be located at the sampling points.
Shifting the window along t may vary the phase in the .nu.-domain.
If each is weighted by the factor e.sup.-i2.pi.l.tau., the sum may
be a delta shifted by l and multiplied by the magnitude of (A) at
.nu.+l.
[0120] To recover a complete Fourier spectrum from the FST, the
.nu.-domain windows may be weighted so they still sum to a delta,
but one that is displaced from .nu.. This can be done by using a
complex sinusoid as the weighting function. Multiplying a .nu.-line
of the FST by a weighting sinusoid with frequency l, then summing,
produces the Fourier coefficient for frequency .nu.+l. The
recovered Fourier coefficients may be windowed by the Fourier
transform of the time domain window (another Gaussian) meaning the
magnitude of the recovered Fourier coefficients is attenuated by an
amount A.sub..nu.l, equal to the magnitude of the Gaussian window
in the .nu.-domain at frequency .nu.+l. Dividing by this value
provides the correct Fourier coefficient. Note when l=0,
A.sub..nu.l=1 since the Gaussian windows are normalized.
[0121] Certain limits may be observed in this inversion process.
Each .nu.-line has a number of samples that is less than or equal
to the number of samples in the original signal, with fewer samples
at lower frequencies. The sinusoidal weighting function may satisfy
the Nyquist criterion so, for example, the line at .nu.=4, with
eight samples, can be used to calculate only Fourier coefficients
.nu.={0, 1, 2, . . . , 7}. For eight samples, the Nyquist frequency
is 4 so for these frequencies -3.ltoreq.l.ltoreq.3 and the sampling
rate of the .nu.-line is sufficient. In practice, the range is
somewhat less as the .nu.-space windows decay quickly for higher
values of l and may even possess zeros since the signal-space
windows are a Gaussian multiplied by a boxcar, which is a Gaussian
convolved with a sinc in the .nu.-domain. From a noise perspective,
if two .nu.-lines both allow the calculation of a particular
Fourier coefficient, it is best to use the one that has been
attenuated least, i.e., the one with the largest value of
A.sub..nu.l.
[0122] This inversion process is equivalent to performing a Fourier
transform on each .nu.-line in the FST, so it can be accomplished
efficiently by fast Fourier transforming each line and weighting
the resulting coefficients appropriately. Where there is overlap,
the minimum noise coefficients can then be chosen to make up the
final Fourier spectrum.
[0123] FIG. 14 shows the process of transform and inversion for a
sample biomedical signal, using both the continuous ST and the FST.
The signal (FIG. 14A) is a short electrocardiogram (ECG) recording
from a publicly available subset of the European ST-T Database,
consisting of the first 512 samples from the V4 lead of record
e0103, a 62 year old male subject complaining of mixed angina. This
signal was transformed using both the ST and FST, then inverted.
The left column shows the results of the ST while the right is from
the FST. In FIG. 14B are the magnitudes of the positive frequency
halves of the transforms, FIG. 14C shows the magnitude of the
Fourier spectra resulting from inversion and FIG. 14D is the
signals recovered with the inverse FT. Note that in FIG. 14B
frequency is on the vertical axis while in FIG. 14C it is on the
horizontal axis. The results of the inversion of both transforms
are identical to each other and to the original signal, within
machine precision. Though the continuous ST produces a smoother
looking spectrum, this difference is illusory--since the same
information is contained in both the ST and FST spectra, the
smoothness is produced by an interpolation process inherent in the
continuous ST oversampling. This is demonstrated by the fact that
the original function can be recovered by inverting either
dataset.
The Discrete Orthonormal S-Transform
[0124] In one embodiment, the basis function of the DOST may be a
sum of complex sinusoids multiplied by "frequency shift" and "phase
shift" terms. In one embodiment, the frequency shift term can be
understood as the kernel and the phase shift term as the .tau.
dependence of the window function. The sum of sinusoids term is
equivalent to a boxcar function in the discrete frequency domain.
Taking the Fourier transform of this results in an aliased sinc
function, i.e. one that contains only frequencies corresponding to
the discrete Fourier coefficients. This sinc is the window function
for the DOST. In a continuous formulation the boxcar function in
the frequency domain is a continuous boxcar and so the window in
the time domain is a true sinc function. Though the kernel function
of the DOST is the same as that of the ST and FST, the sinc window
is a significant departure. It may be better to consider the DOST
not as a discrete formulation of the ST, which is defined as using
a Gaussian window, but rather as a discrete form of a new
transform, the continuous sinc transform (CST), named for the
window function.
[0125] The use of a sinc window has several consequences. The sinc
is infinite in extent but its Fourier pair, the boxcar function, is
finite, so DOST basis functions can be constructed which are
orthonormal. Interestingly, Stockwell suggests that some redundancy
may be desirable and identifies the possibility for non-orthogonal
basis functions as a strength of the DOST. Though the DOST window
is a very desirable boxcar function in the frequency domain, it is
a much less desirable sinc function in the time domain. The sinc
has considerably larger side lobes than the Gaussian, and since
each sample in the DOST spectrum represents a local region defined
by the window, any filtering operations performed on the spectrum
will exhibit classic Gibb's ringing in the time-domain. Similarly,
localization on the .tau.-axis of the transform spectrum will
suffer, again because of the large sidelobes of the sinc. These may
prove to be serious drawbacks in medical applications such as
filtering and techniques that examine restricted areas of
pathology, such as brain tumors.
[0126] The algorithm proposed by Stockwell for calculating the DOST
uses the full signal length and resolution for every calculation so
its computational complexity is O(N.sup.2) compared to O(N log N)
for the FST. The storage requirements of both the DOST and FST are
the same, O(N).
[0127] The examples of the FST and DOST suggest the construction of
other transforms using the general transform framework. Many
functions have been suggested for use as windows in Fourier
analysis, and the properties that make these windows desirable in
that application are also advantageous in the construction of new
time-frequency transforms. For example, a FST analogue can be
constructed that replaces the Gaussian window with a
Blackmann-Harris window or Hamming window. Windows have been
constructed that are specifically designed to be finite but
minimize the ringing associated with the boxcar-sinc pair, which
allows an orthogonal DOST-like transform with fewer time domain
ringing effects. Prolate spheroidal wavefunctions, which can be
approximated by Kaiser-Bessel functions, have the smallest
simultaneous regions of support in both domains and so may provide
maximum simultaneous localization in both the frequency and time
domains.
GFT for a Multi-Dimensional Signal
GFT for a 2D Signal
[0128] When the signal is two-dimensional, its resulting GFT is
four-dimensional. Let the signal be written as s(x.sub.1, x.sub.2)
and its GFT as .PSI.(.chi..sub.1,.chi..sub.2,k.sub.1,k.sub.2), with
.chi..sub.m the analogue of x.sub.m, and k.sub.m the Fourier
conjugate of x.sub.m for m={1,2}, such that
.PSI. ( .chi. 1 , .chi. 2 , k 1 , k 2 ) = .intg. - .infin. +
.infin. .intg. - .infin. + .infin. s ( x 1 , x 2 ) w ( .chi. 1 - x
1 , .chi. 2 - x 2 , k 1 , k 2 ) - 2 .pi. ( k 1 x 1 + k 2 x 2 ) x 1
x 2 ( 10 ) ##EQU00032##
[0129] It may be assumed that the general window function w(*) is
separable, and by extension orthogonal, in
(.chi..sub.m-x.sub.m,k.sub.m), i.e.,
w(.chi..sub.1-x.sub.1,.chi..sub.2-x.sub.2,k.sub.1,k.sub.2).fwdarw.w.sub.-
1(.chi..sub.1-x.sub.1,k.sub.1)w.sub.2(.chi..sub.2-x.sub.2,k.sub.2)
(11)
[0130] Note that this is actually not overly restrictive; rather,
it ensures that the multi-dimensional GFT be separable along its
Cartesian axes, just like the multi-dimensional FT. Analogous to
its 1D counterpart, a-space is now synthesized via the 2D-FT of
.PSI.(*) along the .chi..sub.m-axes. Then, by (i) letting k.sub.m '
be the Fourier conjugate of .chi..sub.m, (ii) substituting Eq. (11)
into Eq. (, and (iii) re-ordering terms, we obtain
.alpha. ( k 1 ' , k 2 ' , k 1 , k 2 ) = .intg. - .infin. + .infin.
.intg. - .infin. + .infin. .PSI. ( .chi. 1 , .chi. 2 , k 1 , k 2 )
- 2 .pi. ( k 1 ' .chi. 1 + k 2 ' .chi. 2 ) .chi. 1 .chi. 2 = .intg.
- .infin. + .infin. .intg. - .infin. + .infin. w 2 ( .chi. 2 - x 2
, k 2 ) - 2 .pi. ( k 2 x 2 + k 2 ' .chi. 2 ) [ .intg. - .infin. +
.infin. .intg. - .infin. + .infin. s ( x 1 , x 2 ) w 1 ( .chi. 1 -
x 1 , k 1 ) - 2 .pi. ( k 1 x 1 + k 1 ' .chi. 1 ) x 1 .chi. 1 ] x 2
.chi. 2 ( 12 ) ##EQU00033##
[0131] In one embodiment, the square-bracketed term may be
expressed, such that
.alpha. ( k 1 ' , k 2 ' , k 1 , k 2 ) = .intg. - .infin. + .infin.
.intg. - .infin. + .infin. w 2 ( .chi. 2 - x 2 , k 2 ) - 2 .pi. ( k
2 x 2 + k 2 ' .chi. 2 ) [ W 1 ( k 1 ' , k 1 ) S 1 ( k 1 ' + k 1 , x
2 ) ] x 2 .chi. 2 = W 1 ( k 1 ' , k 1 ) .intg. - .infin. + .infin.
.intg. - .infin. + .infin. S 1 ( k 1 ' + k 1 , x 2 ) w 2 ( .chi. 2
- x 2 , k 2 ) - 2 .pi. ( k 2 x 2 + k 2 ' .chi. 2 ) x 2 .chi. 2 ( 13
) ##EQU00034##
where W.sub.1(k.sub.1,k.sub.1) is the FT of
w.sub.1(.chi..sub.1,k.sub.1) along the .chi..sub.1-axis, and
S.sub.1(k.sub.1,x.sub.2) is the FT of s(x.sub.1,x.sub.2) along the
x.sub.1-axis. Furthermore, Eq. (13) may also be expressed such
that
.alpha.(k.sub.1',k.sub.2',k.sub.1,k.sub.2)=W.sub.1(k.sub.1',k.sub.1)[W.s-
ub.2(k.sub.2',k.sub.2)S.sub.2(k.sub.1'+k.sub.1,k.sub.2'+k.sub.2)]
(14)
where W.sub.2(k.sub.2',k.sub.2) is the FT of
w.sub.2(.chi..sub.2,k.sub.2) along the .chi..sub.2-axis, while
S.sub.2(k.sub.1'+k.sub.1,k.sub.2) is the FT of
S.sub.1(k.sub.1'+k.sub.1,x.sub.2) along the x.sub.2-axis. Of
course, S.sub.2(*).fwdarw.S(*) since S.sub.2(*) can effectively be
recognized as the shifted 2D-FT of s(x.sub.1,x.sub.2), from which
we obtain the compact expression described in Eq. (15).
.alpha.(k.sub.1',k.sub.2',k.sub.1,k.sub.2)=W.sub.1(k.sub.1',k.sub.1)W.su-
b.2(k.sub.2',k.sub.2)S(k.sub.1'+k.sub.1,k.sub.2'+k.sub.2) (15)
GFT for an M-D Signal
[0132] When the signal is M-dimensional, say s(x.sub.1, . . . ,
x.sub.M), its GFT is 2M-dimensional, namely .PSI.(.chi..sub.1, . .
. , .chi..sub.M,k.sub.1, . . . , k.sub.M), where .chi..sub.m and
k.sub.m are still the analogue and Fourier conjugate of x.sub.m,
respectively, but with m={1, 2, . . . , M}. In one embodiment, the
M-D window function w(*) may imply that
w(.chi..sub.1-x.sub.1, . . . , .chi..sub.M-x.sub.M,k.sub.1, . . . ,
k.sub.M).fwdarw.w.sub.1(.chi..sub.1-x.sub.1,k.sub.1) . . .
w.sub.M(.chi..sub.M-x.sub.M,k.sub.M) (16)
[0133] As demonstrated in the previous subsection, a-space for a 2D
signal is a simple extension of its 1D counterpart: effectively the
terms are re-arranged to form nested loops of a-space integrals of
effective 1D signals. This re-ordering and nesting may be possible
because (i) the Fourier transform is both linear and orthogonal,
and (ii) the window functions are separable along their orthogonal
Cartesian axes. Extending this to M-dimensions, then, it may be
shown that the general result
.alpha.(k.sub.1', . . . , k.sub.M',k.sub.1, . . . ,
k.sub.M)=W.sub.1(k.sub.1',k.sub.1) . . .
W.sub.M(k.sub.M',k.sub.M)S(k.sub.1'+k.sub.1, . . . ,
k.sub.M'+k.sub.M) (17)
where k.sub.m ' is the Fourier conjugate of .chi..sub.m,
W.sub.m(k.sub.m',k.sub.m) is the FT of w.sub.m(.chi..sub.m,k.sub.m)
along the .chi..sub.m-axis, and S(k.sub.1, . . . , k.sub.M) is the
M-dimensional FT of s(x.sub.1, . . . , x.sub.M).sub..
The Fast GFT Algorithm
[0134] One inefficiency of the conventional S-transform (and GFT)
algorithm is that it is not a one-to-one transform nor does it
provide multi-resolution or multi-scale information, in
contradistinction to the WT algorithm. This is most apparent from
visualizing .alpha.-space for a 1D signal whereby each k-line is a
shifted version of S(k') multiplied by W(k',k). In such an
embodiment, the width of W(k',k) is small at low frequencies, and
increases as |k| increases. However, in comparing successive
k-lines of a-space, the relative shift of S(k') is just one
k'-pixel, and the W(k',k) widths are similar (but not equal).
[0135] In one embodiment, W(k',k) may effectively provide a
progressive bandpass width over a "block" of S(k'). .PSI.(.chi.,k)
may be calculated via the IFT for blocks whose length is less than
that of S(k'). In an embodiment having dyadic sampling in |k|, the
block sizes are themselves dyadic, such that the total number of
points remains N. In one embodiment, the k'-resolution in each
block may remain constant. Therefore, upon IFT, each .PSI.(.chi.,k)
.chi.-line at a specific k-value has the same extent (i.e.,
field-of-view) as the original signal, but with different
resolution.
[0136] For M-dimensional signals, the blocks will of course be
M-dimensional, and the dyadic scale will occur in all directions.
This "blocking" is easily visualized for M=2, where the centre
(k'.sub.1,k'.sub.2)-value of each block may also be its
(k.sub.1,k.sub.2)-value. In one embodiment, the effective 2D
weighting function may be
W.sub.1(k'.sub.1,k.sub.1)W.sub.2(k'.sub.2,k.sub.2). Extending this
beyond 2D may be possible. In one embodiment, the efficient FGFT
algorithm involves 4 steps: [0137] 1. Take the M-dimensional FT of
signal s(x.sub.1, . . . , x.sub.M) to produce S(k.sub.1, . . . ,
k.sub.M) [0138] 2. Synthesize the M-dimensional, "blocked"
weighting function [0139] 3. Multiply S(k.sub.1, . . . , k.sub.M)
from Step 1 with the weighting function of Step 2 to get .alpha.(*)
[0140] 4. Perform the in-place, M-dimensional IFTs for each block
of .alpha.(*) in Step 3 to generate the desired transform,
.PSI.(*)
Computational Burden
[0141] The computational order for the multi-resolution, one-to-one
GFT of an M-dimensional signal of size
N.sub.1.times.N.sub.2.times., . . . , .times.N.sub.M may be
calculated by investigating the order for each of the four steps
required in the algorithm. Step 1 is the M-dimensional FT of the
signal, whose order may be known to be N.sub.1N.sub.2 . . . N.sub.M
log(N.sub.1N.sub.2 . . . N.sub.M). Step 2 can be pre-calculated,
and involves the usual memory allocation and related overhead. Step
3 is a point-by-point multiplication of complex function S(*) by
the real weighting function, such that its order is 2N.sub.1N.sub.2
. . . N.sub.M. Step 4 involves block-wise IFTs, which can be
calculated in-place. To calculate the effective number of
operations required, let the number of points in the jth M-block be
P.sub.1,jP.sub.2,j . . . P.sub.M,j. Since we have a one-to-one
transform, then
j = 1 B m P m , j = N m ##EQU00035##
for m={1, 2, . . . , M}, where B.sub.m is the number of blocks
along the mth dimension. And, let P.sub.m,j=f.sub.m,jN.sub.m with
f.sub.m,j.epsilon.(0,1) for all m and j, where f.sub.m,j represents
the fraction of points along dimension m for block j. Thus, we
have
j = 1 B m P m , j = N m j = 1 B m f m , j N m = N m j = 1 B m f m ,
j = 1 , .A-inverted. m ( 18 ) ##EQU00036##
[0142] Since the IFT of each P.sub.1,j . . . P.sub.M,j-block is of
computational order P.sub.1,j . . . P.sub.M,j log(P.sub.1,j . . .
P.sub.M,j), and Step 4 requires an IFT for all M-dimensional dyadic
blocks, then
j 1 = 1 B 1 j 2 = 1 B 2 j M = 1 B M P 1 , j 1 , P 2 , j 2 P m , j M
log ( P 1 , 1 P 2 , j 2 P m , j M ) = j 1 = 1 B 1 j 2 = 1 B 2 j M =
1 B M f 1 , j 1 f 2 , j 2 f M , j M N 1 N 2 N M log ( f 1 , j 1 f 2
, j 2 f M , j M N 1 N 2 N M ) = N 1 N 2 N M j 1 = 1 B 1 j 2 = 1 B 2
j M = 1 B M f 1 , j 1 f 2 , j 2 f M , j M log ( f 1 , j 1 f 2 , j 2
f M , j M ) + N 1 N 2 N M log ( N 1 N 2 N M ) j 1 = 1 B 1 j 2 = 1 B
2 j M = 1 B M f 1 , j 1 f 2 , j 2 f M , j M ( 19 ) ##EQU00037##
[0143] The fractional terms f.sub.m,j are within 0 to 1, so the
first multi-summation term .ltoreq.0, whereas the second
multi-summation is unity by invoking Eq. (18) repeatedly. Thus, the
order of Step 4 is necessarily .ltoreq.N.sub.1N.sub.2 . . . N.sub.M
log(N.sub.1N.sub.2 . . . N.sub.M). Since the order of the algorithm
is the highest order for the 4 steps (namely Step 1), then the
computational order of the FGFT algorithm is the same as that for
the M-dimensional FFT, i.e., N.sub.1N.sub.2 . . . N.sub.M
log(N.sub.1N.sub.2 . . . N.sub.M).
Inverse FGFT Algorithm
[0144] From the calculated multi-resolution FGFT matrix .PSI.(*),
the M-dimensional signal can be easily recovered by performing
these simple steps: [0145] 1. Perform the in-place, M-dimensional
FTs for each block of .PSI.(*) to generate .alpha.(*) [0146] 2.
Synthesize the M-dimensional, "blocked" weighting function [0147]
3. Divide .alpha.(*) from Step 1 by the weighting function in Step
2 to get S(k.sub.1, . . . , k.sub.M) [0148] 4. Take the
M-dimensional IFT of S(k.sub.1, . . . , k.sub.M) in Step 3 to
produce signal s(x.sub.1, . . . , x.sub.M)
[0149] Since the steps are completely equivalent to the forward
FGFT, albeit in reverse order, the computational order/burden of
the inverse FGFT is identical to its forward counterpart.
Examples
Texture Analysis Images
[0150] Image analysis techniques usually operate based on the
intensity or color of pixels in an image. Texture analysis is a
complementary technique that instead measures patterns of intensity
or color. One way to measure texture in an image is via the
frequency spectrum. For example, fine, detailed textures have
frequency spectra with more power in the high frequencies compared
to smooth, homogenous textures. The Fast GFT offers a very powerful
means of measuring texture as it efficiently produces a frequency
spectrum for each pixel in the image.
[0151] FIG. 15 and FIG. 16 show two applications of texture
analysis of medical images. In FIG. 15, patients with a type of
brain tumor called oligodendroglioma were scanned using magnetic
resonance (MR) imaging. The resulting images were analyzed with the
fast GFT to determine any texture differences. Significant
differences were discovered between tumors with and without a
particular genetic abnormality: deletion of the 1p and 19q
chromosomal arms. The presence of this abnormality is known to
predict tumors that respond to chemotherapy and are associated with
a good response to radiotherapy and good patient prognosis. Fast
GFT texture analysis was able to detect 1p/19q deleted tumors with
a sensitivity of 93% and specificity of 96%. Radiologists and
neurologists examining the same series of images were only able to
identify 1p/19q deletion with a sensitivity of 70% and specificity
of 63%. By using texture analysis of MR images this important tumor
feature can be identified non-invasively, without the need for
brain biopsy surgery to obtain tissue samples for genetic
analysis.
[0152] FIG. 16 demonstrates the same fast GFT texture technique
applied to computed tomography (CT) images of the wrist bones of
normal, age matched controls and osteoporosis patients. The texture
analysis technique was able to identify subjects who would suffer
at least one bone fracture before a follow up doctor's appointment
with a sensitivity of 82% and specificity of 88%. On the same
dataset the gold standard for diagnosis of osteoporosis, dual x-ray
absorbtiometry, predicted fractures with a sensitivity of 60% and
specificity of 71%.
Sound Analysis
[0153] The fast GFT algorithm produces local frequency spectra very
efficiently. The speed of the algorithm makes real-time analysis of
many signals possible. An example of this is voice analysis,
demonstrated in FIG. 17. In this example, a volunteer was recorded
speaking the English vowel sounds a, e, i, o and u. Compact disc
(CD) quality audio recording produces 44,100 samples per second. In
this example, the fast GFT of each second of CD-quality audio took
less than 27 ms, which is more than adequate for real-time
analysis. In comparison, the equivalent S-transform algorithm took
more than 17 minutes to perform for every second of input
recording. Real time spectral audio analysis is critical for
applications including voice recognition and automatic
dictation.
General Signal Analysis and Filtering
[0154] Most acquired signals are non-stationary: their frequency
content changes with time. Examples include climate data such as
ocean level and wind patterns, seismic signals used in oil and gas
exploration and audio and electrical signals acquired for medical
diagnosis. Many of these signals have either a high temporal
sampling rate or a very long sampling period, both of which result
in a great number of data points to be analyzed. In most cases,
practical examples of these signals are impossible to transform
with the standard S-transform technique, but are easily handled by
the fast GFT algorithm. An example is included in FIGS. 8A-C, which
shows an electrocardiogram (ECG) signal and its spectrum, obtained
via the fast GFT algorithm.
[0155] ECG monitoring is routinely conducted over a period of
hours. In many cases, as with sound analysis, a real-time spectrum
is useful. In others, offline analysis of the entire acquisition is
desired. The full ECG acquisition from FIGS. 8A-C is two hours long
at a sampling rate of 250 Hz. This signal took less than 1.5 s to
transform with the fast GFT algorithm, and the transform occupies 3
MB of memory. Using the standard S-transform algorithm the
computation would take more than 27 days and occupy 47 TB of
storage space, which is clearly impractical.
Radio Communications
[0156] Simple radio transmissions, such as familiar AM and FM
radio, begin with a fixed carrier frequency and encode information
through modification of that carrier. In the case of AM radio, the
amplitude of the carrier is modified. In FM radio, the frequency is
changed slightly. In either case, since the frequency does not vary
much from the basic carrier frequency, a short time Fourier
transform is often adequate to characterize the radio transmission.
In contrast, many modern systems use frequency hopping or spread
spectrum techniques either to more efficiently utilize the
available radio spectrum or to make radio transmissions resistant
to eavesdropping. In this case, the frequency content of the
transmission varies significantly and a technique offering
progressive resolution, such as the fast GFT, may be more useful
for analysis of the transmission.
[0157] All of the methods disclosed and claimed herein can be made
and executed without undue experimentation in light of the present
disclosure. While the apparatus and methods of this invention have
been described in terms of preferred embodiments, it will be
apparent to those of skill in the art that variations may be
applied to the methods and in the steps or in the sequence of steps
of the method described herein without departing from the concept,
spirit and scope of the invention. In addition, modifications may
be made to the disclosed apparatus and components may be eliminated
or substituted for the components described herein where the same
or similar results would be achieved. All such similar substitutes
and modifications apparent to those skilled in the art are deemed
to be within the spirit, scope, and concept of the invention as
defined by the appended claims.
APPENDIX A
Derivation of Equation (6)
[0158] The .alpha.-domain is defined in Eq. (5) as the Fourier
transform of S(.tau.,.nu.) along the .tau. axis. That is:
.alpha. ( v ' , v ) = .tau. { S ( .tau. , v ) } = .intg. - .infin.
+ .infin. S ( .tau. , v ) - 2 .pi. .tau.v ' .tau. = .intg. -
.infin. + .infin. { .intg. - .infin. + .infin. s ( t ) w ( .tau. -
t , .sigma. ) - 2 .pi. tv t } - 2 .pi. .tau.v ' .tau. = .intg. -
.infin. + .infin. s ( t ) - 2 .pi. tv t .intg. - .infin. + .infin.
w ( .tau. - t , .sigma. ) - 2 .pi. .tau.v ' t = .intg. - .infin. +
.infin. s ( t ) - 2 .pi. tv t [ W ( v ' , .sigma. ) - 2 .pi. tv ' ]
( Fourier Shift Theorem ) = W ( v ' , .sigma. ) .intg. + .infin. s
( t ) - 2 .pi. t ( v + v ' ) t = F ( v ' + v ) W ( v ' , .sigma. )
( 8 ) ##EQU00038##
where I.sub..tau. indicates the 1D Fourier transform over the .tau.
variable and F(.nu.) and W(.nu.',.sigma.) are the Fourier
transforms of the signal, s(t), and the window function,
w(.tau.,.sigma.), respectively.
APPENDIX B
TABLE-US-00004 [0159] TABLE 1 Window Kernel Transform w (t, .tau.,
v) .phi.(t, .tau., v) Fourier 1 e.sup.-i2.pi.vt (FT) (Allen and
Mills, 2004) Short-Time Fourier (STFT) (Allen and Mills, 2004) .PI.
[ C ( t - .tau. ) ] , .sigma. 2 .pi. e - ( t - .tau. ) 2 .sigma. 2
2 ##EQU00039## e.sup.-i2.pi.vt S- Transform (ST) (Stockwell et al.,
1996) v 2 .pi. e - ( t - .tau. ) 2 v 2 2 ##EQU00040##
e.sup.-i2.pi.vt Sinc Transform sin ( vt ) t ##EQU00041##
e.sup.-i2.pi.vt Morlet Wavelet (Goupillaud et al., 1984) v .PI. [ v
( t - .tau. ) ] ( 1 + e - .sigma. 2 - 2 e 3 4 .sigma. 2 ) 1 2 .pi.
1 4 e 1 2 v 2 ( t - .tau. ) 2 ##EQU00042## e i .sigma. v ( t -
.tau. ) - e 1 2 .sigma. 2 ##EQU00043## Haar Wavelet (Allen and
Mills, 2004) {square root over (|v|)}.PI.[v(t - .tau.)] { 1 0
.ltoreq. v ( t - .tau. ) < 1 2 , - 1 1 2 .ltoreq. v ( t - .tau.
) < 1 , 0 otherwise . } comb ( t , vt v max ) ##EQU00044## Six
continuous transforms and their corresponding kernels and windows
in the general transform (Eq. 9). The first four rows belong to the
Fourier family, using the kernel e.sup.-i2.pi..mu.t, while the
fifth and sixth rows are examples of common wavelet transforms. The
Haar wavelet kernel is normally finite in extent but here, for
consistency, it is convolved with a comb function to make it
periodic. Applying a boxcar window selects one particular period. C
and .sigma. are constants that determine the width of the window in
the STFT and Morlet wavelet, respectively. is the convolution
operator. The comb function is defined as comb ( t , T ) = n = -
.infin. .infin. .delta. ( t - nT ) . ##EQU00045##
* * * * *