Нужно быстро и часто вычислять свертку. Использую код приведенный ниже.
Есть ли возможность ускорить вычисления?
Может кто-нибудь сумеет это записать с помощью стандартных алгоритмов?

class Exception {};

typedef std::vector<double> data;

void ConvolutionImpl(const data& shortSample, const data& longSample, size_t offset, data& result);

void Convolution(const data& sampleOne, const data& sampleTwo, size_t offset, data& result)
{
    if (sampleOne.size() <= sampleTwo.size())
        return ConvolutionImpl(sampleOne, sampleTwo, offset, result);
    else
        return ConvolutionImpl(sampleTwo, sampleOne, offset, result);
}

void ConvolutionImpl(const data& shortSample, const data& longSample, size_t offset, data& result)
{
    if (shortSample.empty())
        throw Exception();

    const size_t shortSize = shortSample.size();
    const size_t longSize = longSample.size();
    result = data(longSize, .0);

    for (size_t i = 0; i < longSize; ++i)
    {
        const size_t resultIndex = (i + longSize - offset) % longSize;
        for (size_t j = 0; j < shortSize; ++j)
        {
            const size_t longIndex = (i + longSize - j) % longSize;                    //Эта строка съедает ~15% времени
            result[resultIndex] += shortSample[j] * longSample[longIndex];    //Эта строка съедает ~80% времени
        }
    }
}

Используется так:

int main()
{
        //...

        data x;
        data y;

        //...

        Convolution(x, y, x.size()/2, result);

        //...
}

Здравствуйте, Аноним, Вы писали:

А>Нужно быстро и часто вычислять свертку. Использую код приведенный ниже.
А>Есть ли возможность ускорить вычисления?

Можно не вычислять (i+longSize-j)%longSize,
а сделать

resultIndex = (i<offset) ? longSize+i-offset : i-offset;
longIndex   = (i<j)      ? longSize+i-j      : i-j;

А>Может кто-нибудь сумеет это записать с помощью стандартных алгоритмов?

std::inner_product, std::transform, реверсные итераторы.

Здравствуйте, Кодт, Вы писали:

А>>Есть ли возможность ускорить вычисления?

К>Можно не вычислять (i+longSize-j)%longSize,
К>а сделать
К>

К>resultIndex = (i<offset) ? longSize+i-offset : i-offset;
К>longIndex   = (i<j)      ? longSize+i-j      : i-j;
К>

Спасибо, помогло — эта строка стала быстрее на 25%.

А>>Может кто-нибудь сумеет это записать с помощью стандартных алгоритмов?

К>std::inner_product, std::transform, реверсные итераторы.
Попробовал, но читабельности это не добавило. Плюнул и оставил как есть.

Здравствуйте, Аноним, Вы писали:

А>Нужно быстро и часто вычислять свертку. Использую код приведенный ниже.
А>Есть ли возможность ускорить вычисления?

Ага, использовать FFT Выигрыш в N/log2(N) раз c точностью до мультипликативной констаныты.

Здравствуйте, Аноним, Вы писали:

А>Нужно быстро и часто вычислять свертку. Использую код приведенный ниже.
А>Есть ли возможность ускорить вычисления?

угу...

void ConvolutionImpl(const data& shortSample, const data& longSample, size_t offset, data& result)
{
    if (shortSample.empty())
        throw Exception();

    const size_t shortSize = shortSample.size();
    const size_t longSize = longSample.size();

    //result = data(longSize, .0);
    
    result.resize(0);    // при частом вызове и повторном использовании result сэкономишь на его переаллокации
    result.resize(longSize);

    // доступ по индексу - не самая быстрая операция в STL-векторе, избавляемся по возможности.
    //for (size_t i = 0; i < longSize; ++i)
    data::iterator longBegin = result.begin();
    data::iterator resultBegin = result.begin();
    data::iterator resultEnd = result.end();
    data::iterator resultIndex = resultBegin - offset;

    for(data::const_iterator i = longSample.begin(), ie = longSample.end(); i!=ie; ++i, ++resultIndex)
    {
        if (resultIndex<resultBegin)
            resultIndex += longSize;

        data::const_iterator longIndex = i;
        double sum = 0.0;

        for (data::const_iterator j = shortSample.begin(), je = shortSample.end(); j!=je; ++j, --longIndex)
        {
            if (longIndex<longBegin)
                longIndex += longSize;

             sum += *j * *longIndex;    
        }

        *resultIndex = sum;
    }
}

