algorithm
Быстрое преобразование Фурье
Поиск…
Вступление
Вещественная и комплексная форма ДПФА (D iscrete F ourier T ransforms) может быть использована для выполнения частотного анализа или синтеза для любых дискретных и периодических сигналов. FFT ( F ast F ourier T ransform) представляет собой реализацию DFT, которая может быть быстро выполнена на современных процессорах.
Radix 2 FFT
Простейшим и, возможно, самым известным методом вычисления БПФ является алгоритм Decix-2 Decimation in Time. Radix-2 FFT работает, разлагая N-точечный сигнал во временной области в N сигналов временной области, каждый из которых состоит из одной точки ,
Разложение сигнала или «прореживание во времени» достигается путем бина, изменяющего индексы для массива данных временной области. Таким образом, для шестнадцатиточечного сигнала образец 1 (двоичный 0001) заменяется образцом 8 (1000), образец 2 (0010) заменяется на 4 (0100) и так далее. Образцовая свопинг с использованием метода обратного бита может быть достигнута просто в программном обеспечении, но ограничивает использование Radix 2 FFT сигналами длиной N = 2 ^ M.
Значение одноточечного сигнала во временной области равно его значению в частотной области, поэтому этот массив разложенных единичных точек временной области не требует преобразования, чтобы стать массивом точек частотной области. N отдельных точек; однако необходимо перестроить в один N-точечный частотный спектр. Оптимальная реконструкция полного частотного спектра выполняется с использованием расчетов бабочки. Каждый этап реконструкции в Fix Radix-2 выполняет множество двухточечных бабочек, используя аналогичный набор экспоненциальных весовых функций Wn ^ R.
БПФ удаляет избыточные вычисления в Дискретном преобразовании Фурье, используя периодичность Wn ^ R. Спектральная реконструкция завершена в log2 (N) этапах расчетов бабочек, дающих X [K]; реальные и мнимые данные частотной области в прямоугольной форме. Для преобразования в амплитуду и фазу (полярные координаты) требуется найти абсолютное значение, √ (Re2 + Im2) и аргумент tan-1 (Im / Re).
Ниже показана полная диаграмма потока бабочек для восьмиступенчатого Fix Radix 2. Обратите внимание, что входные сигналы ранее были переупорядочены в соответствии с процедурой прореживания во времени, описанной ранее.
БПФ обычно работает на сложных входах и производит сложный вывод. Для реальных сигналов мнимая часть может быть установлена равной нулю, а действительная часть установлена на входной сигнал x [n], однако возможны многие оптимизации, связанные с преобразованием данных реального времени. Значения Wn ^ R, используемые во время реконструкции, могут быть определены с использованием экспоненциального весового уравнения.
Значение R (экспоненциальная весовая мощность) определяется текущим этапом в спектральной реконструкции и вычислением тока в конкретной бабочке.
Пример кода (C / C ++)
Пример кода AC / C ++ для вычисления Radix 2 FFT можно найти ниже. Это простая реализация, которая работает для любого размера N, где N - мощность 2. Она примерно на 3 раза медленнее, чем самая быстрая реализация FFTw, но все же является хорошей основой для будущей оптимизации или для изучения того, как работает этот алгоритм.
#include <math.h>
#define PI 3.1415926535897932384626433832795 // PI for sine/cos calculations
#define TWOPI 6.283185307179586476925286766559 // 2*PI for sine/cos calculations
#define Deg2Rad 0.017453292519943295769236907684886 // Degrees to Radians factor
#define Rad2Deg 57.295779513082320876798154814105 // Radians to Degrees factor
#define log10_2 0.30102999566398119521373889472449 // Log10 of 2
#define log10_2_INV 3.3219280948873623478703194294948 // 1/Log10(2)
// complex variable structure (double precision)
struct complex
{
public:
double Re, Im; // Not so complicated after all
};
// Returns true if N is a power of 2
bool isPwrTwo(int N, int *M)
{
*M = (int)ceil(log10((double)N) * log10_2_INV);// M is number of stages to perform. 2^M = N
int NN = (int)pow(2.0, *M);
if ((NN != N) || (NN == 0)) // Check N is a power of 2.
return false;
return true;
}
void rad2FFT(int N, complex *x, complex *DFT)
{
int M = 0;
// Check if power of two. If not, exit
if (!isPwrTwo(N, &M))
throw "Rad2FFT(): N must be a power of 2 for Radix FFT";
// Integer Variables
int BSep; // BSep is memory spacing between butterflies
int BWidth; // BWidth is memory spacing of opposite ends of the butterfly
int P; // P is number of similar Wn's to be used in that stage
int j; // j is used in a loop to perform all calculations in each stage
int stage = 1; // stage is the stage number of the FFT. There are M stages in total (1 to M).
int HiIndex; // HiIndex is the index of the DFT array for the top value of each butterfly calc
unsigned int iaddr; // bitmask for bit reversal
int ii; // Integer bitfield for bit reversal (Decimation in Time)
int MM1 = M - 1;
unsigned int i;
int l;
unsigned int nMax = (unsigned int)N;
// Double Precision Variables
double TwoPi_N = TWOPI / (double)N; // constant to save computational time. = 2*PI / N
double TwoPi_NP;
// complex Variables (See 'struct complex')
complex WN; // Wn is the exponential weighting function in the form a + jb
complex TEMP; // TEMP is used to save computation in the butterfly calc
complex *pDFT = DFT; // Pointer to first elements in DFT array
complex *pLo; // Pointer for lo / hi value of butterfly calcs
complex *pHi;
complex *pX; // Pointer to x[n]
// Decimation In Time - x[n] sample sorting
for (i = 0; i < nMax; i++, DFT++)
{
pX = x + i; // Calculate current x[n] from base address *x and index i.
ii = 0; // Reset new address for DFT[n]
iaddr = i; // Copy i for manipulations
for (l = 0; l < M; l++) // Bit reverse i and store in ii...
{
if (iaddr & 0x01) // Detemine least significant bit
ii += (1 << (MM1 - l)); // Increment ii by 2^(M-1-l) if lsb was 1
iaddr >>= 1; // right shift iaddr to test next bit. Use logical operations for speed increase
if (!iaddr)
break;
}
DFT = pDFT + ii; // Calculate current DFT[n] from base address *pDFT and bit reversed index ii
DFT->Re = pX->Re; // Update the complex array with address sorted time domain signal x[n]
DFT->Im = pX->Im; // NB: Imaginary is always zero
}
// FFT Computation by butterfly calculation
for (stage = 1; stage <= M; stage++) // Loop for M stages, where 2^M = N
{
BSep = (int)(pow(2, stage)); // Separation between butterflies = 2^stage
P = N / BSep; // Similar Wn's in this stage = N/Bsep
BWidth = BSep / 2; // Butterfly width (spacing between opposite points) = Separation / 2.
TwoPi_NP = TwoPi_N*P;
for (j = 0; j < BWidth; j++) // Loop for j calculations per butterfly
{
if (j != 0) // Save on calculation if R = 0, as WN^0 = (1 + j0)
{
//WN.Re = cos(TwoPi_NP*j)
WN.Re = cos(TwoPi_N*P*j); // Calculate Wn (Real and Imaginary)
WN.Im = -sin(TwoPi_N*P*j);
}
for (HiIndex = j; HiIndex < N; HiIndex += BSep) // Loop for HiIndex Step BSep butterflies per stage
{
pHi = pDFT + HiIndex; // Point to higher value
pLo = pHi + BWidth; // Point to lower value (Note VC++ adjusts for spacing between elements)
if (j != 0) // If exponential power is not zero...
{
//CMult(pLo, &WN, &TEMP); // Perform complex multiplication of Lovalue with Wn
TEMP.Re = (pLo->Re * WN.Re) - (pLo->Im * WN.Im);
TEMP.Im = (pLo->Re * WN.Im) + (pLo->Im * WN.Re);
//CSub (pHi, &TEMP, pLo);
pLo->Re = pHi->Re - TEMP.Re; // Find new Lovalue (complex subtraction)
pLo->Im = pHi->Im - TEMP.Im;
//CAdd (pHi, &TEMP, pHi); // Find new Hivalue (complex addition)
pHi->Re = (pHi->Re + TEMP.Re);
pHi->Im = (pHi->Im + TEMP.Im);
}
else
{
TEMP.Re = pLo->Re;
TEMP.Im = pLo->Im;
//CSub (pHi, &TEMP, pLo);
pLo->Re = pHi->Re - TEMP.Re; // Find new Lovalue (complex subtraction)
pLo->Im = pHi->Im - TEMP.Im;
//CAdd (pHi, &TEMP, pHi); // Find new Hivalue (complex addition)
pHi->Re = (pHi->Re + TEMP.Re);
pHi->Im = (pHi->Im + TEMP.Im);
}
}
}
}
pLo = 0; // Null all pointers
pHi = 0;
pDFT = 0;
DFT = 0;
pX = 0;
}
Радикс 2 Обратный БПФ
Из-за сильной двойственности преобразования Фурье, настройка выхода прямого преобразования может привести к обратному БПФ. Данные в частотной области могут быть преобразованы во временную область следующим образом:
- Найти комплексное сопряжение данных частотной области путем инвертирования мнимой составляющей для всех экземпляров K.
- Выполнение прямого БПФ на данных сопряженной частотной области.
- Разделите каждый вывод результата этого FFT на N, чтобы указать значение истинной временной области.
- Найти комплексное сопряжение вывода путем инвертирования мнимой составляющей данных во временной области для всех экземпляров n.
Примечание : данные частоты и временной области являются комплексными переменными. Обычно мнимая составляющая сигнала временной области после обратного БПФ либо равна нулю, либо игнорируется как ошибка округления. Повышение точности переменных с 32-битного поплавка до 64-битного двойного или 128-битного удвоения значительно снижает ошибки округления, возникающие в результате нескольких последовательных операций БПФ.
Пример кода (C / C ++)
#include <math.h>
#define PI 3.1415926535897932384626433832795 // PI for sine/cos calculations
#define TWOPI 6.283185307179586476925286766559 // 2*PI for sine/cos calculations
#define Deg2Rad 0.017453292519943295769236907684886 // Degrees to Radians factor
#define Rad2Deg 57.295779513082320876798154814105 // Radians to Degrees factor
#define log10_2 0.30102999566398119521373889472449 // Log10 of 2
#define log10_2_INV 3.3219280948873623478703194294948 // 1/Log10(2)
// complex variable structure (double precision)
struct complex
{
public:
double Re, Im; // Not so complicated after all
};
void rad2InverseFFT(int N, complex *x, complex *DFT)
{
// M is number of stages to perform. 2^M = N
double Mx = (log10((double)N) / log10((double)2));
int a = (int)(ceil(pow(2.0, Mx)));
int status = 0;
if (a != N) // Check N is a power of 2
{
x = 0;
DFT = 0;
throw "rad2InverseFFT(): N must be a power of 2 for Radix 2 Inverse FFT";
}
complex *pDFT = DFT; // Reset vector for DFT pointers
complex *pX = x; // Reset vector for x[n] pointer
double NN = 1 / (double)N; // Scaling factor for the inverse FFT
for (int i = 0; i < N; i++, DFT++)
DFT->Im *= -1; // Find the complex conjugate of the Frequency Spectrum
DFT = pDFT; // Reset Freq Domain Pointer
rad2FFT(N, DFT, x); // Calculate the forward FFT with variables switched (time & freq)
int i;
complex* x;
for ( i = 0, x = pX; i < N; i++, x++){
x->Re *= NN; // Divide time domain by N for correct amplitude scaling
x->Im *= -1; // Change the sign of ImX
}
}