Как сравнить реальную реализацию преобразования Фурье из числовых рецептов в Matlab fft?

Question

Я вычисляю дискретное преобразование Фурье с использованием числовых рецептов и (для подтверждения результата) Matlab. Работа с реальными значениями.

Мой Matlab код

in(1)=0.0; 
in(2)=1.0; 
in(3)=2.0; 
in(4)=3.0; 
in(5)=4.0; 
in(6)=5.0; 
in(7)=6.0; 
in(8)=7.0; 

out = fft(in); 

дает мне

out = 
28.0000 + 0.0000i 
-4.0000 + 9.6569i 
-4.0000 + 4.0000i 
-4.0000 + 1.6569i 
-4.0000 + 0.0000i 
-4.0000 - 1.6569i 
-4.0000 - 4.0000i 
-4.0000 - 9.6569i 

Что data вход мне нужно отправить Числовые рецепты

для того, чтобы получить тот же результат, как Matlab ?

С этой NR forum ссылки, я обнаружил, что вход realft должно быть сдвинут на одно место, поэтому я использую in[0] = 0, и, следовательно, вход 1 элемент больше, чем N.

Код испытания:

#include <stdio.h> 
#define LEN 8 
int main() 
{ 
    float inout[LEN+1]; 

    inout[0] = 0.0; 
    inout[1] = 0.0; 
    inout[2] = 1.0; 
    inout[3] = 2.0; 
    inout[4] = 3.0; 
    inout[5] = 4.0; 
    inout[6] = 5.0; 
    inout[7] = 6.0; 
    inout[8] = 7.0; 

    realft(inout, LEN, 1); 

    for(unsigned int i=0; i<LEN+1; i=i+1) 
     printf("%15.10f \n",inout[i]); 

    return 0; 
} 

Выход из тестового кода:

0.00000000 
28.00000000 
-4.00000000 
-4.00000000 
-9.65685463 
-4.00000000 
-4.00000000 
-4.00000000 
-1.65685427 

который похож, но отличается от Matlab.

realft берется из Numerical Recipes (использует four1 рутина):

Answer 1

Первый аспект, который отличается от реализаций MatLab и Numerical Recipes', кроме индексации вне по-одному, что вы уже заметили, и что они основаны на несколько другом определении БПФ. Более конкретно, Matlab uses a negative complex exponential for the forward transform, в то время как числовые рецепты используют положительную комплексную экспоненту. Соответственно, реализация числовых рецептов даст результат, который является комплексно сопряженным с результатом из Matlab.

Другое дело, что эти реализации приводят к результату в другом порядке упаковки, а реализация Numerical Recipes выводит только не избыточную нижнюю половину спектра. Графически это отображение может быть представлено (после изменения знака мнимой части в соответствии с предыдущим пунктом относительно комплексного сопряжения) следующим графиком:

Answer 2

Я думаю, вам нужно обратиться к исходный код realft(). Кажется, что в использовании этой функции что-то не так. Причина в том, что хотя вход находится в реальных значениях, значения спектра, как правило, сложны. Однако в вашей программе результаты возвращаются из исходного inout, который находится в реальном типе данных, и мнимая часть отбрасывается.

Answer 3

Для двух примеров у вас разные входы (аналогичные, но разные). Для примера Matlab ваш вход был рампой от 0-> 7 включительно. Но из-за различий в индексировании (которые кто-то еще упоминал) в реализации C ваш ввод был 0,0,1,2,3,4,5,6 (в отличие от 0,1,2,3,4 , 5,6,7). Если вы выполняете связку matlab с C, разница в индексировании массива (matlab начинается с 1, C начинается с 0).