Проблема с моей реализацией FFT в C

Стандартный нерекурсивный алгоритм radix-2. Я понимаю, что у меня здесь нет оптимального для производительности (например, повторные триггерные звонки), я просто хочу, чтобы все было правильно, прежде чем оптимизировать и портировать VHDL. Я использую заголовок complex.h.Проблема с моей реализацией FFT в C

Я получаю правильные результаты для массивов длины 8, но не для 16 или 32. ОДНАКО, для реальных сигналов мы все равно получаем правильный результат сопряженной симметрии, только значения ошибочны. Это заставляет меня думать, что сортировка правильная, но, может быть, что-то не так с триггерными вызовами для меньших значений? Может ли кто-нибудь помочь?

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

void FFT(double complex x[], int N){ 

    //bit reverse 

int s = log2(N); 

for(int i = 0; i < N/2; i++){ 
    int h = bitrev(i, s); 

    printf("%u ", h); 
    double complex temp = x[i]; 
    x[i] = x[h]; 
    x[h] = temp; 
} 

unsigned int Np = 2; //NUM POINTS IN EACH BLOCK. INITIALLY 2 
unsigned int Bp = (N/2); 
for(int i = 0; i < s; i++){ 
    int NpP = Np>>1; //num butterflies 
    int BaseT = 0; 
    for (int j = 0; j < Bp; j++){ 

     int BaseB = BaseT + NpP; 

     for(int k = 0; k < NpP; k++){ 

      double complex top = x[BaseT + k]; 
      double complex bot = (ccos(2*pi*k/Np) - I*csin(2*pi*k/Np))*x[BaseB + k]; 
      x[BaseT + k] = top+bot; 
      x[BaseB + k] = top-bot; 
     } 
     BaseT = BaseT + Np; 
    } 
    Bp = Bp>>1; 
    Np = Np<<1; 
}

}

Выход печатается с:

for(int i = 0; i < LENGTH; i++){ 
    printf("%f + %fj\n", creal(x[i]), cimag(x[i])); 
}

Вот ввода длины 8, а правильный выход:

double complex x[LENGTH] = {1.0, -1.0, 1.0, -1.0, 10.0, -1.0, 5.0, -1.0};

Выход:

13.000000 + 0.000000j 
-9.000000 + 4.000000j 
5.000000 + 0.000000j 
-9.000000 + -4.000000j 
21.000000 + 0.000000j 
-9.000000 + 4.000000j 
5.000000 + 0.000000j 
-9.000000 + -4.000000j

Вот длина 16 вход, неправильный вывод, который я получаю, и что это должно быть:

double complex x[SIZE] = {10.0, -1.0, 5.0, -1.0, 4.0, 4.0, 2.0, 6.0, 9.0, 3.0, 8.0, 4.0, 5.0, 4.0, 3.0, 8.0};

выход (неверно):

73.000000 + 0.000000j 
-4.882497 + 10.451032j 
12.535534 + 5.020815j 
6.975351 + 7.165394j 
12.000000 + 7.000000j 
-0.732710 + 0.094326j 
5.464466 + 19.020815j 
2.639856 + 3.379964j 
19.000000 + 0.000000j 
2.639856 + -3.379964j 
5.464466 + -19.020815j 
-0.732710 + -0.094326j 
12.000000 + -7.000000j 
6.975351 + -7.165394j 
12.535534 + -5.020815j 
-4.882497 + -10.451032j

Правильный выход:

73.000000 + 0.000000j 
-4.175390 + 10.743925j 
13.535534 + 4.020815j 
6.268244 + 5.458287j 
10.000000 + 7.000000j 
-1.439817 + 1.801433j 
6.464466 + 20.020815j 
3.346963 + 3.087071j 
19.000000 + 0.000000j 
3.346963 + -3.087071j 
6.464466 + -20.020815j 
-1.439817 + -1.801433j 
10.000000 + -7.000000j 
6.268244 + -5.458287j 
13.535534 + -4.020815j 
-4.175390 + -10.743925j

Для этой длины 16, например, он дает правильный выход (!?):

double complex x[SIZE] = {30.0, -1.0, 4.0, -6.0, 4.0, 9.0, 2.0, 6.0, 9.0, 3.0, -8.0, 4.0, -5.0, 4.0, 3.0, 8.0}; 

66.000000 + 0.000000j 
22.604378 + -9.862676j 
43.535534 + 28.091883j 
24.900438 + 4.851685j 
37.000000 + -3.000000j 
-1.285214 + 2.408035j 
36.464466 + 10.091883j 
37.780399 + 23.693673j 
12.000000 + 0.000000j 
37.780399 + -23.693673j 
36.464466 + -10.091883j 
-1.285214 + -2.408035j 
37.000000 + 3.000000j 
24.900438 + -4.851685j 
43.535534 + -28.091883j 
22.604378 + 9.862676j

EDIT: Я осознал свою проблему: мой цикл обратного хода был совершенно неправильным. Я должен был подумать об этой, казалось бы, простой части, и проверить ее более тщательно. Причина, по которой она работала для последнего случая, состояла в том, что она имела повторяющиеся значения, которые по совпадению дали тот же вектор, что и рабочий поворот. Вот неподвижная петля:

int s = log2(N); 

    for(int i = 0; i < N; i++){ 
     int h = bitrev(i, s); 
     if(i < h){ 
     double complex temp = x[i]; 
     x[i] = x[h]; 
     x[h] = temp; 
     } 
    }

источник

2014-09-29 Matthew Hampsey

I» m не знакомы с кодом FFT. Можете ли вы предоставить примерный набор данных размером 16 или 32, а также полученный результат (и код, который печатает этот ответ), и укажите, какой ответ вы считаете нужным. Это может помочь. –

@JonathanLeffler Я добавил запрошенные объекты –

wild guess: int overflow? –

Эта часть 2*pi*k/Np выпускающая компилятор решить отливку типа. Я бы также проверил, что эти значения вычисляются должным образом. Например, я получаю разные результаты при использовании 2*pi*k/Np; и 2.0*pi*(double) k*(double) Np;

(Модификации: Я скорее добавить это как комментарий, но у меня нет достаточно репутации Не стесняйтесь, чтобы переместить мой ответ :).)

источник

2014-09-29 21:21:29 3150

Спасибо, я изменился на явные приведения. Я понял, какова моя фактическая проблема: мой цикл разворота бит был неправильным. Я редактировал правильный код. –

Проблема с моей реализацией FFT в C

ответ

Смежные вопросы