Frequency-Time Domain Transformation

Any signal can completely be described in time or in frequency domain. As both representations are equivalent, it is possible to transform them into each other. This is done by the so-called Fourier Transformation and the inverse Fourier Transformation, respectively:

$$\qquad \underline{U}(j\omega) = \int\limits_{-\infty}^{\infty} u(t)\cdot e^{-j\omega\cdot t} \; dt$$ Fourier Transformation: (15.178)

$$\qquad u(t) = \frac{1}{2\pi} \cdot \int\limits_{-\infty}^{\infty} \underline{U}(j\omega)\cdot e^{j\omega\cdot t} \; d\omega$$ inverse Fourier Transformation: (15.179)

In digital systems the data $u(t)$ or $\underline{U}(j\omega)$, respectively, consists of a finite number $N$ of samples $u_k$ and $\underline{U}_n$. This leads to the discrete Fourier Transformation (DFT) and its inverse operation (IDFT):

$$\qquad \underline{U}_n = \sum_{k=0}^{N-1} u_k\cdot \exp\left( -j\cdot n\frac{2\pi\cdot k}{N} \right)$$ DFT: (15.180)

$$\qquad u_k = \frac{1}{N} \cdot \sum_{n=0}^{N} \underline{U}_n\cdot \exp\left( j\cdot k\frac{2\pi\cdot n}{N} \right)$$ IDFT: (15.181)

The absolute time and frequency values do not appear anymore in the DFT. They depend on the sample frequency $f_T$ and the number of samples $N$.

$$\Delta f = \frac{1}{N\cdot\Delta t} = \frac{f_T}{N}$$ (15.182)

Where $\Delta t$ is distance between time samples and $\Delta f$ distance between frequency samples.

With DFT the $N$ time samples are transformed into $N$ frequency samples. This also holds if the time data are real numbers, as is always the case in "real life": The complex frequency samples are conjugate complex symmetrical and so equalizing the score:

$$\underline{U}_{N-n} = \underline{U}_n^*$$ (15.183)

That is, knowing the input data has no imaginary part, only half of the Fourier data must be computed.

Fast Fourier Transformation

As can be seen in equation 15.180 the computing time of the DFT rises with $N^2$. This is really huge, so it is very important to reduce the time consumption. Using a strongly optimized algorithm, the so-called Fast Fourier Transformation (FFT), the DFT is reduced to an $N\cdot\log_2 N$ time rise. The following information stems from [61], where the theoretical background is explained comprehensively.

The fundamental trick of the FFT is to cut the DFT into two parts, one with data having even indices and the other with odd indices:

$$\underline{U}_n = \sum_{k=0}^{N-1} u_k\cdot \exp\left( -j\cdot n\frac{2\pi\cdot k}{N} \right)$$ (15.184)

