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

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

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

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

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

C++

    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++

   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++

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.

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

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

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

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

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

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

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

@nonedark2008 1394 / 1023 / 325 Регистрация: 28.07.2012 Сообщений: 2,813
	07.04.2017, 11:41	10
	Сообщение от DimKaKiber Возвожу в квадрат мнимую и вещественную части результата Я в Фурье совсем не разбираюсь, но разве в квадрат возводятся не сами комплексные числа? Добавлено через 4 минуты Сообщение от DimKaKiber ACF[j] += In[index-k]*In[index1-k]; Попробуй массивы пробегать в противоположном направлении, авось кэш-промахов станет меньше. 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 12:01 [ТС]	12
	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	13
	П.С. А вообще, по хорошему, возможно, стоит сеть и под каждую из ваших архитектур написать это дело на ассемблере. 0

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

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

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

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

MrGluck Форумчанин 8215 / 5045 / 1437 Регистрация: 29.11.2010 Сообщений: 13,453
	07.04.2017, 13:45	20
	Сообщение от DimKaKiber но все время приходится память за ним подчищать. он сам все подчищает. В этом как раз его суть - RAII. 0