а вообще, сделай все внутренние вычисления над массивами в памяти, чтобы не извращаться с итераторами,
т.е. процедура верхнего уровня берет вектор, готовит вектор результата и вызывает внутреннюю inline процедуру, ориентированную на обычные С++ массивы, подашь просто адрес первого элемента вектора туда. Это будет самым быстрым решением.

1) Самая торомозная операция у тебя это остаток от деления. дополни массив циклически до степени двойки и используй битовое И вместо остатка.
2) перепиши все на указатели
3) если сможешь на ассемблере применить SSE2 то ускоришь это дело раза в 2-3. фильтр дополни нулями до кратного двойки или четверки
4) можешь взять интеловский компиллер он сам вставит эти инструкции если просечет

Здравствуйте, Андрей Ушаков, Вы писали:

АУ>Здравствуйте, Аноним, Вы писали:

А>>Нужно быстро и часто вычислять свертку. Использую код приведенный ниже.
А>>Есть ли возможность ускорить вычисления?

АУ>Ага, использовать FFT Выигрыш в N/log2(N) раз c точностью до мультипликативной констаныты.

Да ну?

Допустим исходная послед имеет длину N(степень двойки), а фильтр имеет постоянную длин. Тогда вычисление в лоб имеет порядок O(N) -- ну что то вроде с*N*2 операций, не важно.

Вычисление через FFT.
1) найти образ входной последовательности O(N*log(N))
2) найти образ последовательности длины N, первые с членов которой состоят из коэф. фильтра, остальные 0. N*log(N), хотя можно тут уменьшить
3) перемножение последовательностей O(N)
4) обратное преобразование O(N*log(N))
и того O(N) против O(N*log(N))
где же выигрыш?

Здравствуйте, PoM-PoM 40mm, Вы писали:

АУ>>Ага, использовать FFT Выигрыш в N/log2(N) раз c точностью до мультипликативной константы.

PP4>Да ну?

Если длина одной из последовательностей M (много) меньше длины второй N (на этот момент в исходном собщении я внимания действительно не обратил), то выигрыш от FFT будет O(M/log2(N)). Т.е. он все при любом фиксированном отношении M/N при росте N, начиная с некоторого момента, FFT будет выгоднее "непосредственного" решения, т.к. M/log2(N) = (M/N) * (N/log2(N)). Т.е. моя предыдущая оценка остается в силе, ибо я предусмотрительно "заложился" на неизвестный коэффициент

Здравствуйте, PoM-PoM 40mm, Вы писали:

PP4>где же выигрыш?
Согласен.

Плюс борьба с эффектом Гиббса. Если бы FFT было панацеей, то про FIR фильтры я бы ничего не знал.

Всем спасибо.

Здравствуйте, PoM-PoM 40mm, Вы писали:

PP4>3) если сможешь на ассемблере применить SSE2 то ускоришь это дело раза в 2-3. фильтр дополни нулями до кратного двойки или четверки
PP4>4) можешь взять интеловский компиллер он сам вставит эти инструкции если просечет

если делать на floate, то и VC7.1 просечет, если ему в опциях "посоветовать".

Здравствуйте, vdimas, Вы писали:

V>Здравствуйте, PoM-PoM 40mm, Вы писали:

PP4>>3) если сможешь на ассемблере применить SSE2 то ускоришь это дело раза в 2-3. фильтр дополни нулями до кратного двойки или четверки
PP4>>4) можешь взять интеловский компиллер он сам вставит эти инструкции если просечет

V>если делать на floate, то и VC7.1 просечет, если ему в опциях "посоветовать".

Замечательно. А почему он на double не просекает? Там ведь от SSE2 выигрышь более чем в 2 раза (за счет лучшей системы команд)? Вообще я не в курсе, может уже сделали, но пора бы давно всю эту функцию на SSE/SSE2 написать (а так же еще полезные функции при обработке сигналов) и поставлять с компиляторами под x86.

Однажду у меня полином Лагранжа на флоатах падал, я теперь только двойную точность использую

Здравствуйте, PoM-PoM 40mm, Вы писали:

V>>если делать на floate, то и VC7.1 просечет, если ему в опциях "посоветовать".

PP4>Замечательно. А почему он на double не просекает? Там ведь от SSE2 выигрышь более чем в 2 раза (за счет лучшей системы команд)? Вообще я не в курсе, может уже сделали, но пора бы давно всю эту функцию на SSE/SSE2 написать (а так же еще полезные функции при обработке сигналов) и поставлять с компиляторами под x86.

а разве эти команды работают с double?