13,765,976 members
alternative version

#### Stats

24.8K views
58 bookmarked
Posted 9 May 2013
Licenced CPOL

# Quick Fourier Transformation

, 9 May 2013
Some ideas to make the Discrete Fourier Transformation a bit quicker and implemented a lean version of the DFT algorithm.

## Basic DFT algorithm with complex numbers

According to the main clause of Fourier a signal f(n) that is periodic within N and we have N samples of it can be transformed into its Fourier components as follows:

with

Another equation for the complex Fourier components ck looks like

and

That means the complex Fourier component ck is a complex sum of complex products of f(n) and complex vectors spinning around the unit circle. For a complete Fourier transformation we would need N Fourier components (ck) and for each component we need N sine and N cosine values. For 1000 samples that would be 1 million sine and cosine values to be calculated. That consumes a lot of CPU power and time. But it works J

```public void CalcDFT(int N)
{
int k, n;
TKomplex w;
if (N > 0)
{
for (k = 0; k < N; k++)
{
c[k].real = 0;
c[k].imag = 0;
for (n = 0; n < N; n++)
{
w.real = Math.Cos((double)(2.0 * Math.PI * (double)(k * n) / (double)(N)));
w.imag = -Math.Sin((double)(2.0 * Math.PI * (double)(k * n) / (double)(N)));
c[k] = ksum(c[k], kprod(w, y[n]));
}
c[k].real = c[k].real / (double)(N) * 2.0;
c[k].imag = -c[k].imag / (double)(N) * 2.0;
}
}
c[0].real = c[0].real / 2;
c[0].imag = c[0].imag / 2;
}```

Calculated by this standard implementation a Fourier transformation of 1000 samples took 140 ms on my computer (there are certainly faster computers on this world J).

The source of this standard implementation can be found in the Visual C# project DFT_std.

`ksum()` and `kprod()` are simple functions for complex addition and multiplication.

```public TKomplex ksum(TKomplex a, TKomplex b)
{
TKomplex res;
res.real = a.real + b.real;
res.imag = a.imag + b.imag;
return (res);
}
public TKomplex kprod(TKomplex a, TKomplex b)
{
TKomplex res;
res.real = a.real * b.real - a.imag * b.imag;
res.imag = a.real * b.imag + a.imag * b.real;
return (res);
}```

Now, the sine and the cosine function both are periodic to 2π and for k = 1 we have to calculate both, sine and cosine for the values 0 to 2π with step size 2π /N. For k = 2 we have 0 to 4 π with step size 2/N and so on. So we get for k = 2 the same values as we got them already for k = 1. Only, we get each second one and when n gets bigger than N/2 we finished one round of the unit circle and start again at the same values as for n = 0. The bigger k becomes the bigger the steps become and the more rounds of the unit circle we pass. But the sine and the cosine values we get are still the same as we got them for k = 1. That means it’s basically enough just to calculate N sine and N cosine value between 0 and 2π with step size 2π/N and keep them in a look up table for the entire transformation. For 1000 samples we just have to calculate 1000 sine and 1000 cosine values instead of 999000 and that saves quite some time.

We build the look up table with:

```for (k = 1; k < N; k++)
{
we[k].real = Math.Cos((2.0 * Math.PI * (double)(k) / (double)(N)));
we[k].imag = -Math.Sin((2.0 * Math.PI * (double)(k) / (double)(N)));
}```

And the corresponding sine and cosine values we can obtain from the look up table by a simple modulo division of k * n by N. Like this the discrete Fourier transformation becomes not “fast” but quite “quick” and lean J

```public void CalcDFT()   //  Quick and lean Fourier transformation
{
int k, n;
if (N > 0)
{
for (k = 0; k < N; k++)
{
c[k].real = 0;
c[k].imag = 0;
for (n = 0; n < (N - 1); n++)
{
c[k] = ksum(c[k], kprod(we[(k * n) % N], y[n]));
}
c[k].real = c[k].real / N * 2;
c[k].imag = -c[k].imag / N * 2;
}
c[0].real = c[0].real / 2;
c[0].imag = c[0].imag / 2;
}
}```

This implementation takes 53 ms for the same 1000 samples as the above mentioned. Quite a bit quicker.

The code to this implementation is in DFT_lean.

## Complex numbers?

Now there came some more considerations into my mind. In the main clause there are no complex numbers. The input is a sequence of real numbers (at least in my application coming from electric signal analysis J) and we end up with another sequence of real numbers. So, why do we use complex numbers in the middle?

