Fast Algorithms for Multidimensional Signals

A Mathematical Expression is sufficient to represent an algorithm which in most cases represents input/output relationships. The efficiency of such an Algorithm can be evaluated by the amount of computational resources it takes to compute output or the quantity of interest.

Motivation and applications

In an algorithm we can assume that we have a mathematical expression relating the output and input. Example of Algorithms can be Filters, Fourier Transforms, Histograms, Image Enhancements and etc. All of these operations can be expressed as mathematical formulas and can be directly computed as we write the expressions. As one would expect this can be termed as the Obvious Implementation of the Algorithm. Also it is not necessary that this obvious implementation is always efficient.

When people began to compute such things, they began to look for more efficient ways. The wiki page aims at showcasing such efficient and Fast Algorithms for Multidimensional Signals. A multidimensional (M-D) signal can be modeled as a function of M independent variables, where M is greater than or equal to 2. These signals may be categorized as continuous, discrete, or mixed. A continuous signal can be modeled as a function of independent variables which range over a continuum of values, example – an audio wave travelling in space, 3-D space waves measured at different times. A discrete signal, on the other hand, can be modeled as a function defined only on a set of points, such as the set of integers. An Image is the simplest example of 2-D Signal which is spatial in nature.

In the context of Fast Algorithms, consider the below example:

We need to compute A which is given by

A = αγ + αδ + βγ + βδ where α,β,γ and δ are complex variables.

To compute A, we need 4 complex multiplications and 3 complex additions. The above equation can be written in its simplified form as

A = (α + β)(γ + δ)

This form requires only 1 complex multiplication and 2 complex additions.

Thus the second way of computing A is much more efficient and fast compared to the first method of computing A. This is the motivation for the evolution of the Fast Algorithms in the Digital Signal Processing Field. Consequently, many of the real life applications make use of these efficient Algorithms for fast computations.

Problem Statement and Basics

The simplest form of representing a system is through its Impulse response. And the output of such a system is given by the convolution of its input signal and system’s Impulse response. This is mathematically represented as follows:

Convolution of a Discrete Domain Input Signal with Discrete domain Impulse response

$y\left(n_{1},n_{2}\right)=\sum _{l_{1}=-\infty }^{+\infty }\sum _{l_{2}=\infty }^{+\infty }x(l_{1},l_{2})h(n_{1}-l_{1},n_{2}-l_{2})$

Where $h(n_{1},n_{2})$ is the impulse response of the system.

According to the equation above to obtain the Output value at a particular point ( say $y(0,0)$) we need to multiply several values of the input $x(l_{1},l_{2})$ and Impulse Response $h(-l_{1},-l_{2})$. Of course this is dependent on the Region of Support of the input as well as the impulse response. The key point here to be noted is that we need to perform so many complex multiplications and additions to obtain 1 output value.

Assuming a 2-D input signal is of length $M\times M$ and the systems impulse response is of length $N\times N$ we need to perform $M^{2}N^{2}$ number of multiplications to obtain complete output values. This is a huge number and very inefficient way of calculating.

We encounter a similar scenario when we have to compute the Discrete Fourier Transforms of a Signal of interest.

The Direct Calculation of 2-D DFT is simply the evaluation of the double Sum

$X\left(k_{1},k_{2}\right)=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}x(n_{1},n_{2})W_{N_{1}}^{n_{1}k_{1}}W_{N_{2}}^{n_{2}k_{2}}$ ${\text{ for }}0\leq k_{1}\leq N_{1}-1{\text{ and }}0\leq k_{2}\leq N_{2}-1$

${\text{ where }}W_{N}\equiv {\text{exp}}\left(-j{\frac {2\pi }{N}}\right)$

The total number of complex multiplications and complex additions needed to evaluate DFT by direct calculation is $N_{1}^{2}N_{2}^{2}$. This is a naive approach, however, we already know that an N-point 1-D DFT can be computed with far fewer than $N^{2}$ multiplications by using fast Fourier Transform (FFT) algorithm. As described in the next section we can develop fast Fourier transforms for calculating 2-D or higher Dimensional DFTs as well

Existing Approaches

Row Column Decomposition approach for the evaluation of DFT

The DFT sum $X\left(k_{1},k_{2}\right)$ in the previous equation can also be written in the following form

$X\left(k_{1},k_{2}\right)=\sum _{n_{1}=0}^{N_{1}-1}\left[\sum _{n_{2}=0}^{N_{2}-1}x(n_{1},n_{2})W_{N_{2}}^{n_{2}k_{2}}\right]W_{N_{1}}^{n_{1}k_{1}}$

Let $G\left(n_{1},k_{2}\right)$ denote the quantity inside the brackets and is given by:

