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

@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 значений в выходной массив и рисую его.

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

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

Я сам проблему не изучал, но коллеги заморочились.
Например, при вызове

C++

1	Vector.clear();

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

Добавлено через 31 секунду
просто редко работаю с ним

MrGluck · 07.04.2017, 13:52

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

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

Это так. Но вам же это может сыграть на руку, когда захотите заново заполнить вектор. Выделения памяти снова не будет.

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

Здесь полностью согласен.

@Mirmik · 07.04.2017, 13:58

А кстати, господа, что насчет операций с фиксированной точкой? Могут ли они оказаться быстрее плавающей арифметики?

@DimKaKiber · 07.04.2017, 14:04 **[ТС]**

Mirmik, Я пока понять не могу как обойти преобразование числа с плавающей точкой в число с фиксированной. нагрузка то все-равно дополнительная появится.......Хотя она может на самом деле ускорить процесс.
Останавливает еще то, что мне результат этой функции один фиг придется обратно в float переводить, чтобы дальше расчеты продолжить.

@Mirmik · 07.04.2017, 14:36

Зависит от размеров вашего массива.

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

Но будет ли выигрыш, конечно требует проверки.

@nonedark2008 · 07.04.2017, 14:50

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

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

А выложи код с реализацией через ДПФ, авось кто-нибудь поймет и поможет. Можно попробовать сравнить результаты с тем, что выдает какой-нибудь мат-пакет по типу Matlab и прочих.

MrGluck · 07.04.2017, 15:13

в OpenCV были готовые функции для БПФ, сам их пользовал, довольно шустрые. Можно взять для сравнения результата.

@DimKaKiber · 08.04.2017, 15:53 **[ТС]**

nonedark2008, Делаю так:

C++

void fastACF(short *In,int Win,float *&Out)
{
 FastFourier Kernel(2*Win,-1);
 Kernel.ShortToRedat(In,2*Win); //заполняю вещественную часть значениями
 FFT_for_conv(Kernel); //прямое БПФ
 for (long i = 0; i < 2*Win; i++)
  {
   Kernel.Redat[i] *= Kernel.Redat[i];
   Kernel.Redat[i] += Kernel.Imdat[i]*Kernel.Imdat[i];
  } //получил мощность спектра
// memset(Kernel.Redat+Win,0,sizeof(float)*Win);
 memset(Kernel.Imdat,0,sizeof(float)*(2*Win));
 Kernel.SetFlag(1);
 FFT_for_conv(Kernel); //обратное БПФ
 if (Out != NULL) { free(Out); Out = 0; }
 Out = (float*)malloc(sizeof(float)*Win);
 memcpy(Out,Kernel.Redat,sizeof(float)*Win); //заполнил выходной массив
}

Добавлено через 1 час 3 минуты
Само БПФ реализовано следующим образом. Раньше никогда меня не подводили результаты этой функции. Её брал откуда то и "допиливал" под себя.

C++

bool  FFT_for_conv(FastFourier &Params) //FFT преобразование
{
#define  NUMBER_IS_2_POW_K(x)   ((!((x)&((x)-1)))&&((x)>1))  // x is pow(2, k), k=1,2, ...
#define  FT_DIRECT        -1    // Direct transform.
#define  FT_INVERSE        1    // Inverse transform.
    // #define NULL 0
 
    // parameters error check:
    if((Params.Redat == NULL) || (Params.Imdat == NULL))                  return false;
    if((Params.N > 16384) || (Params.N < 1))                            return false;
    if(!NUMBER_IS_2_POW_K(Params.N))                             return false;
    if((Params.LogN < 2) || (Params.LogN > 14))                         return false;
    if((Params.Ft_Flag != FT_DIRECT) && (Params.Ft_Flag != FT_INVERSE)) return false;
    register long  i, j, n, k, io, ie, in, nn;
    float         ru, iu, rtp, itp, rtq, itq, rw, iw, sr;
    static const float Rcoef[14] =
    {  -1.0000000000000000F,  0.0000000000000000F,  0.7071067811865475F,
       0.9238795325112867F,  0.9807852804032304F,  0.9951847266721969F,
       0.9987954562051724F,  0.9996988186962042F,  0.9999247018391445F,
       0.9999811752826011F,  0.9999952938095761F,  0.9999988234517018F,
       0.9999997058628822F,  0.9999999264657178F
    };
    static const float Icoef[14] =
    {   0.0000000000000000F, -1.0000000000000000F, -0.7071067811865474F,
        -0.3826834323650897F, -0.1950903220161282F, -0.0980171403295606F,
        -0.0490676743274180F, -0.0245412285229122F, -0.0122715382857199F,
        -0.0061358846491544F, -0.0030679567629659F, -0.0015339801862847F,
        -0.0007669903187427F, -0.0003834951875714F
    };
    nn = Params.N >> 1;
    ie = Params.N;
    for(n=1; n<=Params.LogN; n++)
    {
        rw = Rcoef[Params.LogN - n];
        iw = Icoef[Params.LogN - n];
        if(Params.Ft_Flag == FT_INVERSE) iw = -iw;
        in = ie >> 1;
        ru = 1.0F;
        iu = 0.0F;
        for(j=0; j<in; j++)
        {
            for(i=j; i< Params.N; i+=ie)
            {
                io       = i + in;
                rtp      = Params.Redat[i]  + Params.Redat[io];
                itp      = Params.Imdat[i]  + Params.Imdat[io];
                rtq      = Params.Redat[i]  - Params.Redat[io];
                itq      = Params.Imdat[i]  - Params.Imdat[io];
                Params.Redat[io] = rtq * ru - itq * iu;
                Params.Imdat[io] = itq * ru + rtq * iu;
                Params.Redat[i]  = rtp;
                Params.Imdat[i]  = itp;
            }
            sr = ru;
            ru = ru * rw - iu * iw;
            iu = iu * rw + sr * iw;
        }
        ie >>= 1;
    }
    for(j=i=1; i< Params.N; i++)
    {
        if(i < j)
        {
            io       = i - 1;
            in       = j - 1;
            rtp      = Params.Redat[in];
            itp      = Params.Imdat[in];
            Params.Redat[in] = Params.Redat[io];
            Params.Imdat[in] = Params.Imdat[io];
            Params.Redat[io] = rtp;
            Params.Imdat[io] = itp;
        }
        k = nn;
        while(k < j)
        {
            j   = j - k;
            k >>= 1;
        }
        j = j + k;
    }
    if(Params.Ft_Flag == FT_DIRECT) return true;
    rw = 1.0F / Params.N;
    for(i=0; i< Params.N; i++)
    {
        Params.Redat[i] *= rw;
        Params.Imdat[i] *= rw;
    }
    return true;
}

