Ускорение функции расчета автокорреляционной функции

@DimKaKiber · Регистрация: 20.11.2013

Студворк — интернет-сервис помощи студентам

Доброго времени суток!
Столкнулся с необходимостью ускорения расчета автокорреляционной функции.

Сейчас она реализуется следующим образом:

C++
1
2
3
4
5
6
7
8
9
    float *ACF = NULL; ACF = (float*)malloc(sizeof(float)*Win);
    for (long j = 0; j < Win; j++)
    {
        ACF[j] = 0.0;
        for (long k = 0; k < Win; k++)
        {
            ACF[j] += In[2*Win-1-k]*In[2*Win-1-k-j];
        }
      }

Реализация функции при помощи быстрого преобразования Фурье не дает идентичных результатов (картинка в приложении).

Делаю таким образом:
1. Заполняю вещественную часть 2*Win значениями вектора. Мнимая нулевая (сигнал сугубо вещественный).
2. Прямое преобразование Фурье.
3. Возвожу в квадрат мнимую и вещественную части результата (2).
4. Обратное преобразование Фурье.
5. Записываю Win значений в выходной массив и рисую его.

Может быть я не понимаю методику или вообще не в ту сторону смотрю?
Буду рад любой помощи и заранее спасибо!

@karaulov6 · 06.04.2017, 18:29

А другое место в коде ускорить не получилось?

@DimKaKiber · 07.04.2017, 05:17 **[ТС]**

Другие места в коде либо уже ускорены, либо ускоряются сейчас.
Это место больше всех временного вклада дает на текущий момент.
По времени - в районе 120-160 мс. на машине с 4 ядрами по 3.1 ГГц (IntelCore). Это примерно 10-15 вызовов операции в обвязке необходимой.
На среднем целевом устройстве время операций повышается примерно в 20 раз (просчитывали специально). Поэтому стараюсь оптимизировать.

@Manowar · 07.04.2017, 06:17

Попробуй использовать функции потокового SIMD расширения. Тем более AVX на таком компе поддерживает. С типом float сразу по 8 чисел обрабатывать будешь.

@DimKaKiber · 07.04.2017, 06:24 **[ТС]**

Есть небольшая проблема, которую я не озвучил (за что дико извиняюсь).
Оптимизируется код, который входит в состав библиотеки, включаемой в приложение, работающее на смартфонах. Пока что целевая ОС Android (4.4 и выше). Требований к минимальном характеристикам железа нет (должно работать на сферическом смарте в вакууме). Поэтому цель - использовать максимально простые функции. Без фанатизма конечно.

@obivan · 07.04.2017, 06:29

может многопоточность подрубить, разбить все циклы на небольшие блоки и пустить вычисления в отдельные потоки

@DimKaKiber · 07.04.2017, 09:50 **[ТС]**

obivan, Тогда, к сожалению, встанет проблема в корректном совмещении результатов. Но это решаемо.....
Просто пока хочется по максимуму оптимизировать для одного потока.

@Mirmik · 07.04.2017, 11:05

Можно посмотреть, какой результат дадут очевидные оптимизации.

Как минимум вы можете поменять местами циклы и выборку 2*Win-1-k делать один раз для прохода по массиву j. 2*Win-1 следует вычислять до входа в цикл. 2*Win-1-k один раз на итерацию внешнего цикла. Не думаю, что оптимизатор это учитывает.

Кроме того, вы, ИМХО (не до конца уверен), 2 раза считаете одно и тоже. Можно ограничиться половиной развёртки, после чего ее отзеркалить в отрицательную область.

Добавлено через 4 минуты
ПС. Почему 2*Win?

P.P.S. Про два раза одно и тоже фигню сказал... Но что-то мне все равно тут не нравится.

Добавлено через 3 минуты
А, понял, почему 2*Win.

Добавлено через 1 минуту
Впрочем, не слушайте меня...

Добавлено через 14 секунд
Впрочем, не слушайте меня...

@DimKaKiber · 07.04.2017, 11:34 **[ТС]**

Mirmik, Действительно даже такие минимальные правки позволили на 5-9 мс выполнение функции ускорить:

C++
1
2
3
4
5
6
7
8
9
10
11
   float *ACF = (float*)calloc(Win,sizeof(float));
    int index = 2*Win-1;
    int index1;
    for (long j = 0; j < Win; j++)
     {
      index1 = index-j;
      for (long k = 0; k < Win; k++)
       {
        ACF[j] += In[index-k]*In[index1-k];
       }
     }

С индексами проглядел. На самом деле вообще никогда так не делаю как раз из соображений быстродействия. Этот код совсем из древнего моего проекта взят был.

Я не понял насчет отзеркаливания, если честно.
Обработке же 2*Win отсчетов вектора подвергается, следовательно на выходе я должен получить Win отсчетов автокорреляционной функции. Что и делается.

