FFT и обратный FFT

Question

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

short FFT(short int dir,long m,double *x,double *y) 
{ 
    long n,i,i1,j,k,i2,l,l1,l2; 
    double c1,c2,tx,ty,t1,t2,u1,u2,z; 

    /* Calculate the number of points */ 
    n = 1; 
    for (i=0;i<m;i++) 
     n *= 2; 

    /* Do the bit reversal */ 
    i2 = n >> 1; 
    j = 0; 
    for (i=0;i<n-1;i++) { 
     if (i < j) { 
     tx = x[i]; 
     ty = y[i]; 
     x[i] = x[j]; 
     y[i] = y[j]; 
     x[j] = tx; 
     y[j] = ty; 
     } 
     k = i2; 
     while (k <= j) { 
     j -= k; 
     k >>= 1; 
     } 
     j += k; 
    } 

    /* Compute the FFT */ 
    c1 = -1.0; 
    c2 = 0.0; 
    l2 = 1; 
    for (l=0;l<m;l++) { 
     l1 = l2; 
     l2 <<= 1; 
     u1 = 1.0; 
     u2 = 0.0; 
     for (j=0;j<l1;j++) { 
     for (i=j;i<n;i+=l2) { 
      i1 = i + l1; 
      t1 = u1 * x[i1] - u2 * y[i1]; 
      t2 = u1 * y[i1] + u2 * x[i1]; 
      x[i1] = x[i] - t1; 
      y[i1] = y[i] - t2; 
      x[i] += t1; 
      y[i] += t2; 
     } 
     z = u1 * c1 - u2 * c2; 
     u2 = u1 * c2 + u2 * c1; 
     u1 = z; 
     } 
     c2 = sqrt((1.0 - c1)/2.0); 
     if (dir == 1) 
     c2 = -c2; 
     c1 = sqrt((1.0 + c1)/2.0); 
    } 

    /* Scaling for forward transform */ 
    if (dir == 1) { 
     for (i=0;i<n;i++) { 
     x[i] /= n; 
     y[i] /= n; 
     } 
    } 

    return(true); 
} 

Но этот код, просто поддерживают массивы размера 2^m .like CLRS книги коды.
Наши массивы, которые должны быть преобразованы, не входят в этот диапазон, и добавление нуля будет дорогостоящим, поэтому я ищу другое решение, которое поможет мне вводить любые размеры.
Что-то вроде того, что IT++ и matlab сделать. но, как мы хотим, чтобы это в чистом C, используя их в impossible.Moreover, IT++ код был заблокирован, как я проверил

Answer 1

Если вы работаете на любой вычислительной платформе массового рынка (Intel с Windows или OS X, iOS и т. Д.), То существует высокопроизводительная реализация FFT, поставляемая поставщиком или производителем.

В противном случае вы должны оценить FFTW.

Написание высокопроизводительных БПФ для размеров, отличных от полномочий двух, является сложной задачей.

Если вы собираетесь использовать свою собственную реализацию, то, что касается только размеров мощности из-два:

Реализация вы показываете вычисляет sqrt во время быстрого преобразования Фурье. Большинство высокопроизводительных реализаций FFT вычисляют константы заблаговременно и сохраняют их в таблице.

Масштабирование содержит операции деления, в x[i] /= n и y[i] /= n. Компилятор, скорее всего, реализует их как инструкции разделения. Подразделение обычно представляет собой медленную инструкцию для обычных процессоров. Было бы лучше вычислить scale = 1./n один раз и умножить на scale вместо деления на n.

Еще лучше было бы опустить масштаб целиком. Преобразование часто полезно без шкалы, или масштаб можно опустить из отдельных преобразований и применять только один раз в виде агрегированной шкалы. (Например., вместо выполнения двух масштабных операций, один в прямом БПФ и один в обратном БПФ, оставьте операцию масштабирования из подпрограмм БПФ и сделайте это только один раз после прямого БПФ и обратного БПФ.)

Вы может быть в состоянии опустить перестановку бит-разворота, если приемлемо иметь данные частотной области в порядке бит-обратного.

Если вы сохраняете перестановку битов, ее можно оптимизировать. Техника для этого зависит от платформы. На некоторых платформах имеется инструкция для изменения битов в целое число (например, ARM имеет rbit). Если вашей платформы нет, вы можете захотеть сохранить индексы разворота бит в таблице или исследовать способы их вычисления быстрее, чем ваш текущий код.

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

В вашем текущем коде используется бабочка radix-2. Radix-4 обычно лучше, потому что он получает некоторое преимущество от того факта, что умножение на i можно сделать только путем изменения того, какая часть данных используется, и изменения некоторых дополнений к вычитаниям и наоборот.

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

Если у вашего целевого процессора есть функции SIMD, вы должны абсолютно использовать их в БПФ; они значительно ускоряют работу FFT.

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

Answer 2

В вашей реализации, я обеспокоен этот бит кода:

 z = u1 * c1 - u2 * c2; 
    u2 = u1 * c2 + u2 * c1; 
    u1 = z; 

(u1, u2) подлежит большому количеству накопленных ошибок округления, когда l1 большой. Вероятно, вы получите неточное преобразование.

Я бы откликнулся на рекомендацию FFTW, но я считаю, что это очень специфично для платформы. (Большинство библиотек FFT есть.) [EDIT: Нет, это на самом деле прямой C89. Это именно то, что вы говорите, что ищете.]

На странице Wikipedia FFT перечислены кучи алгоритмов, которые работают на входные массивы странного размера. Я не знаю, как они работают, но я считаю, что общая идея заключается в том, что вы используете Rader's algorithm или Bluestein's algorithm для входов с простым размером и Cooley-Tukey для уменьшения преобразований композитного размера в кучу простых преобразований.

Answer 3

Для альтернативы FFTW проверьте мой БРП-микс-радиус, который также обрабатывает длины FFT без 2^N. C-источник доступен от http://www.corix.dk/Mix-FFT/mix-fft.html.