Now you know enough to be able to implement DFT and IDFT. The language of choice will be C.
This implementation is made with readability in mind and to be easy to understand. Speed is of no concern here.
All code written here is licensed under Creative Commons Zero license, whose full text can be found at https://creativecommons.org/publicdomain/zero/1.0/. It basically means that you can do with the code whatever you like without asking for any permissions.
Let's start by implementing complex numbers. Altough they are part of C99 standard, their support between compilers differs and is not always available.
First structure to store complex numbers.
typedef struct
{
float r; // real part
float i; // imaginary part
} cnum;
Now time for operations on complex numbers: addition and multiplication
// complex numbers addition a+b
cnum cnum_add(cnum a, cnum b)
{
cnum out;
out.r = a.r + b.r;
out.i = a.i + b.i;
return out;
}
// complex numbers multiplication a*b
cnum cnum_mul(cnum a, cnum b)
{
cnum out;
out.r = (a.r * b.r) - (a.i * b.i);
out.i = (a.r * b.i) + (a.i * b.r);
return out;
}
Now that we have complex numbers we can implement DFT. Let's start with forward DFT. For reminder
$$x = (x_0, x_1, x_2, ..., x_{N-1})$$
$$DFT(x) = y = (y_0, y_1, y_2, ..., y_{N-1})$$
$$y_k = \sum_{n=0}^{N-1}x_n*e^{-j*2\pi*k/N*n} = \sum_{n=0}^{N-1}x_n* (cos(2\pi*k/N*n) - isin(2\pi*k/N*n))$$
So DFT function takes a list x of length N and returns a list y of length N. Both lists contain complex numbers.
#include <math.h>
// C99 standard does not define constant PI
#define PI 3.1415926218032837
void dft(cnum *x, /*out*/cnum *y, int N)
{
// make sure y is zeroed before use
for(int n=0; n<N; ++n)
{
y[n].r = 0.0;
y[n].i = 0.0;
}
float Nf = (float) N;
// for complex sinusoid of k-th frequency
for(int k=0; k<N; ++k)
{
// project signal onto that complex sinusoid
for(int n=0; n<N; ++n)
{
cnum csin;
csin.r = cosf(2*PI*k/Nf*n);
csin.i = -sinf(2*PI*k/Nf*n); // remember about minus
cnum addend = cnum_mul(x[n], csin);
y[k] = cnum_add(y[k], addend);
}
}
}
There are few elements that require some explanation or additional attention.
We include math.h to have sinf() and cosf() functions.
Because C99 standard does not define constant PI (and M_PI is a POSIX or GNU extension) we shall define it ouselves.
Next, before doing anything, values of y are zeroed as we will be using them to store intermediary values of the sum.
Next we explicitly cast N value to type float. This is because we are computing value of cos(2*PI*k/N*n) and sin(2*PI*k/N*n). If N were integer (and compiler did a very strange thing) there could happen integer division, which would ruin the algorithm. To make sure it never happens just force N to be used as a float there and problem solved. Better safe than sorry.
Last thing that you shall remember is to put minus sign before sinf() function. If your code gives results with reversed signs in imaginary numbers (eg. 2+2j instead of 2-2j), then this is where you have made a mistake.
Reminding the formula
$$x_n = \frac{1}{N} \sum_{k=0}^{N-1}y_k * e^{j*2\pi*k/N*n}$$
we see that inverse DFT is almost the same as forward DFT, so their implementations will also be alike.
void idft(cnum *y, /*out*/ cnum *x, int N)
{
// make sure x is zeroed before use
for(int n=0; n<N; ++n)
{
x[n].r = 0.0;
x[n].i = 0.0;
}
float Nf = (float) N;
// for a given point
for(int n=0; n<N; ++n)
{
// in every complex sinusoid
for(int k=0; k<N; ++k)
{
cnum csin;
csin.r = cosf(2*PI*k/Nf*n);
csin.i = sinf(2*PI*k/Nf*n); // no minus this time
cnum addend = cnum_mul(y[k], csin);
x[n] = cnum_add(x[n], addend);
}
}
// divide points by N
for(int n=0; n<N; ++n)
{
x[n].r /= Nf;
x[n].i /= Nf;
}
}
The obvious difference is that this time y is input, x is output. And so x is zeroed before anything happens.
Next internal and external loops are swapped.
In forward DFT projecting signal x onto k-th csin means that k is set and loop is walking over subsequent
points of that csin (by incrementing values of n).
In inverse DFT getting value of a sample n means that n is set and loop is walking over csins with subsequent
frequencies (by incrementing values of k).
In IDFT there is no minus sign in exponential and so none before sinf() function.
Values of IDFT shall be divided by length of a signal to give original signal in time domain. This is dealt with a loop over all the samples at the end of the function.
Time to write main function and make sure that everything works as it should. If your code happens not to work, check places mentioned above.
#include <stdio.h>
int main(void)
{
cnum x[4], y[4], restored_x[4];
// create some input signal
x[0].r = 6.0;
x[0].i = 0.0;
x[1].r = 3.0;
x[1].i = 0.0;
x[2].r = 5.0;
x[2].i = 0.0;
x[3].r = 1.0;
x[3].i = 0.0;
// compute forward dft
dft(x, y, 4);
// display result of forward dft
printf("%f%+fi\n", y[0].r, y[0].i);
printf("%f%+fi\n", y[1].r, y[1].i);
printf("%f%+fi\n", y[2].r, y[2].i);
printf("%f%+fi\n", y[3].r, y[3].i);
// compute inverse dft
idft(y, restored_x, 4);
// separate result of dft and idft with empty line
printf("\n");
// display result of inverse dft
printf("%f%+fi\n", restored_x[0].r, restored_x[0].i);
printf("%f%+fi\n", restored_x[1].r, restored_x[1].i);
printf("%f%+fi\n", restored_x[2].r, restored_x[2].i);
printf("%f%+fi\n", restored_x[3].r, restored_x[3].i);
return 0;
}
Remember that computations done with floating point numbers are often approximate, so don't expect
that x[0] == restored_x[0]. But you shall expect that difference between them to be about $10^{-8}$
in case of float (which is a level of a round off error). To see this you may add following code
at the end of main.
// show the difference between values of original and reconstructed signal
printf("\n");
printf("%d %d\n", x[2].r == restored_x[2].r, x[2].i == restored_x[2].i);
printf("%e, %e\n", x[2].r - restored_x[2].r, x[2].i - restored_x[2].i);
To get a complete code download a source file.
You may also get Makefile, although it does not contain much.
Using cc (gcc) as a compiler you may compile the code using
cc -Wall -Wextra -std=c99 -o dft dft.c -lm
-Wall -Wextra -std=c99 are to make sure we are not doing stupid things and are following c99
standard.
Don't forget to put option -lm to link program with mathematics library. Otherwise
you will get error about unknown symbols sinf, cosf.
Program does not take any parameters, so running it is just simple ./dft
[noniwoo src]$ ./dft 15.000000+0.000000i 1.000000-2.000001i 7.000000+0.000001i 1.000000+2.000000i 6.000000+0.000000i 3.000000-0.000000i 5.000000-0.000000i 1.000000-0.000000i 1 0 0.000000e+00, 8.940697e-08next part