Cześć,
Nie jestem pewien, czy to nie jest bardziej matematyczne pytanie niż programistyczne. Jeśli to nadaje się bardziej na forum matematyków to powiedzcie mi.
Wiem, że jest wiele darmowych i gotowych do użycia algorytmów FFT. Ale po prostu próbuję zrozumieć FFT. Dla zabawy, nauki, z ciekawości
.
Zrobiłem więc oba algorytmy - DFT i FFT, aby je porównać.
I mam problem. Wydaje się, że nie ma dużej różnicy w wydajności między DFT a FFT. Mój FFT jest tylko trochę szybszy niż DFT (w niektórych przypadkach jest dwa razy szybszy, ale jest to maksymalne przyspieszenie)
W większości artykułów o FFT jest coś o odwracaniu bitów. Ale nie widzę powodu, aby używać odwrócenia bitowego.
Pewnie dlatego, że po prostu nie rozumiem tych kwestii. Proszę pomóż mi. Co robię źle?
To jest mój kod: (sprowadzony do main(), gdyż chodzi mi tylko o metodologie przeprowadzenia FFT, ten kod można po prostu skopiować i wkleić do byle jakiego kompilatora online i zadziała, np: https://wandbox.org/permlink/c6NcMzcp78KN7cCT
#include <complex>
#include <iostream>
#include <math.h>
#include <cmath>
#include <vector>
#include <chrono>
#include <ctime>
float _Pi = 3.14159265;
float sampleRate = 44100;
float resolution = 4;
float _SRrange = sampleRate / resolution; // Dziele sample rate aby zrobic swoj for loop mniejszy
// Po prostu aby sprawniej testowac
float bufferSize = 512;
// Clock to klasa, ktora skopiowalem od innego kolegi z tego forum.
// Sluzy do tego aby mierzyc czas wykonania calego for loopa:
class Clock
{
public:
Clock() { start = std::chrono::high_resolution_clock::now(); }
~Clock() {}
float secondsElapsed()
{
auto stop = std::chrono::high_resolution_clock::now();
return std::chrono::duration_cast<std::chrono::microseconds>(stop - start).count();
}
void reset() { start = std::chrono::high_resolution_clock::now(); }
private:
std::chrono::time_point<std::chrono::high_resolution_clock> start;
};
// Funkcja do przeliczania liczb zespolonych na "magnitude":
float _mag_Hf(std::complex<float> sf)
{
float _Re_2;
float _Im_2;
_Re_2 = sf.real() * sf.real();
_Im_2 = sf.imag() * sf.imag();
// Od razu przeliczam magnitude na skale logarytmiczna w dB
return 20*log10(pow(_Re_2 + _Im_2, 0.5f));
}
// Funkcja do przeliczania exp(-j*2*PI*n*k / sampleRate) - gdzie j to liczba urojona:
std::complex<float> _Wnk_Nc(float n, float k)
{
std::complex<float> _Wnk_Ncomp;
_Wnk_Ncomp.real(cosf(-2.0f * _Pi * (float)n * k / sampleRate));
_Wnk_Ncomp.imag(sinf(-2.0f * _Pi * (float)n * k / sampleRate));
return _Wnk_Ncomp;
}
// Funkcja do przeliczenia exp(-j*2*PI*k / sampleRate):
std::complex<float> _Wk_Nc(float k)
{
std::complex<float> _Wk_Ncomp;
_Wk_Ncomp.real(cosf(-2.0f * _Pi * k / sampleRate));
_Wk_Ncomp.imag(sinf(-2.0f * _Pi * k / sampleRate));
return _Wk_Ncomp;
}
int main() {
float scaleFFT = 512; // "Dziel i rzadz" - jezeli wstawimy 1 to bedzie po prostu DFT
// Zastanawiam sie nad maksymalna wartoscia tej zmiennej, zawsze myslalem,
// ze to moze byc max tyle samo co ilosc sampli (buffer size)
// Ale powyzej pewnej wielkosci moj loop zaczyna chodzic wolniej
// Najlepsze osiagi mam przy 16, 32 i 64, przy buffer size 512
std::vector<float> inputSignal; // ciag wartosci wejsciowych
inputSignal.resize(bufferSize); // ustawiam jego wielkosc na rowna ilosci pobieranych sampli
std::vector<std::complex<float>> _Sf; // ciag do przechowywania Transformat Fouriera dla poszczegolnych czestotliwosci
_Sf.resize(scaleFFT); // ilosc pol jest rowna "Dziel i rzadz" czyli scaleFFT.
std::vector<std::complex<float>> _Hf_Db_vect; //Ciag do przechowywania magnitude (w skali logarytmicznej - dB)
_Hf_Db_vect.resize(_SRrange); // ilosc pol musi byc rowna ilosci mierzonych czestotliwosci
std::complex<float> _Sf_I_half; // Transformata Fouriera dla pierwdzej polowy zakresu czestotliwosci
// od 1Hz do Nyquist czyli Sample Rate /2
std::complex<float> _Sf_II_half; // Transformata Fouriera dla drugiej polowy zakresu czestotliwosci
// od Nyquist do Sample rate
for(int i=0; i<(int)_Sf.size(); i++)
inputSignal[i] = cosf((float)i*_Pi); // wypelniam sygnal wejsciowy jakimikolwiek danymi
Clock _time; // rozpoczynam pomiar czasu
/////////////////////////////////////////////////////////////////////////////////
// zaczynam przeliczac wszystkie czestotliwosci (podzielone na pol - Nyquist) ///
/////////////////////////////////////////////////////////////////////////////////
for(int freqBinK=0; freqBinK < _SRrange/2; freqBinK++)
{
// zeruje wszystkie wartosci - wymagane przez nadchodzacego loopa
for(int i=0; i<(int)_Sf.size(); i++)
{
_Sf[i] = 0.0f;
}
// Pobieram wszystkie sample wygnalu wejsciowego i je obrabiam
for (int n=0; n<bufferSize/_Sf.size(); ++n)
{
std::complex<float> _W = _Wnk_Nc(_Sf.size()*(float)n, freqBinK);
// Finalnie tu jest chyba moje "dziel i rzadz"
// I nie widze tutaj potrzeby odwracania bitow, czy cos przeoczylem?
for(int i=0; i<(int)_Sf.size(); i++)
{
_Sf[i] += inputSignal[_Sf.size()*n +i] * _W;
}
}
std::complex<float> _Wk = _Wk_Nc(freqBinK);
_Sf_I_half = 0.0f;
_Sf_II_half = 0.0f;
// przeliczenie transformaty Fouriera dla aktualnej czestotliwosci freqBinK
for(int z=0; z<(int)_Sf.size()/2; z++)
{
_Sf_I_half += _Wk_Nc(2.0f * (float)z * freqBinK) * (_Sf[2*z] + _Wk * _Sf[2*z+1]); // pierwsza polowa - do Nyquist
_Sf_II_half += _Wk_Nc(2.0f * (float)z * freqBinK) * (_Sf[2*z] - _Wk * _Sf[2*z+1]); // druga polowa - do SampleRate
// tutaj rowniez nie widze potrzeby odwracania bitow, to gdzie to do cholery ma byc? :)
}
// Obliczam magnitude w skali logarytmicznej - dB
_Hf_Db_vect[freqBinK] = _mag_Hf(_Sf_I_half); // pierwsza polowa czestotliwosci
_Hf_Db_vect[freqBinK + _SRrange/2] = _mag_Hf(_Sf_II_half); // druga polowa czestotliwosci
}
std::cout << _time.secondsElapsed() << std::endl; // pomiar czasu po przeleceniu wszystkich obliczen
}
Za każdą pomoc z góry dziękuję
Z poważaniem