И что значит

Сообщение от Mirmik

Можно ограничиться половиной развёртки, после чего ее отзеркалить в отрицательную область.

?

Это значит, что вторая половина заполняется теми же значениями, что и первая? Тогда как быть с энергией сигнала, которая в начальных точках автокорреляционной функции неизбежно будет?
Или же надо в нечетные отсчеты запихать результаты расчета?

Добавлено через 3 минуты
Мне вот не нравится то, что при выполнении самой очевидной оптимизации через БПФ результат расчета автокорреляции существенно разнится от результата прямой моей реализации.
Где то накосячить мог, а где уже и непонятно)

@nonedark2008 · 07.04.2017, 11:41

Сообщение от DimKaKiber

Возвожу в квадрат мнимую и вещественную части результата

Я в Фурье совсем не разбираюсь, но разве в квадрат возводятся не сами комплексные числа?

Добавлено через 4 минуты

Сообщение от DimKaKiber

ACF[j] += In[index-k]*In[index1-k];

Попробуй массивы пробегать в противоположном направлении, авось кэш-промахов станет меньше.

@Mirmik · 07.04.2017, 12:01

Про отрицательную часть я фигню сморозил. Не вкурил поначалу в то, как оно работает.

Добавлено через 7 минут
Можете попробовать поменять местами порядок обхода. То есть раньше внутренний цикл обходился по k, то теперь по j. Из-за этого доступ к In[index-k] можно осуществить во внешнем цикле.

Я не уверен, что это как-то ускорит выполнение из-за индексирования ACF[j], но можно потестировать.

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int index = 2*Win-1;
int index_k;
float InK;
 
for (long k = 0; k < Win; k++)
{   
     InK = In[index-k]; 
     index_k = index - k;
     
     for (long j = 0; j < Win; j++)
     {
           ACF[j] += InK * In[index_k - j];
     }
}

P.S. идея с проходом в прямом направлении любопытна.

@DimKaKiber · 07.04.2017, 12:01 **[ТС]**

nonedark2008, По разному уже пробовал, если честно.
Но согласно https://ru.wikipedia.org/wiki/... 0%B8%D1%8F, "Цифровая обработка речевых сигналов" (Рабинер, Шафер) и ряду других талмудов по цифровой обработке сигналов нужно возвести в квадрат по отдельности вещественную и мнимую части результата прямого БПФ, сложить их и обратное преобразование произвести.

Mirmik, Я уже после того как ответ написал увидел, что фигня

@Mirmik · 07.04.2017, 12:08

П.С. А вообще, по хорошему, возможно, стоит сеть и под каждую из ваших архитектур написать это дело на ассемблере.

@DimKaKiber · 07.04.2017, 12:10 **[ТС]**

Mirmik, спасибо большое! Еще 1 мс в среднем удалось выиграть)

Ассемблер, к сожалению, еще изучать придется на ходу......

MrGluck · 07.04.2017, 13:06

C++
1
float *ACF = (float*)calloc(Win, sizeof(float));

и можно убрать строку

Сообщение от DimKaKiber

C++
1
ACF[j] = 0.0;

Добавлено через 1 минуту
А вообще, если Win известен на стадии компиляции/мы знаем что он не будет больше определённого значения и нам не жалко памяти - используйте статический массив. Тогда избежите выделения динамической памяти в рантайме.

@DimKaKiber · 07.04.2017, 13:10 **[ТС]**

MrGluck, Я в 9 посте исправил это)
Памяти жалко) Но на определенном этапе можно ее разок выделить, чтобы, действительно, по 5-12 раз не выделять заново)

MrGluck · 07.04.2017, 13:14

Если есть возможность - посмотрите в сторону кеширования результатов. И мб даже имеет смысл использовать ленивые вычисления (но это надо знать всю вашу архитектуру).

Добавлено через 3 минуты
DimKaKiber, операции выделения/освобождения могут быть достаточно дороги в рантайме + при частом выделении/освобождении вы фрагментируете память и со временем это может сильно замедлить работу.
Если размер (или макс. размер) не известен в момент компиляции, но не меняется со временем - выделите один раз память и не освобождайте до выхода из программы.
Можете, кстати, попробовать использовать vector (тоже сделать так, чтобы он не пересоздавался каждый раз). Он поможет ещё и избавить он лишних выделений памяти если размер нового буфера будет меньше (можете это сделать и вручную, но зачем городить велосипед с возможностью ошибиться). Оверхед от вектора минимален.

@Mirmik · 07.04.2017, 13:31

Не надо использовать контейнеры stl при низкоуровневых оптимизациях!

Выделение памяти тоже крайне сомнительно по части ускорения алгоритма. Выделение памяти - это микросекунды. А мы с миллисекундами боремся. Но, лучше, конечно, один раз выделять.