We have a complex multiplication and a complex addition in our algorithm. So let’s have a look at them. First the multiplication, we[x] * y[n]. That looks like:

`(we[x].re * y[n].re - we[x].im * y[n].im) + j(we[x].re * y[n].im + we[x].im * y[n].re)`

y[n].im is always 0 and therefore the second and the third part are also 0 and our complex multiplication reduces to

`(we[x].re * y[n].re) + j(we[x].im * y[n].re)`

And that means

`y[n].re * cosine(x) + j(y[n].re * sine(x))`

That looks quite like ak and bk of the main clause. In the following addition we just have to sum the real parts and the imaginary parts. So there is no mixing of imaginary and real part in the complete calculation and finally we get the ak part in the real part of ck and the bk part in the imaginary part of ck and we can forget about all the complex stuff J

With the limitation that the input sequence does not contain imaginary numbers we can implement the Fourier transformation without complex numbers:

```public void CalcDFT()
{
int k, n;
if (N > 0)
{
for (k = 0; k < N; k++)
{
a[k] = 0;
b[k] = 0;
for (n = 0; n < (N - 1); n++)
{
a[k] = a[k] + ((cosine[(k * n) % N] * y[n]));
b[k] = b[k] + ((sine[(k * n) % N] * y[n]));
}
a[k] = a[k] / N * 2;
b[k] = b[k] / N * 2;
}
a[0] = a[0] / 2;
b[0] = b[0] / 2;
}
}```

The initialization becomes like this

```public TFftAlgorithm(int order)
{
int k;
N = order;
y = new double[N + 1];
a = new double[N + 1];
b = new double[N + 1];
xw = new double[N + 1];
sine = new double[N + 1];
cosine = new double[N + 1];
cosine[0] = 1.0;    // we don't have to calculate cos(0) = 1
sine[0]   = 0;      //                        and sin(0) = 0
for (k = 1; k < N; k++) //  init vectors of unit circle
{
cosine[k] = Math.Cos((2.0 * Math.PI * (double)(k) / (double)(N)));
sine[k] = Math.Sin((2.0 * Math.PI * (double)(k) / (double)(N)));
}
}```

With this implementation the 1000 samples were transformed in 20 ms. That’s almost 5 times quicker than the first implementation at the top of this article. It’s not an fast Fourier transformation. But it’s at least it’s a quick Fourier transformation and it looks quite lean J

This implementation is in the project DFT_lean_quick.

With the rectangle signal the result should look like this:

The blue shape is the input sequence: A rectangle signal with 1000 samples. The magnitude is 20. The red shape is the inverse Fourier transformation with the first 30 Fourier components.

For more math visit www.mosismath.com :-)

## Share

 Software Developer (Senior) Güdel AG, Langenthal Switzerland
No Biography provided

## You may also be interested in...

 First Prev Next
 I don't know your Formula code Member 139016386-Jul-18 9:57 Member 13901638 6-Jul-18 9:57
 DFT Member 1072698922-Apr-15 5:38 Member 10726989 22-Apr-15 5:38
 Re: DFT Mosi_6222-Oct-15 9:49 Mosi_62 22-Oct-15 9:49
 Could be many times faster D. Christian Ohle7-Mar-15 7:12 D. Christian Ohle 7-Mar-15 7:12
 Re: Could be many times faster Mosi_628-Mar-15 2:51 Mosi_62 8-Mar-15 2:51
 Re: Could be many times faster D. Christian Ohle8-Mar-15 4:35 D. Christian Ohle 8-Mar-15 4:35
 poh asilarslan6-Dec-13 4:21 asilarslan 6-Dec-13 4:21
 My vote of 5 Mihai MOGA13-Jun-13 21:41 Mihai MOGA 13-Jun-13 21:41
 Awesome implementation Jorge Varas13-May-13 10:27 Jorge Varas 13-May-13 10:27
 Awesome implementation, thank you for sharing
 Quicker ? YvesDaoust13-May-13 3:08 YvesDaoust 13-May-13 3:08
 Quick or Fast? Phil Atkin13-May-13 1:49 Phil Atkin 13-May-13 1:49
 My vote of 5 11-May-13 15:40 11-May-13 15:40
 My vote of 5 Orlando Selenu9-May-13 23:54 Orlando Selenu 9-May-13 23:54
 Last Visit: 15-Nov-18 7:14     Last Update: 15-Nov-18 7:14 Refresh 1