$G\left(n_{1},k_{2}\right)=\sum _{n_{2}=0}^{N_{2}-1}x(n_{1},n_{2})W_{N_{2}}^{n_{2}k_{2}}$

$X\left(k_{1},k_{2}\right)=\sum _{n_{1}=0}^{N_{1}-1}G\left(n_{1},k_{2}\right)W_{N_{1}}^{n_{1}k_{2}}$

Employing this method DFT of $X$ can be computed as multiple 1-D DFTs. That is, each column of $G$ can be considered as 1-D DFT of the corresponding column of $x$($n_{1}$ = constant). And each row of $X$ is the 1-DFT of the corresponding row of the $G$($n_{2}$ = constant). Hence we are computing 2-D DFT by decomposing into Row and Column DFTs.

The same principle is employed for evaluating a M - dimensional signal.

Now let's talk about the computational savings we get using this approach. It is observed that we require $N_{1}N_{2}(N_{1}+N_{2})$ complex additions and multiplications. Further, if each of these 1-D DFTs are computed using 1-D FFT number of complex multiplications can be further reduced to $N_{1}N_{2}{\frac {\log _{2}N_{1}N_{2}}{2}}$

Vector Radix Fast Fourier Transform

Just like 1-D FFT, decimation in time can be achived in case of 2-D Signals. In 1-D DFT a signal of length of power of 2,the corresponding DFT can be expressed in terms of two half length DFTs, each of these can again be expressed as combination of quarter length DFTs and so on.

In case of 2-D Signals we can express $(N_{1}{\text{ x }}N_{2})$ DFT in terms of four ${\frac {N_{1}}{2}}{\text{ x }}{\frac {N_{2}}{2}}$ DFTs (assuming $N_{1}$ and $N_{2}$ are of length of power of 2). For the sake of simplicity let us assume that $N_{1}=N_{2}=N$. The DFT double sum can be decomposed into four separate summations, one over those samples of $x$ for which both $n_{1}$ and $n_{2}$ are even, one for which $n_{1}$ is even and $n_{2}$ is odd, one for which $n_{1}$ is odd and $n_{2}$ is even and the last one for which $n_{1}$ and $n_{2}$ are odd.

This is written as :

$X\left(k_{1},k_{2}\right)=S_{00}(k_{1},k_{2})+S_{01}(k_{1},k_{2})W_{N}^{k_{2}}+S_{10}(k_{1},k_{2})W_{N}^{k_{1}}+S_{11}(k_{1},k_{2})W_{N}^{k_{1}+k_{2}}$

where

$S_{00}(k_{1},k_{2})\equiv \sum _{m_{1}=0}^{N/2-1}\sum _{m_{2}=0}^{N/2-1}x(2m_{1},2m_{2})W_{N}^{2m_{1}k_{1}+2m_{2}k_{2}}$

$S_{01}(k_{1},k_{2})\equiv \sum _{m_{1}=0}^{N/2-1}\sum _{m_{2}=0}^{N/2-1}x(2m_{1},2m_{2}+1)W_{N}^{2m_{1}k_{1}+2m_{2}k_{2}}$

$S_{10}(k_{1},k_{2})\equiv \sum _{m_{1}=0}^{N/2-1}\sum _{m_{2}=0}^{N/2-1}x(2m_{1}+1,2m_{2})W_{N}^{2m_{1}k_{1}+2m_{2}k_{2}}$

$S_{11}(k_{1},k_{2})\equiv \sum _{m_{1}=0}^{N/2-1}\sum _{m_{2}=0}^{N/2-1}x(2m_{1}+1,2m_{2}+1)W_{N}^{2m_{1}k_{1}+2m_{2}k_{2}}$

All the arrays $S_{00}S_{01}S_{10}$ and $S_{11}$ are the each periodic in $(k_{1},k_{2})$ with horizontal and vertical periods ${\frac {N}{2}}$. Using this fact and also the fact that $W_{N}^{N/2}=-1,$ we can obtain the following identities :

$X\left(k_{1},k_{2}\right)=S_{00}(k_{1},k_{2})+S_{01}(k_{1},k_{2})W_{N}^{k_{2}}+S_{10}(k_{1},k_{2})W_{N}^{k_{1}}+S_{11}(k_{1},k_{2})W_{N}^{k_{1}+k_{2}}$

$X\left(k_{1}+{\frac {N}{2}},k_{2}\right)=S_{00}(k_{1},k_{2})+S_{01}(k_{1},k_{2})W_{N}^{k_{2}}-S_{10}(k_{1},k_{2})W_{N}^{k_{1}}-S_{11}(k_{1},k_{2})W_{N}^{k_{1}+k_{2}}$

$X\left(k_{1},k_{2}+{\frac {N}{2}}\right)=S_{00}(k_{1},k_{2})-S_{01}(k_{1},k_{2})W_{N}^{k_{2}}+S_{10}(k_{1},k_{2})W_{N}^{k_{1}}-S_{11}(k_{1},k_{2})W_{N}^{k_{1}+k_{2}}$

