fftw - loosing accuracy in cosine transform

From: Kamaraju Kusumanchi (kk288_at_cornell.edu)
Date: 09/23/04 

Date: Thu, 23 Sep 2004 00:01:29 -0400

Hi all,
    Sorry if this is a wrong place for posting fftw related question.
But on the www.fftw.org, I could not find any email address other than
fftw@fftw.org which I think is meant to leave feedback rahter than
asking questions. In the past I got very good response on fft related
issues from this news group, so I am posting it here.

I am using fftw 3.0.1, absoft 8.0, fedora core 1. An example probably
explains my problem better. The output and sample code is given at the end.

My questions
(1) why I am getting only half the accuracy when doing a cosine
transform as opposed to the fouerier transform?
For example in fourier transform the round off error is O(1e-18) whereas
in cosine transform it is O(1e-9).
(2) Is it something that can be avoided or is it inherent in the cosine
transform algorithm?
Better still, can I improve the code to reduce the round-off errors in
the cosine transform?

  output of the fourier transform
   -4.92871531348332D-017 -6.22399809642460D-018
    1.00000000000000 -1.03683365493055D-016
    2.93123004154575D-017 1.25603286447767D-017
    7.51060464024562D-018 7.65378971138986D-018
    1.31628920003318D-017 1.45926836152971D-017
   -8.53809210832335D-018 2.77146745197110D-017
    1.54345126076430D-017 2.74725077800306D-018
    1.81322035841033D-017 7.65378971138986D-018
    6.22399809642460D-018 7.65378971138986D-018
    1.84314369322536D-018 -1.11783997592600D-018
   -1.26933184863306D-019 2.74725077800306D-018
   -7.51060464024562D-018 7.65378971138986D-018
   -7.14895807482630D-019 1.45926836152971D-017
   -5.98208548668877D-017 -1.47589455874086D-017
   -1.40047209926778D-017 1.25603286447767D-017
   -1.81322035841033D-017 7.65378971138986D-018
  output of the cosine/chebyshev transform
   0.499999985996426
   0.999999976990354
   0.500000025327541
    1.54879742012126D-008
   -5.75090005171598D-009
    3.17029677751287D-009
   -2.08134160024177D-009
    1.51237007364410D-009
   -1.17638936092401D-009
    9.62117874234755D-010
   -8.18547139920208D-010
    7.19508711332980D-010
   -6.50485120223127D-010
    6.02994678319307D-010
   -5.71958809144919D-010
    5.54385367040194D-010
   -2.74346247976232D-010

program end_of_spectrum
   use nrtype, WPC => DPC, WP => DP
   implicit none
   include "fftw3.f"

   integer, parameter :: N=16
   integer :: j, npo = N+1 ! npo = n plus one
   integer(I8B) :: plan_fou, plan_cheb
   complex(WPC) :: sample(0:N-1), output(0:N-1)
   real(WP) :: cheb_sample(0:N), cheb_output(0:N)
   ! sample array contains input to the FFTW program
   ! The output of fftw is stored in output array
   ! Initialize sample to a simple function such that the actual fourier
   ! coeffficients will be zero after some inital modes.
   ! Then by observing the end values of output we know what the
round-off errors
   ! from the fftw program are
   !
   ! cheb_sample and cheb_output contain the input and output to the
chebyshev
   ! transform.

   call dfftw_plan_dft_1d( plan_fou, N, sample, output, FFTW_FORWARD,
FFTW_PATIENT)
   call dfftw_plan_r2r_1d( plan_cheb, npo, cheb_sample, cheb_output, &
   FFTW_REDFT00, FFTW_PATIENT )

   ! initialize sample
   do j=0,N-1
     sample(j) = cmplx(cos(2.0 * PI_D * j/N), sin(2.0 * PI_D * j /N))
   end do

   call dfftw_execute(plan_fou)
   output = output/N

   write(*,*) 'output of the fourier transform'
   do j=0,N-1
     write(*,*) real(output(j)), aimag(output(j))
   end do

   do j=0,N
     ! initialize the sample to x^2 + x
     cheb_sample(j) = cos(j*PI/N) * (cos(j*PI/N) + 1.0)
   end do

   call dfftw_execute(plan_cheb)
   cheb_output(1:N-1) = cheb_output(1:N-1)/N
   cheb_output(0) = cheb_output(0)/(2.0 * N)
   cheb_output(N) = cheb_output(N)/(2.0 * N)

   write(*,*) 'output of the cosine/chebyshev transform'
   do j=0,N
     write(*,*) cheb_output(j)
   end do

   ! conclusions:
   ! for 1D fourier dft using double precision, the ratio of the largest
to the
   ! smallest coefficients is of the order of 1e-16
   ! for 1D cosine transform the ratio of the largest to the smallest is
of the
   ! order of 1e-8
end program end_of_spectrum

Relevant Pages

  • FFTW: successive real2complex 1-D FFTs fail
    ... I've been experimenting lately with the great FFTW library in C, ... Perform 1-D FFTs on every row of the same image array ... I first create the plans, then I load the image into fftw_in, and then I ... The strange thing is that the second transform produces wrong results. ...
    (comp.dsp)
  • Re: FFTW versus NRC FFT
    ... Before exchanging the FFT routine known from the NRC ... Applying both codes for that function the FFTW algorithm ... what the FFT NRC code does. ... In particular, if I remember correctly, the "sinft" transform that NRC ...
    (sci.math.num-analysis)
  • Re: FFTW Fortran Real2Complex R2C output array structure? fftw_plan_dft_r2c_3d?
    ... So, in the C FFTW ... for an r2c transform the last dimension is cut ...
    (comp.dsp)
  • Re: Confused by FFTW output
    ... > I am trying to get acquainted with the latest FFTW and wrote a small ... the FFT computes the discrete Fourier ... transform, which is *not the same* as the continuous Fourier ... (This periodicity is also why you ...
    (sci.math.num-analysis)
  • Fourier *sine* and *cosine* integral
    ... while deriving sine and cosine transform from a Fourier ... we will assume that fis odd function for ...
    (sci.math)