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