Ошибки округления, дающие неверные результаты в ДПФ?

Question

Я бил головой о стену на этом ДПФ. Он должен распечатать: 8,0,0,0,0,0,0,0, но вместо этого я получаю 8, а затем очень очень маленькие цифры. Являются ли эти ошибки округления? Я могу что-нибудь сделать? Мой Radix2 FFT дает правильные результаты, кажется глупым, что DFT не может работать.

Я начал с сложных чисел, поэтому я знаю, что есть недостаток, я попытался разбить его, чтобы проиллюстрировать проблему.

#include <cstdlib> 
#include <math.h> 
#include <iostream> 
#include <complex> 
#include <cassert> 

#define SIZE 8 
#define M_PI 3.14159265358979323846 

void fft(const double src[], double dst[], const unsigned int n) 
{ 
    for(int i=0; i < SIZE; i++) 
    { 
     const double ph = -(2*M_PI)/n; 
     const int gid = i; 

     double res = 0.0f; 
     for (int k = 0; k < n; k++) { 

      double t = src[k]; 

      const double val = ph * k * gid; 
      double cs = cos(val); 
      double sn = sin(val); 

      res += ((t * cs) - (t * sn)); 
      int a = 1; 
     } 

     dst[i] = res; 
     std::cout << dst[i] << std::endl; 
    } 
} 

int main(void) 
{ 
    double array1[SIZE]; 
    double array2[SIZE]; 

    for(int i=0; i < SIZE; i++){ 
     array1[i] = 1; 
     array2[i] = 0; 
    } 

    fft(array1, array2, SIZE); 

    return 666; 
} 

Answer 1

БПФ может фактически получить более точные результаты, чем прямой DFT расчета, так как меньше арифметические опсы, как правило, позволяют меньше возможностей для ошибок арифметики квантования аккумулировать. В этом разделе есть статья одного из авторов FFTW.

Поскольку DFT/FFT имеет дело с трансцендентальной базой, результаты никогда не будут (за исключением, возможно, в нескольких особых случаях или случайной случайностью) быть точно верными, используя любой несимволический и конечный формат компьютерного номера. Таким образом, значения очень близкие (в пределах нескольких LSB) к нулю должны просто игнорироваться как шум или считаться одинаковыми с нулем.