Ostatnio robię sobie program, który ma obrazować fragment szerokiego widma sygnału z duża częstotliwością próbkowania.
Na przykład, chciałbym wyświetlić któreś 500 wartości 65536-punktowego widma (np. 500 pierwszych, 500 ostatnich lub 500 gdzieś ze środka). Ogólnie mówiąc, chodzi o obliczenie i zobrazowanie maksymalnie 5% widma o bardzo dużej liczbie punktów.
W tym celu zastosowałem algorytm FFT wzięty z https://www.nayuki.io/page/free-small-fft-in-multiple-languages (wersja "C (C99 and above") i własny algorytm DFT wprost z definicji, ale z możliwością obliczenia tylko fragmentu widma zamiast całego. Oba algorytmy bazują na tablicy float[] i podzieliłem na część inicjalizującą (część, którą wystarczy wykonać tylko raz w przypadku wielokrotnego obliczania DFT o tych samych parametrach) oraz wykonawczą (obliczanie widma konkretnego sygnału).
Okazuje się, że czasowo algorytm FFT obliczający całe widmo działa wyraźnie szybciej od DFT obliczającego fragment. Dlaczego tak się dzieje nawet wtedy, gdy chcę obliczyć 1% widma, wtedy złożoność staje się bardziej zbliżona do O(n) zamiast do O(n^2), a ściślej złożoność wynosi O(n*m), gdzie n to szerokość widma, a m to liczba wartości widma w potrzebnym zakresie (DFT trwa 2-3 razy dłużej niż FFT)?
W jaki sposób w algorytmie FFT można wymusić pominięcie obliczania części widma (chociażby obliczenie połowy lub ćwiartki widma zamiast całego)? W obu algorytmach, w funkcji inicjującej należy podać wielkość widma, jako parametr (wtedy do obliczenia potrzeba podać tablicę próbek dwa razy większą).
Zamieszczam kod DFT bez obliczania argumentu, bo potrzebuję tylko moduł z wartości transformaty.
#ifndef DFT_F_H
#define DFT_F_H
#include <cmath>
#define SignalValueType_F float
#define M_PI 3.1415926535897932384626433832795
class DFT_F
{
public:
DFT_F();
~DFT_F();
SignalValueType_F * FourierC;
SignalValueType_F * FourierS;
void Init(int N);
void Fourier(SignalValueType_F * SignalRAW, SignalValueType_F * SignalA, int RAWCount, int TransCount, int TransStep, int TransOffset);
};
#endif // DFT_F_H
#include "dft_f.h"
DFT_F::DFT_F()
{
FourierC = new SignalValueType_F[2100000];
FourierS = new SignalValueType_F[2100000];
}
DFT_F::~DFT_F()
{
delete[] FourierC;
delete[] FourierS;
}
///
/// \brief DFT_F::Init - Tworzenie tablic wartosci funkcji sinus i cosinus
/// \param N
///
void DFT_F::Init(int N)
{
int i;
float iX;
float NX;
NX = N;
for (i = 0; i < N; i++)
{
iX = i;
FourierC[i] = cosf(2.0 * M_PI * iX / NX);
FourierS[i] = sinf(2.0 * M_PI * iX / NX);
}
}
///
/// \brief DFT_F::Fourier
/// \param SignalRAW - Tablica probek sygnalu surowego
/// \param SignalA - Tablica amplitudy widma
/// \param RAWCount - Sygnal surowy - liczba probek
/// \param TransCount - Widmo - Liczba wartosci
/// \param TransStep - Widmo - Krok obliczania wartosci (decymacja)
/// \param TransOffset - Widmo - Wartosc poczatkowa
///
void DFT_F::Fourier(SignalValueType_F * SignalRAW, SignalValueType_F * SignalA, int RAWCount, int TransCount, int TransStep, int TransOffset)
{
int i, ii;
long iX, LUT;
float SignalRAWC;
float SignalRAWS;
float RAWCountX = RAWCount;
iX = TransOffset % RAWCount;
for (i = 0; i < TransCount; i++)
{
SignalRAWC = 0;
SignalRAWS = 0;
LUT = iX;
for (ii = 0; ii < RAWCount; ii++)
{
while (LUT >= RAWCount)
{
LUT -= RAWCount;
}
SignalRAWC += (SignalRAW[ii] * FourierC[LUT]);
SignalRAWS += (SignalRAW[ii] * FourierS[LUT]);
LUT = LUT + iX;
}
SignalRAWC = (SignalRAWC / RAWCountX);
SignalRAWS = (SignalRAWS / RAWCountX);
SignalA[i] = sqrtf((SignalRAWC * SignalRAWC) + (SignalRAWS * SignalRAWS));
iX += TransStep;
}
}
Algorytm FFT po odpowiednim dostosowaniu. Funkcja Init inicjalizuje obiekt, FFT oblicza transformatę, a IFFT oblicza transformatę odwrotną. W tym przypadku, należy podać dwie tablice, jedna to część rzeczywista, a druga to część urojona. Po przeliczeniu można otrzymać moduł i argument za pomocą funkcji Magnitude i Phase.
#ifndef FFT_H
#define FFT_H
#include <cmath>
#define FloatingType float
class FFTX
{
public:
FFTX();
~FFTX();
void Init(size_t FourierBase, int FourierWindow, int WinFactor);
void FFT(FloatingType real[], FloatingType imag[]);
void FFTW(FloatingType real[], FloatingType imag[]);
void IFFT(FloatingType real[], FloatingType imag[]);
FloatingType Magnitude(FloatingType real, FloatingType imag);
FloatingType Phase(FloatingType real, FloatingType imag);
private:
int FourierWindow_;
size_t FourierBase_;
FloatingType * cos_table;
FloatingType * sin_table;
unsigned int Fourier_levels;
void transform_radix2_init(size_t n);
void transform_radix2(FloatingType real[], FloatingType imag[], size_t n);
};
#endif // FFT_H
#include "fftx.h"
FFTX::FFTX()
{
cos_table = new FloatingType[4000000];
sin_table = new FloatingType[4000000];
}
FFTX::~FFTX()
{
delete[] cos_table;
delete[] sin_table;
}
#define FloatingType float
void FFTX::Init(size_t FourierBase)
{
FloatingType M_PI = 3.1415926535897932384626433832795;
size_t I;
FourierBase_ = FourierBase;
transform_radix2_init(FourierBase);
}
void FFTX::transform_radix2_init(size_t n)
{
size_t i;
FloatingType M_PI = 3.1415926535897932384626433832795;
// Compute levels = floor(log2(n))
{
size_t temp = n;
Fourier_levels = 0;
while (temp > 1)
{
Fourier_levels++;
temp >>= 1;
}
}
// Trignometric tables
for (i = 0; i < n / 2; i++)
{
cos_table[i] = cosf(2 * M_PI * i / n);
sin_table[i] = sinf(2 * M_PI * i / n);
}
}
void FFTX::transform_radix2(FloatingType real[], FloatingType imag[], size_t n)
{
size_t size;
size_t i;
size_t x;
// Bit-reversed addressing permutation
for (i = 0; i < n; i++)
{
unsigned int k;
size_t j = 0;
x = i;
for (k = 0; k < Fourier_levels; k++, x >>= 1)
{
j = (j << 1) | (x & 1);
}
if (j > i)
{
FloatingType temp = real[i];
real[i] = real[j];
real[j] = temp;
temp = imag[i];
imag[i] = imag[j];
imag[j] = temp;
}
}
// Cooley-Tukey decimation-in-time radix-2 FFT
for (size = 2; size <= n; size *= 2)
{
size_t halfsize = size / 2;
size_t tablestep = n / size;
for (i = 0; i < n; i += size)
{
size_t j;
size_t k;
for (j = i, k = 0; j < i + halfsize; j++, k += tablestep)
{
FloatingType tpre = real[j+halfsize] * cos_table[k] + imag[j+halfsize] * sin_table[k];
FloatingType tpim = -real[j+halfsize] * sin_table[k] + imag[j+halfsize] * cos_table[k];
real[j + halfsize] = real[j] - tpre;
imag[j + halfsize] = imag[j] - tpim;
real[j] += tpre;
imag[j] += tpim;
}
}
if (size == n) // Prevent overflow in 'size *= 2'
{
break;
}
}
}
void FFTX::FFT(FloatingType real[], FloatingType imag[])
{
unsigned int I;
transform_radix2(real, imag, FourierBase_);
for (I = 0; I < FourierBase_; I++)
{
real[I] = real[I] / FourierBase_;
imag[I] = imag[I] / FourierBase_;
}
}
void FFTX::IFFT(FloatingType real[], FloatingType imag[])
{
transform_radix2(imag, real, FourierBase_);
}
FloatingType FFTX::Magnitude(FloatingType real, FloatingType imag)
{
return sqrt((real * real) + (imag * imag));
}
FloatingType FFTX::Phase(FloatingType real, FloatingType imag)
{
if (((real * real) + (imag * imag)) < 0.00001)
{
return 0;
}
else
{
return atan2(real, imag);
}
}