P.S. Кеширование результатов - это хорошо. Возможно вам не нужно каждый раз пересчитывать весь массив. Например, если у вас реалтайм и скользящее окно, вы можете на лету добавлять вклад от новых точек и вычитать вклад от отбрасываемых. (правда, могут ошибки округлений накапливаться).
Это можно полностью решить перейдя от вычислений с плавающей точкой к вычислениям с фиксированной.

@DimKaKiber · 07.04.2017, 13:44 **[ТС]**

MrGluck, Спасибо большое за пояснение. Я поразмышлял и решил все-таки выделить память один раз и не суетиться. При чем действительно после 25 минут работы приложение начинает "подтупливать".
Опять таки приходится постоянно заполнять нулями вектор ACF для расчетов.
Но здесь

C++
1
memset(ACF,0,sizeof(float)*Win);

в помощь)

Сейчас по всему коду пройдусь и в возможных местах память статично выделю под массивы известной размерности

Добавлено через 4 минуты
Mirmik, К сожалению именно в этой функции пересчет постоянный происходит. Другие функции, где это возможно, именно так изначально и проектировались.
А vector штука удобная, но все время приходится память за ним подчищать. После освобождения, или при добавлении новых элементов, например.

MrGluck · 07.04.2017, 13:45

Сообщение от DimKaKiber

но все время приходится память за ним подчищать.

он сам все подчищает. В этом как раз его суть - RAII.

Новые блоги и статьи Все статьи Все блоги /
сукцессия микоризы: основная теория в виде двух уравнений. anaschu 11.01.2026 https:/ / rutube. ru/ video/ 7a537f578d808e67a3c6fd818a44a5c4/	WordPad для Windows 11 Jel 10.01.2026 WordPad для Windows 11 — это приложение, которое восстанавливает классический текстовый редактор WordPad в операционной системе Windows 11. После того как Microsoft исключила WordPad из. . .	Classic Notepad for Windows 11 Jel 10.01.2026 Old Classic Notepad for Windows 11 Приложение для Windows 11, позволяющее пользователям вернуть классическую версию текстового редактора «Блокнот» из Windows 10. Программа предоставляет более. . .	Почему дизайн решает? Neotwalker 09.01.2026 В современном мире, где конкуренция за внимание потребителя достигла пика, дизайн становится мощным инструментом для успеха бренда. Это не просто красивый внешний вид продукта или сайта — это. . .
Модель микоризы: классовый агентный подход 3 anaschu 06.01.2026 aa0a7f55b50dd51c5ec569d2d10c54f6/ O1rJuneU_ls https:/ / vkvideo. ru/ video-115721503_456239114	Owen Logic: О недопустимости использования связки «аналоговый ПИД» + RegKZR ФедосеевПавел 06.01.2026 Owen Logic: О недопустимости использования связки «аналоговый ПИД» + RegKZR ВВЕДЕНИЕ Введу сокращения: аналоговый ПИД — ПИД регулятор с управляющим выходом в виде числа в диапазоне от 0% до. . .	Модель микоризы: классовый агентный подход 2 anaschu 06.01.2026 репозиторий https:/ / github. com/ shumilovas/ fungi ветка по-частям. коммит Create переделка под биомассу. txt вход sc, но sm считается внутри мицелия. кстати, обьем тоже должен там считаться. . . .	Расчёт токов в цепи постоянного тока igorrr37 05.01.2026 / * Дана цепь постоянного тока с сопротивлениями и напряжениями. Надо найти токи в ветвях. Программа составляет систему уравнений по 1 и 2 законам Кирхгофа и решает её. Последовательность действий:. . .

@karaulov6 23 / 23 / 6 Регистрация: 23.03.2013 Сообщений: 245
	06.04.2017, 18:29
	А другое место в коде ускорить не получилось? 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 05:17 [ТС]
	Другие места в коде либо уже ускорены, либо ускоряются сейчас. Это место больше всех временного вклада дает на текущий момент. По времени - в районе 120-160 мс. на машине с 4 ядрами по 3.1 ГГц (IntelCore). Это примерно 10-15 вызовов операции в обвязке необходимой. На среднем целевом устройстве время операций повышается примерно в 20 раз (просчитывали специально). Поэтому стараюсь оптимизировать. 0

@Manowar 1719 / 568 / 187 Регистрация: 12.03.2016 Сообщений: 2,169
	07.04.2017, 06:17
	Попробуй использовать функции потокового SIMD расширения. Тем более AVX на таком компе поддерживает. С типом float сразу по 8 чисел обрабатывать будешь. 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 06:24 [ТС]
	Есть небольшая проблема, которую я не озвучил (за что дико извиняюсь). Оптимизируется код, который входит в состав библиотеки, включаемой в приложение, работающее на смартфонах. Пока что целевая ОС Android (4.4 и выше). Требований к минимальном характеристикам железа нет (должно работать на сферическом смарте в вакууме). Поэтому цель - использовать максимально простые функции. Без фанатизма конечно. 0

