# Thread: Discrete Fourier Transform implementation

1. ## Discrete Fourier Transform implementation

Hi,

I have defined a filter, X, with unit magnitude and a phase of 45deg. Its real and imaginary parts are:

real{X} = {0.70701, 0.70701, 0.70701, 0.70701}
imag{X} = {0.70701, 0.70701, 0.70701, 0.70701}

So:
magnitude{X} = {1, 1, 1, 1}
phase{X} = {45deg, 45deg, 45deg, 45deg}

Now, if I take the inverse DFT of the real/imagninary parts, and then take the DFT of that result, I should end up with the same real/imaginary parts:

X --> [IDFT] --> h[n] --> [DFT] --> X (h[n] is always real).

However, this is not what I'm getting when I perform the IDFT/DFT. What I get back is the correct real parts, but the imaginary parts are all zero.

The algorithms (see below) I'm using for the IDFT/DFT is are direct implementations of the exponential forms.

What am I not doing right?

Thanks.

DFT algorithm:
Code:
omega = 2pi / N

do while k < N
realX[k] = 0
imagX[k] = 0

do while n < N
realX[k] = realX[k] + h[n] * cos(k * omega * n)
imagX[k] = imagX[k] - h[n] * sin(k * omega * n)
n = n + 1

realX[k] = realX[k] / N;
imagX[k] = imagX[k] / N;
k = k + 1
IDFT algorithm:
Code:
do while n < N
do while k < N
h[n] = h[n] + (realX[k] * cos(k * omega * n) - imagX[k] * sin(k * omega * n))

2. Originally Posted by algorithm
Hi,

I have defined a filter, X, with unit magnitude and a phase of 45deg. Its real and imaginary parts are:

real{X} = {0.70701, 0.70701, 0.70701, 0.70701}
imag{X} = {0.70701, 0.70701, 0.70701, 0.70701}

So:
magnitude{X} = {1, 1, 1, 1}
phase{X} = {45deg, 45deg, 45deg, 45deg}

Now, if I take the inverse DFT of the real/imagninary parts, and then take the DFT of that result, I should end up with the same real/imaginary parts:

X --> [IDFT] --> h[n] --> [DFT] --> X (h[n] is always real).

However, this is not what I'm getting when I perform the IDFT/DFT. What I get back is the correct real parts, but the imaginary parts are all zero.

The algorithms (see below) I'm using for the IDFT/DFT is are direct implementations of the exponential forms.

What am I not doing right?

Thanks.

DFT algorithm:
Code:
omega = 2pi / N

do while k < N
realX[k] = 0
imagX[k] = 0

do while n < N
realX[k] = realX[k] + h[n] * cos(k * omega * n)
imagX[k] = imagX[k] - h[n] * sin(k * omega * n)
n = n + 1

realX[k] = realX[k] / N;
imagX[k] = imagX[k] / N;
k = k + 1
IDFT algorithm:
Code:
do while n < N
do while k < N
h[n] = h[n] + (realX[k] * cos(k * omega * n) - imagX[k] * sin(k * omega * n))
Without disecting what you are doing in detail it has an error somewhere:

Code:
>
>
>sp=0.70701*[1+I,1+I,1+I,1+I]
Column 1 to 2:
0.70701+0.70701i            0.70701+0.70701i
Column 3 to 4:
0.70701+0.70701i            0.70701+0.70701i
>help fftw1
function fftw1 (data,frwrd)
## Default for frwrd :
1
##
##  Wrapper function for FFTW.
##
##  This operates on row vectors only, there seems to
##  be a problem about converting an Euler complex matrix to
##  an array of fftw complex. Maybe one day this can be sorted out
##  there does not appear to be a large penalty in feeding one
##  row at a time to fftw.
##
##  frwrd=1 gives the forward transform in the Euler sense
##  all other values give the normalised backward transform
##
##  note fftw defines the forward and backward transforms
##  with the opposite exponents to Euler. in the DLL this
##  is corrected so that the calls have the same sense.
##
##  on complex data this is a factor of 2 down on Matlabs
##  fftw routines, for real data this is closer to a factor
##  of 4 due to the use of rfftw in Matlab.
##
##  RL 2002-06
##
>xx=fftw1(sp,0)
Column 1 to 2:
0.70701+0.70701i                        0+0i
Column 3 to 4:
0+0i                        0+0i
>yy=fftw1(xx,1)
Column 1 to 2:
0.70701+0.70701i            0.70701+0.70701i
Column 3 to 4:
0.70701+0.70701i            0.70701+0.70701i
>
The FFT algorithm is a fast implementation of the DFT, so what do you have for h ?)

(still supprised that the normalisation is such that fftw1(fftw1(data,0),1)=data )

CB

3. Hi,

Thanks for taking the fft to confirm the correct result.

The reason I'm using the unoptimised DFT is so I know exactly what's happening (this helps me locate the problem). When it's working I hope to use the FFT.

When I take the IDFT of the real/imaginary vectors, I obtain the time-signal h[n], whose values are:

h[n] = {2.8284, -2.220446049250313E-16, -1.1102230246251565E-16, 0.0}.

Those E-16 could probably be treated at zero, so:

h[n] = {2.8284, 0, 0, 0}.

It looks like my last three values agree with your result...but the first one (2.8284) is way off.

The algorithms I'm using are direct implementations of the discrete exponential forms. That makes me think the algorithm has to be correct, but it seems like it isn't.

Do you have any other suggestions?

4. Originally Posted by algorithm
Hi,

Thanks for taking the fft to confirm the correct result.

The reason I'm using the unoptimised DFT is so I know exactly what's happening (this helps me locate the problem). When it's working I hope to use the FFT.

When I take the IDFT of the real/imaginary vectors, I obtain the time-signal h[n], whose values are:

h[n] = {2.8284, -2.220446049250313E-16, -1.1102230246251565E-16, 0.0}.

Those E-16 could probably be treated at zero, so:

h[n] = {2.8284, 0, 0, 0}.

It looks like my last three values agree with your result...but the first one (2.8284) is way off.

The algorithms I'm using are direct implementations of the discrete exponential forms. That makes me think the algorithm has to be correct, but it seems like it isn't.

Do you have any other suggestions?
1. You have only computed the real component of the inverse transform (and the inverse transform is not real in this case as the frequency domain data does not have the required conjugate symmetry)

2. There is an arbitrary normalisation in the DFT/IDFT pair, and a factor of n (there size of the transform) difference in the magnitude is often noted (that is normalisation only applied to one of the transform pair, so a factor of 4 difference is not surprising).

CB