Comments 22
В итоге реализация fft_radix4_readConjSwap_simd128 обогнала родную EML (FFTW) на 18% (пусть это и 2 мкс). «Безумству храбрых поём мы песню!»
Это после https://github.com/FFTW/fftw3/pull/405 ?
Тоже с стэнда Эльбрус? :)
Можно получить доступ к Эльбрус‑8СВ по первому запросу на почту office@serpas.ru
Можно скачать QEMU и кросс‑компилятор x86‑e2k и у себя на Linux (x86) спокойно компилировать и запускать e2k‑код: https://dev.mcst.ru/download/
Ё! Loop unrolling в исходном коде. Жуть. А условный #pragma unroll(8) ?
А с GNUC vector extensions как будет производительность? Насколько тайлинг или SoA могут ускорить?
Раскрутка цикла компилятором (#pragma unroll(8)) дублирует всё содержимое тела цикла, включая инструкцию bitrevd (раскрутили 8 раз - получили 8 инструкций bitrevd).
При ручной раскрутке можно оставить только одну инструкцию bitrevd.
Чем меньше инструкций в раскрученном цикле, тем лучше можно раскручивать дальше.
При прочих равных, лучше использовать #pragma unroll , т.к. и код понятнее, и учёт некратных случаев компилятором добавляется.
Этот момент недостаточно акцентирован в статье.
Надо добавить некие размышления об этом перед 3. reverse_radix2_x2_bad.
function fft(input, inverse = false) {
const result = _fft(input, inverse);
// Normalize only once, at the top level
if (inverse) {
const N = result.length;
for (let i = 0; i < N; i++) {
result[i].re /= N;
result[i].im /= N;
}
}
return result;
}
function _fft(input, inverse) {
const N = input.length;
if (N <= 1) return input.map(x =>
typeof x === 'number' ? { re: x, im: 0 } : { re: x.re, im: x.im }
);
if ((N & (N - 1)) !== 0)
throw new Error('Размер массива должен быть степенью двойки');
const complexInput = input.map(x =>
typeof x === 'number' ? { re: x, im: 0 } : { re: x.re, im: x.im }
);
const even = complexInput.filter((_, i) => i % 2 === 0);
const odd = complexInput.filter((_, i) => i % 2 === 1);
// ✅ Pass inverse down recursively
const evenFFT = _fft(even, inverse);
const oddFFT = _fft(odd, inverse);
const output = new Array(N);
for (let k = 0; k < N / 2; k++) {
const angle = (inverse ? 1 : -1) * 2 * Math.PI * k / N;
const cos = Math.cos(angle);
const sin = Math.sin(angle);
const tRe = oddFFT[k].re * cos - oddFFT[k].im * sin;
const tIm = oddFFT[k].re * sin + oddFFT[k].im * cos;
output[k] = { re: evenFFT[k].re + tRe, im: evenFFT[k].im + tIm };
output[k + N / 2] = { re: evenFFT[k].re - tRe, im: evenFFT[k].im - tIm };
}
// ✅ No normalization here — handled by the public wrapper
return output;
}
// Создаем тестовый сигнал
const N = 8;
const signal = [];
for (let n = 0; n < N; n++) {
const value = 2 * Math.sin(2 * Math.PI * 1 * n / N);
signal.push(value);
}
console.log('\nSignal');
console.log(signal.map(v => v.toFixed(4)));
console.log('\n=== FFT (N=8) ===');
const start = performance.now();
const fftResult = fft(signal, false);
fftResult.forEach((c, i) => {
console.log(`[${i}] ${c.re.toFixed(4)} + ${c.im.toFixed(4)}i`);
});
console.log('\n=== Back FFT (N=8) ===');
const bfft = fft(fftResult, true);
bfft.forEach((c, i) => {
const diff = Math.abs(c.re - signal[i]);
console.log(`[${i}] ${c.re.toFixed(4)} (было ${signal[i].toFixed(4)}) ${diff < 0.001 ? '✅' : '❌'}`);
});
const time = performance.now() - start;
console.log(`Time: ${time.toFixed(4)} ms`);
Сделал на javascript, мне тоже печатать пост?
Скрытый текст
За сколько тактов выполняется?
Учитывая то что это делается на более высоком уровне программирования то больше чем на сишке, а учитывая то что я могу этим комментарием словить еще кучу минусов в свою карму, то да я в более плохом положении чем думал, когда решил по фану выложить данный комментарий, но я точно не ожидал что будет столько минусов к нему даже без единого плюса - хотя я не смог осилить данный пост, потому что в нём нет короткого единого кода как сделать fft и я понимаю что точнее не понимаю почему мой код который очень просто открыть в браузере и изменять не понравился местной специфичной публике которой якобы проще из данной портянки вытащить код для преобразования фурье если они хотя бы смогут это сделать - не самый быстрый , а хотя бы какой-то работающий
А почему в ассемблерном выхлопе операции над двойными и квадро-числами делаются на одинарных регистрах? Или это (очередно) баг вывода lcc -S ?
Обычно для оптимизации работы с данными делают именно чтение линейным, а запись можно раскидывать по адресам, иначе какая раскрутка и какой префетч когда у нас по большому счету 4ре массива в кэше толкаются.
Ну и вместо всех вот этих бесполезных прагм лучше указывать restrict const ✶data_in, restrict ✶data_out чтоб компилятор мог примерно понять что мы делаем с данными в рантайме.
Так выводит lcc -S (я только табы на пробелы заменил).
В Reverse мы либо читаем подряд и пишем хаотично, либо читаем хаотично и пишем подряд.
Таково свойство перестановки Reverse («bit reversal»).
Чтение подряд приводит к тормозам при записи.
Хаотичное чтение не даёт использовать APB.
Компромиссом стало чтение подряд из нескольких мест.
Так и APB работает, и запись делается большими кусками.
Хотя эти куски кладутся в память всё так же хаотично.
Данные читаются из памяти в буфер APB.
Не понятно, используется ли тут кэш L1 вообще.
Получается APB умеет по четырем адресам в память ходить. Скорей всего это только в пределах одного общего массива работает, но все равно любопытно.
Про кэш я конечно же ересь сказал, оно же не префетчер в обычном понимании, и само никуда (кроме своего буфера) ничего не подкачивает, так что ничего там не засоряется и не толкается.
Стены инструкций в конце это жесть конечно, при этом удивляет что библиотека EML местами умеет еще быстрее. Так что плотность кода это далеко не все.
Правильно ли я понял, что лучший результат для БПФ 4К комплексных float это 15078 тактов? Сколько тактов такое же БПФ на С66?
Когда-то я работал с DSP К1967ВН44. Для ориентира , у него БПФ с такими же числами для 4К- 23420 тактов, 16К - 117930, 64К-635000. Судя по всему у Вас еще не самый лучший результат получился.
Да, на данный момент, вариант fft_radix4_readConjSwap является самым лучшим.
Это развитие варианта fft_radix4 путём уменьшения числа инструкций в цикле за счёт большего числа чтений из памяти. Хорошо работает, если данные помещаются в кэш.
Конкретное число тактов плавает от запуска к запуску.
Здесь ещё не показан вариант прямой векторизации, который тоже сделан..
Там получилось похуже из-за необходимости конвертации входных данных в удобную форму для векторизации. Если данные исходно давать в удобной форме, то скорость должна быть примерно такой же (не проверял).
Надо будет добавить его в статью когда-нибудь. Логически он тоже развивается из fft_radix4, но другим путём.
Про C66 не могу сказать, потому что там не делалось отдельно FFT, там сразу делалась задача поиска картинки в картинке для 512*512 точек (т.е. было FFT для 2^18 и 2^9 точек).
И замеры тактов там не производились, сразу считались микросекунды.
Помню, что там была реализована обработка 4 комплексных чисел за 5 тактов (32/5 байтов за такт на каждом DSP). Там не был реализован вариант 2x. Был сделан аналог fft_radix2.
Там оптимизировалась вся задача целиком и ускорение этой части не было в приоритете.
Все время забываю про ИИ. Он пишет, что 4К бпф 66-й делает за 26000-29000 тактов при самом лучшем расположении данных. Но меня больше сразила цифра тактов для 16К - 328000. Или 131000 если использовать реализацию на все ядра. Но , конечно, веры в эти цифры нет :)
Есть статья про C66, в которой авторы утверждают, что можно сделать FFT для 4K точек за 27500 тактов.
Есть исходники?
Выложил у себя на Google Drive.
Компилирую так:./compile.sh */*.c
Запускаю проверку корректности так:./check_all.sh 12 1000 x
Вывод всех результатов сразу:
taskset 1 ./speed_all_stage.sh 12 10000 x
taskset 1 ./speed_all_fft.sh 12 10000 x
Данные для графиков генерирую так:taskset 1 ./data_for_diagram.py ./fft_radix2/fft_radix2_speedtest 4 21 1 1000 -1 > out.txt
Пишу алгоритм FFT на Си для процессора Эльбрус