$X\left(k_{1}+{\frac {N}{2}},k_{2}+{\frac {N}{2}}\right)=S_{00}(k_{1},k_{2})-S_{01}(k_{1},k_{2})W_{N}^{k_{2}}-S_{10}(k_{1},k_{2})W_{N}^{k_{1}}+S_{11}(k_{1},k_{2})W_{N}^{k_{1}+k_{2}}$

The above equation tell us how to compute the four DFT points $X\left(k_{1},k_{2}\right)X\left(k_{1}+{\frac {N}{2}},k_{2}\right)X\left(k_{1},k_{2}+{\frac {N}{2}}\right){\text{ and }}X\left(k_{1}+{\frac {N}{2}},k_{2}+{\frac {N}{2}}\right)$ for a particular value of $(k_{1},k_{2})$ from the four points $S_{00}(k_{1},k_{2})S_{01}(k_{1},k_{2})S_{10}(k_{1},k_{2}){\text{ and }}S_{11}(k_{1},k_{2})$. And these corresponding samples of $S_{00}(k_{1},k_{2})$ can be obtained by evaluating the $\left({\frac {N}{2}}\times {\frac {N}{2}}\right)$ -point DFT (similarly other $S_{ij}$ arrays can be obtained).

Thus we see that the $N\times N$ DFT can be expressed in terms of four ${\frac {N}{2}}\times {\frac {N}{2}}$ DFTs.

By analogy from 1-D case, the computation depicted in the below figure is called a ${\text{butterfly}}$ or more precisely $radix-(2\times 2){\text{ butterfly}}$ .

Each butterfly requires three complex multiplications and eight complex additions for the calculation of outputs from the inputs. And to compute all the samples of $X$ from $S_{00}(k_{1},k_{2})S_{01}(k_{1},k_{2})S_{10}(k_{1},k_{2}){\text{ and }}S_{11}(k_{1},k_{2})$ requires calculations of ${\frac {N^{2}}{4}}$ butterflies.

This decimation procedure is performed $\log _{2}N$ times when $N$ is a power of 2. Each stage of decimation consists of ${\frac {N^{2}}{4}}$ butterflies, and each butterfly involves three complex multiplications and eight complex additions and hence the number of complex multiplications that needs to be performed during the computation of an $(N\times N)$ -point radix $(2\times 2)$ FFT is given by

$C_{vectorRadix(2\times 2)}={\frac {3N^{2}}{4}}\log _{2}N$

Small Radix Cooley-Tukey algorithms

A two-dimensional discrete Fourier transform can be computed by applying the Cooley–Tukey FFT first to each row and then to each column. This can be regarded simply as a parenthesization of the equation of the two-dimensional Fourier transform with either the row sum or the column sum within the inner parentheses. Let $v$be an array of elements $v_{i'i''}$ for $i'=0,...n'-1$ and $i''=0,...n''-1$ from the field F. The two-dimensional Fourier transform of v is another two dimensional array V of elements of F, given by

$V_{k'k''}=\sum _{i'=0}^{n'-1}\sum _{i''=0}^{n''-1}w^{i'k'}\mu ^{i''k''}v_{i'i''}$

$k'=0,...,n'-1{\text{ and }}k''=0,...,n''-1$

where ω is an nth root of unity in the field F and μ is an n"th root of unity in the field F, which ordinarily would be chosen equal to ω when n' equals n". It follows that

$V_{k'k''}=\sum _{i'=0}^{n'-1}w^{i'k'}\left[\sum _{i''=0}^{n''-1}\mu ^{i''k''}v_{i'i''}\right]$

$V_{k'k''}=\sum _{i''=0}^{n''-1}w^{i''k''}\left[\sum _{i'=0}^{n'-1}\mu ^{i'k'}v_{i'i''}\right]$

Hence we see that one can compute the two-dimensional Fourier transform by computing a one-dimensional Fourier transform, either first along every column then along every row, or the other way around. Any FFT algorithm can be used on the rows and any FFT algorithm, possibly a different one, can be used on the columns. Many good FFT algorithms are available, and any one of them can be chosen for the purpose of reducing the number of multiplications and the number of addition.

Let $M(n\times n)$ be the number of multiplications in the field F needed by this algorithm to compute an n by n-point Fourier transform with n a power of two. It satisfies the recursion

$M(n\times n)=4M\left({\frac {n}{2}}\times {\frac {n}{2}}\right)+{\frac {3}{4}}n^{2}$

And this recursion is solved by

$M(n\times n)={\frac {3}{4}}n^{2}\left(\log _{2}n-C\right)$, where C is a constant chosen to fit the number of multiplications