inverse fft, производящий закрытый, но неправильный вывод

Question

Я выполняю обратный фрейм-2 diff. Я использую свойства сопряжения и масштабирования, чтобы вернуть результат. Я сопрягаю свой входной вектор, выполняю регулярный radix-2 fft (не ifft), сопрягаю результаты, затем масштабируюсь на 1.0/N. Тем не менее, я не получаю правильные результаты:

int main(){ 
    const int n = 4; 
    complex<double> x[n]; 

     // Test signal 

    x[0] = complex<double>(10,0); 
    x[1] = complex<double>(-2,0); 
    x[2] = complex<double>(-2,2); 
    x[3] = complex<double>(-2,-2); 

     print(x,n); 

    fft_inverse(x,n); 

    print(x,n); 

} 
//dif fft. works 
void fft(complex<double> X[], int N){ 
     if(N == 1){return;} 

    complex<double> *temp = new complex<double>[N]; 
    for(int i=0; i<N; i++){ 
     temp[i]=X[i]; 
     } 
     for(int i = 0; i<N/2; i++){ 
     complex<double> tw(cos(-2*M_PI*i/N),sin(-2*M_PI*i/N)); 
      X[i] = temp[i] + temp[i+N/2]; 
      X[i+N/2] = temp[i]-temp[i+N/2]; 
      X[i+N/2] = X[i+N/2]*tw; 
     } 

     fft(X,N/2); 
     fft(X+N/2,N/2); 
} 
void fft_inverse(complex<double> X[], int N){ 
    //conjugate 
    for(int i = 0; i<N;i++){ 
     X[i] = conj(X[i]); 
    } 
    //perform fft 
    fft(X,N); 
    //conjugate again 
    for(int i = 0; i<N;i++){ 
     X[i] = conj(X[i]); 
    } 
    //scale by 1.0/N 
    double norm_N = 1.0/N; 
    for(int i = 0; i<N;i++){ 
     X[i] *= norm_N; 
    } 
} 

Вот мои результаты: ввода:

(10,0) (-2,0) (-2,2) (-2, -2)

Выход:

(1, -0) (3,1) (2,5, -0,5) (3,5, -0,5)

Вывод должен быть:

(1,0) (2,0) (3,0) (4,0)

Что происходит? Я проверил свой fft на том, что должно быть на выходе, и получил правильный результат, поэтому я не уверен, в чем проблема.

Answer 1

Похоже, что ваш код дает правильный вывод, но бункеры в неправильном порядке:

octave> X = [ 10, -2, -2+2i, -2-2i ] 
X = 

    10 + 0i -2 + 0i -2 + 2i -2 - 2i 

octave> x = ifft(X) 
x = 

    1.00000 + 0.00000i 2.50000 - 0.50000i 3.00000 + 1.00000i 3.50000 - 0.50000i 

fft_inverse() выглядит нормально, так что я подозреваю, что fft() может использовать некоторые дополнительные тестирования/отладки - это может быть, что вы необходимо сделать некоторую разворот бит для вашей индексации.

Answer 2

Перейдем к рассмотрению основ: Обратный fft длины 2N оценивает полином p(z) степени 2N у корней единства q^i где q^N = -1. Чтобы сделать это быстро, полином разбивается на четные и нечетные части, p(z)=pe(z^2)+z*po(z^2). Затем

p(q^i) =pe(q^2i)+q^i*po(q^2i) 
p(q^(N+i))=pe(q^2i)-q^i*po(q^2i) 

В реализациях, значение ре и ро получается рекурсивными вызовами ifft длиной N.

Передняя FFT теперь обратная операция, которая, вычисление коэффициентов полинома от его значения. Для этого, значения p(q^i) преобразуются в значения pe(q^2i) и po(q^2i) путем изменения приведенных выше формул,

pe(q^2i)=0.5*(p(q^i)+p(q^(N+i))) 
po(q^2i)=0.5*(p(q^i)-p(q^(N+i)))*q^(-i) 

который вы должны признать из вашего кода. Тогда fft рекурсивно вызывается для определения коэффициентов pe и po из их теперь вычисленных значений. Факторы 0,5 обычно опускаются и применяются как деление на длину массива вне подпрограмм fft.

Если вы ошиблись сейчас, то в том, что коэффициенты pe и po чередуются, чередуются, в последовательности коэффициентов p. Знаменитый рисунок бабочки.

void fft(complex<double> X[], int N){ 
    fft_internal(X,N,1); 
} 

void fft_internal(complex<double> X[], int N, int step){ 
    if(N == 1){return;} 

    complex<double> *temp = new complex<double>[N]; 
    for(int i=0; i<N; i++){ 
     temp[i]=X[step*i]; 
    } 
    for(int i = 0; i<N/2; i++){ 
    complex<double> tw(cos(-2*M_PI*i/N),sin(-2*M_PI*i/N)); 
     X[2*step*i  ] = temp[i] + temp[i+N/2]; 
     X[2*step*i+step] = (temp[i] - temp[i+N/2])*tw; 
    } 

    fft_internal(X  ,N/2, 2*step); 
    fft_internal(X+step,N/2, 2*step); 
} 

С правильно переупорядоченных входов ([2] < -> [1]),

// Test signal 

x[0] = complex<double>(10,0); 
x[1] = complex<double>(-2,2); 
x[2] = complex<double>(-2,0); 
x[3] = complex<double>(-2,-2); 

это возвращает ожидаемый результат.