Results 1 to 4 of 4

Math Help - Discrete Fourier Transform implementation

  1. #1
    Member
    Joined
    Nov 2008
    Posts
    92

    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))
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by algorithm View Post
    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
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Member
    Joined
    Nov 2008
    Posts
    92
    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?
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by algorithm View Post
    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
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Discrete Fourier Transform (DFT).
    Posted in the Advanced Math Topics Forum
    Replies: 1
    Last Post: November 27th 2010, 10:45 PM
  2. [SOLVED] Discrete fourier transform
    Posted in the Differential Geometry Forum
    Replies: 12
    Last Post: November 9th 2010, 10:19 AM
  3. Discrete Fourier Transform
    Posted in the Calculus Forum
    Replies: 0
    Last Post: September 10th 2010, 06:36 AM
  4. 2D discrete Fourier transform
    Posted in the Calculus Forum
    Replies: 7
    Last Post: May 13th 2010, 03:17 AM
  5. Eigenvectors for the discrete fourier transform
    Posted in the Advanced Applied Math Forum
    Replies: 1
    Last Post: December 2nd 2009, 03:16 AM

Search Tags


/mathhelpforum @mathhelpforum