-1

this is the first time I'm asking a question here, so i hope you'll understand my problem.

The Thing is, that I want to do my own fft(), without using the given one in R. So far it works well for a series like seq(1,5).

But for c(1,1) something strange happen. As far as I could point it out it seems that x - x is not 0 in that case. Here the lines of code:

    series <- c(1,1)                  # defining the Serie
        nr_samples <- length(series)      # getting the length

#################
# Calculating the harmonic frequncy
#################

        harmonic <- seq(0,(nr_samples-1))
        harmonic <- 2*pi*harmonic
        harmonic <- harmonic/nr_samples

#################
# Exponential funktion needed for summing up
#################

        exponential <- function(index, omega){

          result <- exp(-((0+1i)*omega*index))

          return(result)

        }

#################
# The sum for calculating the fft
#################

        my_fft <- function(time_series, omega){

          nr_samples <- length(time_series)
          summand <- 0

   # In the next loop the mistakes Happens       
   # While running this loop for harmonic[2]
   # the rseult should be 0 because
   # summand = - exp_factor
   # The result for summand + exp_factor 
   # is 0-1.22464679914735e-16i

     for (i in 1:nr_samples){

            exp_factor <- exponential((i-1), omega)            
            summand <- summand + time_series[i]*exp_factor     
            print(paste("Summand", summand, "Exp", exp_factor))
          }                                                    

          return(summand)                                      
        }


    transform <- sapply(harmonic, function(x){my_fft(series,x)})
    fft_transform <- fft(series)
    df <- data.frame(transform, fft_transform)
    print(df)

Could anyone tell me, why summand + exp_factor, for harmonic[2] is not zero??

Thomas
  • 43,637
  • 12
  • 109
  • 140
Benjamin Mohn
  • 301
  • 3
  • 12

1 Answers1

3

This is commonly referred to FAQ 7.31 which says:

The only numbers that can be represented exactly in R’s numeric type are integers and fractions whose denominator is a power of 2. Other numbers have to be rounded to (typically) 53 binary digits accuracy. As a result, two floating point numbers will not reliably be equal unless they have been computed by the same algorithm, and not always even then. For example

R> a <- sqrt(2)
R> a * a == 2
[1] FALSE
R> a * a - 2
[1] 4.440892e-16

The function all.equal() compares two objects using a numeric tolerance of .Machine$double.eps ^ 0.5. If you want much greater accuracy than this you will need to consider error propagation carefully.

For more information, see e.g. David Goldberg (1991), “What Every Computer Scientist Should Know About Floating-Point Arithmetic”, ACM Computing Surveys, 23/1, 5–48, also available via http://www.validlab.com/goldberg/paper.pdf.

To quote from “The Elements of Programming Style” by Kernighan and Plauger:

10.0 times 0.1 is hardly ever 1.0.

(End of quote)

The Goldberg paper is legendary, and you may want to read it. This is a property of all floating point computation and not specific to R.

Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725