Вот на одном уважаемом сайте нашел такой исходник для вычисления быстрой корреляции через БПФ. Там в API многие операторы перегружены, так что будем считать, что синтаксис правильный... Вопрос по теории.
Известно, что корреляция вычисляется также как свертка, но для получения ядра необходимо развернуть паттерн-сигнал вокруг нулевой точки.
Смущает, отсутсвие этого преобразования в помеченном месте.
Это ошибка или Быстрая корреляция вычисляется именно так, т.е без преобразования в ядро свертки? И почему?
Код:
/*************************************************************************
Корелляция с использованием БПФ
На входе:
Signal - сигнал, с которым проводим корелляцию.
Нумерация элементов от 0 до SignalLen-1
SignalLen - длина сигнала.
Pattern - образец, корелляцию сигнала с которым мы ищем
Нумерация элементов от 0 до PatternLen-1
PatternLen - длина образца
На выходе:
Signal - значения корелляции в точках от 0 до
SignalLen-1.
*************************************************************************/
void fastcorrelation(ap::real_1d_array& signal,
int signallen,
const ap::real_1d_array& pattern,
int patternlen)
{
ap::real_1d_array a1;
ap::real_1d_array a2;
int nl;
int i;
double t1;
double t2;
// начиная отсюда
nl = signallen+patternlen;
i = 1;
while(i<nl)
{
i = i*2;
}
nl = i;
//до сюда, ИМХО, тоже ошибка, т.к. nl должно быть >= signallen+patternlen-1
//Например, если signal = 3, pattern=2, nl = signallen+patternlen-1 = 4 что является
// допустимым значением. Но данный алгоритм найдет 8, а не 4
// А интервал должен быть от 0 до nl-1 , т.е. от 0 до 3
// Или я не прав ???
a1.setbounds(0, nl-1);
a2.setbounds(0, nl-1);
for(i = 0; i <= signallen-1; i++)
{
a1(i) = signal(i);
}
for(i = signallen; i <= nl-1; i++)
{
a1(i) = 0;
}
for(i = 0; i <= patternlen-1; i++)
{
a2(i) = pattern(i);
}
for(i = patternlen; i <= nl-1; i++)
{
a2(i) = 0;
}
//
вот тут, по идее, должно быть преобразование паттерн-сигнала в ядро свертки,
или сразу надо было переворачивать при копировании.
//
realfastfouriertransform(a1, nl, false);
realfastfouriertransform(a2, nl, false);
a1(0) = a1(0)*a2(0);
a1(1) = a1(1)*a2(1);
for(i = 1; i <= nl/2-1; i++)
{
t1 = a1(2*i);
t2 = a1(2*i+1);
a1(2*i) = t1*a2(2*i)+t2*a2(2*i+1);
a1(2*i+1) = t2*a2(2*i)-t1*a2(2*i+1);
}
realfastfouriertransform(a1, nl, true);
for(i = 0; i <= signallen-1; i++)
{
signal(i) = a1(i);
}
}