Изучение БПФ - Почему это не быстро?
Я не уверен, если это больше математика или больше вопросов программирования. Если это математика, пожалуйста, скажи мне.
Я знаю, что есть много готовых к использованию бесплатных проектов FFT. Но я пытаюсь понять метод БПФ. Просто для удовольствия и для изучения. Поэтому я сделал оба алгоритма - DFT и FFT, чтобы сравнить их.
Но у меня проблема с моим БПФ. Кажется, нет большой разницы в эффективности. Мой FFT только немного быстрее, чем DFT (в некоторых случаях это в два раза быстрее, но это максимальное ускорение)
В большинстве статей о БПФ есть кое-что о инверсии битов. Но я не вижу причины использовать бит-реверс. Вероятно, это так. Я не понимаю это Пожалуйста, помогите мне. Что я делаю не так?
Это мой код (вы можете скопировать его здесь и посмотреть, как он работает - онлайн-компилятор):
#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; // I devide Sample Rate to make the loop smaller,
//just to perform tests faster
float bufferSize = 512;
// Clock class is for measure time to execute whole loop:
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;
};
// Function to calculate magnitude of complex number:
float _mag_Hf(std::complex<float> sf);
// Function to calculate exp(-j*2*PI*n*k / sampleRate) - where "j" is imaginary number:
std::complex<float> _Wnk_Nc(float n, float k);
// Function to calculate exp(-j*2*PI*k / sampleRate):
std::complex<float> _Wk_Nc(float k);
int main() {
float scaleFFT = 512; // devide and conquere - if it's "1" then whole algorhitm is just simply DFT
// I wonder what is the maximum of that value. I alvays thought it should be equal to
// buffer size (number o samples) but above some value it start to work slower then DFT
std::vector<float> inputSignal; // array of input signal
inputSignal.resize(bufferSize); // how many sample we will use to calculate Fourier Transform
std::vector<std::complex<float>> _Sf; // array to store Fourier Transform value for each measured frequency bin
_Sf.resize(scaleFFT); // resize it to size which we need.
std::vector<std::complex<float>> _Hf_Db_vect; //array to store magnitude (in logarythmic dB scale)
//for each measured frequency bin
_Hf_Db_vect.resize(_SRrange); //resize it to make it able to store value for each measured freq value
std::complex<float> _Sf_I_half; // complex to calculate first half of freq range
// from 1 to Nyquist (sampleRate/2)
std::complex<float> _Sf_II_half; // complex to calculate second half of freq range
//from Nyquist to sampleRate
for(int i=0; i<(int)_Sf.size(); i++)
inputSignal[i] = cosf((float)i/_Pi); // fill the input signal with some data, no matter
Clock _time; // Start measure time
for(int freqBinK=0; freqBinK < _SRrange/2; freqBinK++) // start calculate all freq (devide by 2 for two halves)
{
for(int i=0; i<(int)_Sf.size(); i++) _Sf[i] = 0.0f; // clean all values, for next loop we need all values to be zero
for (int n=0; n<bufferSize/_Sf.size(); ++n) // Here I take all samples in buffer
{
std::complex<float> _W = _Wnk_Nc(_Sf.size()*(float)n, freqBinK);
for(int i=0; i<(int)_Sf.size(); i++) // Finally here is my devide and conquer
_Sf[i] += inputSignal[_Sf.size()*n +i] * _W; // And I see no reason to use any bit reversal, how it shoul be????
}
std::complex<float> _Wk = _Wk_Nc(freqBinK);
_Sf_I_half = 0.0f;
_Sf_II_half = 0.0f;
for(int z=0; z<(int)_Sf.size()/2; z++) // here I calculate Fourier transform for each freq
{
_Sf_I_half += _Wk_Nc(2.0f * (float)z * freqBinK) * (_Sf[2*z] + _Wk * _Sf[2*z+1]); // First half - to Nyquist
_Sf_II_half += _Wk_Nc(2.0f * (float)z *freqBinK) * (_Sf[2*z] - _Wk * _Sf[2*z+1]); // Second half - to SampleRate
// also don't see need to use reversal bit, where it shoul be??? :)
}
// Calculate magnitude in dB scale
_Hf_Db_vect[freqBinK] = _mag_Hf(_Sf_I_half); // First half
_Hf_Db_vect[freqBinK + _SRrange/2] = _mag_Hf(_Sf_II_half); // Second half
}
std::cout << _time.secondsElapsed() << std::endl; // time measuer after execution of whole loop
}
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();
return 20*log10(pow(_Re_2 + _Im_2, 0.5f)); //transform magnitude to logarhytmic dB scale
}
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;
}
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;
}
3 ответа
Одна огромная ошибка, которую вы делаете, это вычисление веса бабочки (которое включает в себя sin а также cos) на лету (в _Wnk_Nc()). sin а также cos как правило, стоимость от 10 до 100 с тактов, в то время как другие операции типа "бабочка" - просто "множитель" и "сложение", которые занимают всего несколько циклов, отсюда и необходимость их учитывать. Все быстрые реализации FFT делают это как часть шага инициализации (обычно называемого "созданием плана" или подобным). Смотрите, например, FFTW и KissFFT.
Кроме вышеупомянутой оптимизации "предварительного расчета весов бабочек", большинство реализаций FFT также используют SIMD инструкции по векторизации кода.
// также не вижу необходимости использовать бит разворота, где он должен быть?
Самая первая петля бабочки должна быть проиндексирована обратным битом. Эти индексы обычно рассчитываются внутри рекурсии, но для расчета по циклу эти индексы также являются дорогостоящими, поэтому их лучше предварительно рассчитать в плане.
Сочетание этих подходов к оптимизации приводит к ускорению примерно в 100 раз
В большинстве быстрых реализаций FFT либо используется таблица поиска предварительно вычисленных коэффициентов твида, либо простая рекурсия для вращения коэффициентов твида на лету вместо вызова тригонометрических функций библиотеки математики внутри внутреннего цикла FFT.
Для больших БПФ использование рекурсивной формулы триггера с меньшей вероятностью разрушает кэш данных на современных процессорах.