Это продолжение быстрого преобразования Фурье C ++. В итоге я немного изменил алгоритм, чтобы уменьшить использование памяти.
#include <vector>
#include <iostream>
#include <complex>
#include <cmath>
#include <algorithm>
#define N 1048576 // Problem size, must be a power of 2
#define PI 3.14159265358979323846
typedef std::vector<std::complex<double>> complexVec;
/*
Calculates the N'th roots of unity.
*/
complexVec rootsOfUnityCalculator() {
complexVec table;
table.resize(N);
for (size_t k = 0; k < N; k++) {
table[k] = std::complex<double>(std::cos(-2.0 * PI * k / N), std::sin(-2.0 * PI * k / N));
}
return table;
}
/*
Permutes an input vector of size length - a power of 2 - according
to the bit reversal permutation.
TESTED: SUCCESS
*/
void bitReversal(complexVec& input) {
std::vector<size_t> permutation;
permutation.resize(N);
permutation[0] = 0;
size_t k = 2;
// Calculating the permutation according to recursive relation
while (k <= N) {
for (size_t j = 0; j < k / 2; j++) {
permutation[j] *= 2; // permutation_k(j) = 2 * permutation_(k/2)(j)
permutation[j + k / 2] = permutation[j] + 1; // permutation_k(j + k/2) = 2 * permutation_(k/2)(j) + 1
}
k *= 2;
}
// Permute the input
for (size_t j = 0; j < N; j++) {
if (j < permutation[j]) {
std::swap(input[j], input[ permutation[j] ]);
}
}
}
void unorderedFFT(complexVec& input, const complexVec& table) {
size_t k = 2;
while (k <= N) {
for (size_t r = 0; r < N / k; r++) {
for (size_t j = 0; j < k / 2; j++) {
std::complex<double> plusMinus = input[r * k + j + k / 2] * table[N / k * j]; // omega_k^j = (omega_N^(N/k))^j = omega_N^(Nj/k)
input[r * k + j + k / 2] = input[r * k + j] - plusMinus;
input[r * k + j] += plusMinus;
}
}
k *= 2;
}
}
/*
Calculates the FFT
TESTED: SUCCESS
*/
void FFT(complexVec& input, const complexVec& table) {
bitReversal(input);
unorderedFFT(input, table);
}
// Prints array of N complex numbers
void printArray(complexVec& array) {
for (size_t k = 0; k < N - 1; k++) {
std::cout << array[k] << ", ";
}
std::cout << array[N - 1] << std::endl;
}
int main() {
complexVec input;
input.resize(N);
// Initialize the input to all 0, 1, ..., N
for (size_t k = 0; k < N; k++) {
input[k] = std::complex<double>(k, 0.0);
}
const complexVec table = rootsOfUnityCalculator();
FFT(input, table);
}
1 ответ
Результаты профилирования
Используя профилировщик в Visual Studio, на моем компьютере (Intel Haswell, 4770K) я получаю примерно такую разбивку:
unorderedFFT107 мсrootsOfUnityCalculator21 мсbitReversal17 мс
Так что давайте улучшим все это.
Корни единства
Корни N из единицы на комплексной плоскости — это набор равноотстоящих точек на единичной окружности. Их координаты могут быть вычислены с использованием синуса и косинуса, но альтернатива начинается с 1 + 0i и последовательно применяя к нему небольшое вращение вектора (альтернативный вид: вычисляется первый нетривиальный корень, а затем итеративно возводится в степень 0, 1, 2 .. N-1). Таким образом, необходима только одна пара sin / cos (кстати, std::polar может использоваться для замены такой пары sin / cos более простым кодом). Обратной стороной является то, что он может быть менее точным. Обычно это нормально, но это немного искажает результаты.
Например, вот так:
complexVec rootsOfUnityCalculator() {
complexVec table;
table.resize(N);
auto step = std::polar(1.0, 1.0 / N);
std::complex<double> root = 1.0;
for (std::size_t k = 0; k < N; k++) {
table[k] = root;
root *= step;
}
return table;
}
При этом время, которое требуется, сократилось до 7,5 мс, то есть коэффициент 2,3, что, безусловно, является значительным. Это можно сделать даже быстрее, чем это сделать, развернув цикл и удалив зависимость между копиями тела цикла, например, если развернуть на 2, тогда будет отдельный root1 и root2 каждый из которых умножается на stepSquared так что их вычисления независимы (улучшение параллелизма на уровне команд). Это сократило время до 6 мс.
Перестановка с обращением битов
Это хорошо изученная перестановка, потому что она представляет практический интерес, а именно это точное приложение. Известны различные более эффективные методы, включая некоторые варианты «увеличения обратного целого числа». Если целое число i и его обратное rev доступны оба, вычисление следующих значений для них обоих возможно эффективно.
Конечно i можно просто увеличить. rev требуется дополнительная работа, но есть надежда: XOR целого числа и его преемника — это «хорошая маска», состоящая из непрерывного диапазона единиц и непрерывного диапазона нулей без «беспорядка». Это является результатом того, как при приращении перенос проходит через конечные единицы, переворачивая каждую из них, затем, наконец, первый найденный ноль также переворачивается, и перенос «поглощается» этим нулем — так что биты, перевернутые, начинают наименьший значащий бит и образуют непрерывный диапазон. Это актуально, потому что «красивую маску» можно отменить, сдвинув ее, а затем применив XOR. rev это приводит к «инкременту с инвертированием битов». Вариант этого трюка заключается в использовании команды отсчета нулевого значения, чтобы определить длину непрерывного диапазона перевернутых битов, а затем основывать маску непосредственно на этом. Пойдем с этим:
void bitReversal(complexVec& input) {
std::size_t maskLen = std::countr_zero(unsigned(N));
// Permute the input
for (std::size_t j = 0, rev = 0; j < N; j++) {
if (j < rev)
std::swap(input[j], input[rev]);
std::size_t maskLen = std::countr_zero(j + 1) + 1;
rev ^= N - (N >> maskLen);
}
}
Здесь используется очень новый std::countr_zero от <bit> заголовок, функция C ++ 20. Есть и другие способы написать это, например _tzcnt_u64. unsigned(N) немного прискорбно, но требуется для countr_zero, который отказывается работать с подписанными типами. Желательно избежать этого, сделав N константа типа size_t вместо #define-ин его, что я все равно рекомендую.
Используя этот трюк, bitReversal снизился с 17 мс до 14 мс (кроме того, память, используемая permutation вектор сохраняется). Это немного лучше, но ничего лучше, чем улучшение вычисления корней единства.
Это не лучший способ сделать это. В if соответствует плохо спрогнозированной ветви, а шаблон доступа к памяти является полуслучайным. Существуют различные полезные статьи о технике, которую я использовал здесь, и о дальнейших улучшениях, например практически эффективные методы выполнения перестановки с обратным битом в C ++ 11 на архитектуре x86-64.
Фактическое БПФ
Настоящее мясо алгоритма. Поскольку в предыдущем вопросе вы просили не беспокоиться о предложениях по различным алгоритмам, я не буду. Однако я все же могу предложить улучшение производительности: используйте SSE3. SSE2 на самом деле достаточно, но SSE3 добавляет ADDSUBPD инструкция, удобная для сложного умножения:
__m128d multiplyComplex(__m128d A, __m128d B)
{
__m128d ARealReal = _mm_shuffle_pd(A, A, 0);
__m128d AImagImag = _mm_shuffle_pd(A, A, 3);
__m128d BRealImag = B;
__m128d BImagReal = _mm_shuffle_pd(B, B, 1);
return _mm_addsub_pd(_mm_mul_pd(ARealReal, BRealImag), _mm_mul_pd(AImagImag, BImagReal));
}
Что можно было бы использовать так:
void unorderedFFT2(complexVec& input, const complexVec& table) {
std::size_t k = 2;
while (k <= N) {
for (std::size_t r = 0; r < N / k; r++) {
for (std::size_t j = 0; j < k / 2; j++) {
__m128d input0 = _mm_loadu_pd((double*)&input[r * k + j + k / 2]);
__m128d input1 = _mm_loadu_pd((double*)&input[r * k + j]);
__m128d twiddle = _mm_loadu_pd((double*)&table[N / k * j]);
__m128d product = multiplyComplex(input0, twiddle);
_mm_storeu_pd((double*)&input[r * k + j + k / 2], _mm_sub_pd(input1, product));
_mm_storeu_pd((double*)&input[r * k + j], _mm_add_pd(input1, product));
}
}
k *= 2;
}
}
Приведение этих указателей выглядит устрашающе, но, согласно параграф о доступе, ориентированном на массив из std::complex. reinterpret_cast можно использовать, если вам не нравятся приведения в стиле C.
По общему признанию, использование SSE не очень удобно для новичков. Но это действительно помогает, время unorderedFFT снизился до 70 мс, что является самым большим улучшением в абсолютном выражении.
Использование таких встроенных функций требует включения соответствующего заголовка, например #include <tmmintrin.h> (или новее) в этом случае. В зависимости от компилятора может также потребоваться передача определенных параметров командной строки (GCC и Clang нуждаются в этом, а MSVC — нет).