$$= \sum_{k=0}^{N/2-1} u_{2k}\cdot \exp\left( -j\cdot n\frac{2\pi\cdo... ...0}^{N/2-1} u_{2k+1}\cdot \exp\left( -j\cdot n\frac{2\pi\cdot (2k+1)}{N} \right)$$ (15.185)

$$= \underbrace{ \sum_{k=0}^{N/2-1} u_{2k}\cdot \exp\left( -j\cdot n\... ...} u_{2k+1}\cdot \exp\left( -j\cdot n\frac{2\pi\cdot k}{N/2} \right) }_{F_{odd}}$$ (15.186)

$$W_{n,N} = \exp\left( 2\pi\cdot j\cdot \frac{n}{N} \right)$$ with (15.187)

The new formula shows no speed advantages. The important thing is that the even as well as the odd part each is again a Fourier series. Thus the same procedure can be repeated again and again until the equation consists of $N$ terms. Then, each term contains only one data $u_k$ with factor $e^0=1$. This works if the number of data is a power of two (2, 4, 8, 16, 32, ...). So finally, the FFT method performs $\log_2 N$ times the operation

$$u_{k1,even} + W_{n,x}\cdot u_{k2,odd}$$ (15.188)

to get one data of $\underline{U}_n$. This is called the Danielson-Lanzcos algorithm. The question now arises which data values of $u_k$ needs to be combined according to equation (15.188). The answer is quite easy. The data array must be reordered by the bit-reversal method. This means the value at index $k_1$ is swapped with the value at index $k_2$ where $k_2$ is obtained by mirroring the binary number $k_1$, i.e. the most significant bit becomes the least significant one and so on. Example for $N=8$:

000	$\leftrightarrow$	000		011	$\leftrightarrow$	110		110	$\leftrightarrow$	011
001	$\leftrightarrow$	100		100	$\leftrightarrow$	001		111	$\leftrightarrow$	111
010	$\leftrightarrow$	010		101	$\leftrightarrow$	101

Having this new indexing, the values to combine according to equation 15.188 are the adjacent values. So, performing the Danielson-Lanzcos algorithm has now become very easy.

Figure 15.4 illustrates the whole FFT algorithm starting with the input data $u_k$ and ending with one value of the output data $\underline{U}_n$.

Figure 15.4: principle of a FFT with data length 8

This scheme alone gives no advantage. But it can compute all output values within, i.e. no temporary memory is needed and the periodicity of $W_{n,N}$ is best exploited. To understand this, let's have a look on the first Danielson-Lanczos step in figure 15.4. The four multiplications and additions have to be performed for each output value (here 8 times!). But indeed this is not true, because $W_{n,2}$ is 2-periodical in $n$ and furthermore $W_{n,2} = -W_{n+1,2}$. So now, $u_0 + W_{0,2}\cdot u_4$ replaces the old $u_0$ value and $u_0 - W_{0,2}\cdot u_4$ replaces the old $u_4$ value. Doing this for all values, four multiplications and eight additions were performed in order to calculate the first Danielson-Lanczos step for all (!!!) output values. This goes on, as $W_{n,4}$ is 4-periodical in $n$ and $W_{n,4} = -W_{n+2,4}$. So this time, two loop iterations (for $W_{n,4}$ and for $W_{n+1,4}$) are necessary to compute the current Danielson-Lanczos step for all output values. This concept continues until the last step.

Finally, a complete FFT source code in C should be presented. The original version was taken from [61]. It is a radix-2 algorithm, known as the Cooley-Tukey Algorithm. Here, several changes were made that gain about 10% speed improvement.

There are many other FFT algorithms mainly aiming at higher speed (radix-8 FFT, split-radix FFT, Winograd FFT). These algorithms are much more complex, but on modern processors with numerical co-processors they gain no or hardly no speed advantage, because the reduced FLOPS are equalled by the far more complex indexing.

Real-Valued FFT

All physical systems are real-valued in time domain. As already mentioned above, this fact leads to a symmetry in frequency domain, which can be exploited to save 50% memory usage and about 30% computation time. Rewriting the C listing from above to a real-valued FFT routine creates the following function. As this scheme is not symmetric anymore, an extra procedure for the inverse transformation is needed. It is also depicted below.

More-Dimensional FFT

A standard Fourier Transformation is not useful in harmonic balance methods, because with multi-tone excitation many mixing products appear. The best way to cope with this problem is to use multi-dimensional FFT.

Fourier Transformations in more than one dimension soon become very time consuming. Using FFT mechanisms is therefore mandatory. A more-dimensional Fourier Transformation consists of many one-dimensional Fourier Transformations (1D-FFT). First, 1D-FFTs are performed for the data of the first dimension at every index of the second dimension. The results are used as input data for the second dimension that is performed the same way with respect to the third dimension. This procedure is continued until all dimensions are calculated. The following equations shows this for two dimensions.

$$\underline{U}_{n1,n2} = \sum_{k_2=0}^{N_2-1} \sum_{k_1=0}^{N_1-1} u_{k_1,k_2}\cdot \exp\l... ...k_1}{N_1} \right) \cdot \exp\left( -j\cdot n_2\frac{2\pi\cdot k_2}{N_2} \right)$$ (15.189)

$$= \sum_{k_2=0}^{N_2-1} \exp\left( -j\cdot n_2\frac{2\pi\cdot k_2}{N... ...2}\cdot \exp\left( -j\cdot n_1\frac{2\pi\cdot k_1}{N_1} \right) }_\text{1D-FFT}$$ (15.190)

Finally, a complete $n$-dimensional FFT source code should be presented. It was taken from [61] and somewhat speed improved.

Parameters:

ndim	-	number of dimensions
num[]	-	array containing the number of complex samples for every dimension
data[]	-	array containing the data samples,
		real and imaginary part in alternating order (length: 2*sum of num[]),
		going through the array, the first dimension changes least rapidly !
		all subscripts range from 1 to maximum value !
isign	-	is 1 to calculate FFT and -1 to calculate inverse FFT