@obivan Падаван С++ 447 / 261 / 89 Регистрация: 11.11.2014 Сообщений: 916
	07.04.2017, 06:29
	может многопоточность подрубить, разбить все циклы на небольшие блоки и пустить вычисления в отдельные потоки 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 09:50 [ТС]
	obivan, Тогда, к сожалению, встанет проблема в корректном совмещении результатов. Но это решаемо..... Просто пока хочется по максимуму оптимизировать для одного потока. 0

@Mirmik techpriest 634 / 213 / 57 Регистрация: 27.02.2014 Сообщений: 1,180
	07.04.2017, 11:05
	Можно посмотреть, какой результат дадут очевидные оптимизации. Как минимум вы можете поменять местами циклы и выборку 2Win-1-k делать один раз для прохода по массиву j. 2Win-1 следует вычислять до входа в цикл. 2Win-1-k один раз на итерацию внешнего цикла. Не думаю, что оптимизатор это учитывает. Кроме того, вы, ИМХО (не до конца уверен), 2 раза считаете одно и тоже. Можно ограничиться половиной развёртки, после чего ее отзеркалить в отрицательную область. Добавлено через 4 минуты* ПС. Почему 2Win? P.P.S. Про два раза одно и тоже фигню сказал... Но что-то мне все равно тут не нравится. Добавлено через 3 минуты* А, понял, почему 2Win. Добавлено через 1 минуту* Впрочем, не слушайте меня... Добавлено через 14 секунд Впрочем, не слушайте меня... 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 12:01 [ТС]
	nonedark2008, По разному уже пробовал, если честно. Но согласно https://ru.wikipedia.org/wiki/... 0%B8%D1%8F, "Цифровая обработка речевых сигналов" (Рабинер, Шафер) и ряду других талмудов по цифровой обработке сигналов нужно возвести в квадрат по отдельности вещественную и мнимую части результата прямого БПФ, сложить их и обратное преобразование произвести. Mirmik, Я уже после того как ответ написал увидел, что фигня 0

@Mirmik techpriest 634 / 213 / 57 Регистрация: 27.02.2014 Сообщений: 1,180
	07.04.2017, 12:08
	П.С. А вообще, по хорошему, возможно, стоит сеть и под каждую из ваших архитектур написать это дело на ассемблере. 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 12:10 [ТС]
	Mirmik, спасибо большое! Еще 1 мс в среднем удалось выиграть) Ассемблер, к сожалению, еще изучать придется на ходу...... 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 13:10 [ТС]
	MrGluck, Я в 9 посте исправил это) Памяти жалко) Но на определенном этапе можно ее разок выделить, чтобы, действительно, по 5-12 раз не выделять заново) 0

MrGluck Форумчанин 8216 / 5047 / 1437 Регистрация: 29.11.2010 Сообщений: 13,453
	07.04.2017, 13:14
	Если есть возможность - посмотрите в сторону кеширования результатов. И мб даже имеет смысл использовать ленивые вычисления (но это надо знать всю вашу архитектуру). Добавлено через 3 минуты DimKaKiber, операции выделения/освобождения могут быть достаточно дороги в рантайме + при частом выделении/освобождении вы фрагментируете память и со временем это может сильно замедлить работу. Если размер (или макс. размер) не известен в момент компиляции, но не меняется со временем - выделите один раз память и не освобождайте до выхода из программы. Можете, кстати, попробовать использовать vector (тоже сделать так, чтобы он не пересоздавался каждый раз). Он поможет ещё и избавить он лишних выделений памяти если размер нового буфера будет меньше (можете это сделать и вручную, но зачем городить велосипед с возможностью ошибиться). Оверхед от вектора минимален. 1

@Mirmik techpriest 634 / 213 / 57 Регистрация: 27.02.2014 Сообщений: 1,180
	07.04.2017, 13:31
	Не надо использовать контейнеры stl при низкоуровневых оптимизациях! Выделение памяти тоже крайне сомнительно по части ускорения алгоритма. Выделение памяти - это микросекунды. А мы с миллисекундами боремся. Но, лучше, конечно, один раз выделять. P.S. Кеширование результатов - это хорошо. Возможно вам не нужно каждый раз пересчитывать весь массив. Например, если у вас реалтайм и скользящее окно, вы можете на лету добавлять вклад от новых точек и вычитать вклад от отбрасываемых. (правда, могут ошибки округлений накапливаться). Это можно полностью решить перейдя от вычислений с плавающей точкой к вычислениям с фиксированной. 0