@DimKaKiber · 10.04.2017, 18:12 **[ТС]**

Хрень какая то твориться........
Уже по методичке сделал:

$https://www.cyberforum.ru/cgi-bin/latex.cgi?{F}^{-1}\left[F(X)*{F}^{*}(X) \right]$

Т.е. обратное БПФ от призведения спектра на комплексно-сопряженный спектр - все равно хрень получается((

Уже тупо все в лоб закодил....

C++

void fastACF(short *In,int Win,float *&Out)
{
 FastFourier Kernel(2*Win,-1);
 Kernel.ShortToRedat(In,2*Win);
 FFT_for_conv(Kernel);
 FastFourier X(2*Win,-1);
 X.Zeros();
 for (long i = 0; i < 2*Win; i++)
  {
   X.Redat[i] = Kernel.Redat[i];
   X.Imdat[i] *=-1;
  }
 float a_re,b_im;
 for (long i = 0; i < Kernel.N; i++)
 {
     a_re = Kernel.Redat[i]*X.Redat[i]-Kernel.Imdat[i]*X.Imdat[i];
     b_im = Kernel.Imdat[i]*X.Redat[i]+Kernel.Redat[i]*X.Imdat[i];
     Kernel.Redat[i] = a_re;
     Kernel.Imdat[i] = b_im;
 }
 Kernel.SetFlag(1);
 FFT_for_conv(Kernel);
 Out = (float*)malloc(sizeof(float)*Win);
 memcpy(Out,Kernel.Redat,sizeof(float)*Win);
}

Может я туплю и функция у меня вообще не АКФ считает? Может заблуждался все это время.... Уже ни в чем не сомневаюсь, если честно((((

@DimKaKiber · 10.04.2017, 18:21 **[ТС]**

опечатался - не

C++

1	X.Imdat[i] *=-1;

,

а

C++

1	X.Imdat[i] =-1Kernel.Imdat[i];

@nonedark2008 · 10.04.2017, 23:53

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

Т.е. обратное БПФ от призведения спектра на комплексно-сопряженный спектр - все равно хрень получается

Ну, в формуле проблем не было. Умножение числа на комплексно-сопряженное и дает в результате сумму квадратов мнимой и действительной частей. Я уже давно советовал проверить расчеты в какой-нибудь мат-библиотеке: Matlab, MathCAD, Maple...
Может вообще оказаться, что у тебя не ответ, а одни сплошные ошибки округления. Какой из результатов вообще больше похож на правду?

@DimKaKiber · 11.04.2017, 12:47 **[ТС]**

nonedark2008, Проверял. результаты немного другие получаются. Сейчас разбираюсь сижу где все-таки не так сделал.

Вот результат именно расчета автокорреляции. Только работает он в 2 раза быстрее исходного алгоритма.

C++

     for (int i = 0; i < Win; i++)
       {
        for (int j = 0; j < Win-i; j++)
          {
           ACF[i]+=In[j]*In[j+i];
          }
        }

@DimKaKiber · 11.04.2017, 16:26 **[ТС]**

В последнем посте у меня корреляция куска с самим собой на участке Win считалась. Для куска размером 2*Win поправил так:

C++

    for (int i = 0; i < Win; i++)
    {
        for (int j = 0; j < Win; j++)
        {
            ACF[i] += In[j]*In[j+i];
        }
    }

При этом удалось все это дело ускорить при помощи быстрого преобразования Фурье.
Время работы функции с обвязкой упало с 25-140 мс на моей машине до 1-25 мс. в зависимости от количества необходимых операций (постоянно плавает). Результат вычислений испытания на модельных и "полевых" наборах данных выдержал. Результатом доволен - теперь полезу в обвязку

Ускорение через быстрое преобразование Фурье реализовал как свертку комплексно-сопряженного спектра части сигнала размером Win со спектром куска размером 2*Win.

C++

void fastACF(short *In,int Win,float *&Out)
{
 if (Out != NULL) { free(Out); Out = 0; }
 Out = (float*)malloc(sizeof(float)*Win);
 FastFourier Kernel(2*Win,-1);
 Kernel.ShortToRedat(In,2*Win);
 FFT_for_conv(Kernel);
 FastFourier S(2*Win,-1);
 S.ShortToRedat(In,Win);
 FFT_for_conv(S);
 float a_re,b_im;
 for (long i = 0; i < Kernel.N; i++)
 {
     S.Imdat[i] *= -1;
     a_re = Kernel.Redat[i]*S.Redat[i]-Kernel.Imdat[i]*S.Imdat[i];
     b_im = Kernel.Imdat[i]*S.Redat[i]+Kernel.Redat[i]*S.Imdat[i];
     S.Redat[i] = a_re;
     S.Imdat[i] = b_im;
 }
 S.SetFlag(1);
 FFT_for_conv(S);
 memcpy(Out,S.Redat,sizeof(float)*Win);
}

Всем большое спасибо за помощь и дельные советы!

MrGluck Форумчанин 8215 / 5045 / 1437 Регистрация: 29.11.2010 Сообщений: 13,453
	07.04.2017, 13:52	22
	Сообщение от DimKaKiber количество элементов становится нулевым, но размер выделенной памяти остается прежним. Возможно это и не так. Я, повторяюсь, в это не углублялся Это так. Но вам же это может сыграть на руку, когда захотите заново заполнить вектор. Выделения памяти снова не будет. 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 13:56 [ТС]	23
	Здесь полностью согласен. 0

@Mirmik techpriest 634 / 213 / 57 Регистрация: 27.02.2014 Сообщений: 1,180
	07.04.2017, 13:58	24
	А кстати, господа, что насчет операций с фиксированной точкой? Могут ли они оказаться быстрее плавающей арифметики? 0

@DimKaKiber 12 / 12 / 0 Регистрация: 20.11.2013 Сообщений: 167
	07.04.2017, 14:04 [ТС]	25
	Mirmik, Я пока понять не могу как обойти преобразование числа с плавающей точкой в число с фиксированной. нагрузка то все-равно дополнительная появится.......Хотя она может на самом деле ускорить процесс. Останавливает еще то, что мне результат этой функции один фиг придется обратно в float переводить, чтобы дальше расчеты продолжить. 0

@Mirmik techpriest 634 / 213 / 57 Регистрация: 27.02.2014 Сообщений: 1,180
	07.04.2017, 14:36	26
	Зависит от размеров вашего массива. Ваш алгоритм имеет квадратичную сложность, а преобразование - линейную. Чем больше массив, тем несопоставимей временные интервалы. Если у вас хотя-бы пятьдесят точек в этом массиве, вы затрат времени на преобразования даже незаметите по сравнению со временем вычисления развёртк Но будет ли выигрыш, конечно требует проверки. 0

@nonedark2008 1394 / 1023 / 325 Регистрация: 28.07.2012 Сообщений: 2,813
	07.04.2017, 14:50	27
	Сообщение от DimKaKiber Но согласно https://ru.wikipedia.org/wiki/... 0%B8%D1%8F, "Цифровая обработка речевых сигналов" (Рабинер, Шафер) и ряду других талмудов по цифровой обработке сигналов нужно возвести в квадрат по отдельности вещественную и мнимую части результата прямого БПФ, сложить их и обратное преобразование произвести. А выложи код с реализацией через ДПФ, авось кто-нибудь поймет и поможет. Можно попробовать сравнить результаты с тем, что выдает какой-нибудь мат-пакет по типу Matlab и прочих. 0

MrGluck Форумчанин 8215 / 5045 / 1437 Регистрация: 29.11.2010 Сообщений: 13,453
	07.04.2017, 15:13	28
	в OpenCV были готовые функции для БПФ, сам их пользовал, довольно шустрые. Можно взять для сравнения результата. 0

@nonedark2008 1394 / 1023 / 325 Регистрация: 28.07.2012 Сообщений: 2,813
	10.04.2017, 23:53	32
	Сообщение от DimKaKiber Т.е. обратное БПФ от призведения спектра на комплексно-сопряженный спектр - все равно хрень получается Ну, в формуле проблем не было. Умножение числа на комплексно-сопряженное и дает в результате сумму квадратов мнимой и действительной частей. Я уже давно советовал проверить расчеты в какой-нибудь мат-библиотеке: Matlab, MathCAD, Maple... Может вообще оказаться, что у тебя не ответ, а одни сплошные ошибки округления. Какой из результатов вообще больше похож на правду? 0