Корреляция

Для определения степени похожести двух сигналов используется понятие корреляции.

Совокупность отсчетов сигнала можно считать вектором в N-мерном пространстве. Тогда можно определить корреляцию как угол между двумя векторами в N-мерном пространстве. Чем меньше угол, тем больше вектора похожи. На практике проще всего вычислять косинус угла. Вспоминая геометрию, получаем формулу:

Удобство этой формулы в том, что она даёт нормированный результат в диапазоне -1 ... +1 независимо от диапазона входных значений сигналов.

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

Часто проще использовать квадрат косинуса угла, чтобы сократить вычисления.

Если интересует сравнение только переменной составляющей сигналов, необходимо избавиться от постоянной.

Если поэкспериментировать с этой формулой, можно заметить, что сравнение происходит по характеру изменения функции сигналов. Это приводит к тому, что формула будет давать большую степень похожести для даже не очень похожих сигналов. Поэтому на практике, чтобы увеличить надёжность оценки, часто сравнивают не сигналы, а их производные, иногда даже вторые. Как вывод, для синусоидальных сигналов используются сами сигналы, так как их производные - такие же синусоидальные сигналы. Производные находятся простым вычитанием из текущего значения предшествующего. При этом появляется ограничение - сигнал должен быть динамичным. Использование производной также избавляет от необходимости учитывать постоянную составляющую.

Присмотревшись внимательно к этой формуле, можно понять, за счёт чего она дает нормированный результат. В числителе находится взаимная мощность сигналов. В знаменателе произведение мощности каждого из сигналов.

Прямая реализация несложна, однако требует большого объёма вычислений, если надо найти значение корреляционной функции между двумя сигналами, а не просто вычислить значение корреляции в одной точке. Сократить их объём можно воспользовавшись предположением, что если взять произвольное число точек для сравнения, то корреляция сигналов с их использованием будет примерно такая же, как и с использованием полного числа точек. Как вариант, можно делать тестовое предварительное сравнение по меньшему числу точек и если его результат превышает некоторую заданную величину, просчитывать корреляцию с использованием всех точек.

Другим способом увеличения скорости вычисления корреляционной функции является использование быстрого преобразования Фурье.

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

Реализация

/**

* алгоритм взят с http://alglib.sources.ru/fft/

*

*/

 /*************************************************************************

 Корелляция с использованием БПФ

 На входе:

    Signal     -   сигнал, с которым проводим корреляцию.

                    Нумерация элементов от 0 до SignalLen-1

    SignalLen   -   длина сигнала. Больше или равна длине образца

    Pattern     -   образец, корреляцию сигнала с которым мы ищем

                    Нумерация элементов от 0 до PatternLen-1

    PatternLen -   длина образца

 На выходе:

    Correlation -   квадрат корреляции в точках от 0 до SignalLen - PatternLen;

                    полный размер массива 2^int(0.5+Ln2(PatternLen + SignalLen))

 *************************************************************************/

 public static double[] fastCorrelation(final double[] signal,

                                                  int signallen,

                                                  final double[] pattern,

                                                  int patternlen)

 {

         int sumLen = signallen + patternlen;

         int power = 0;

         int i = 1;

         for( ; i < sumLen; power++ ) i = i << 1;

         sumLen = i;

         double[] a1 = new double[sumLen];

         double[] a2 = new double[sumLen];

         double Spower = 0.0;

         for( i = 0; i < patternlen; i++ ) {

                 a1[i] = signal[i];

                 Spower = Spower + signal[i] * signal[i];

         }

         for( ; i < signallen; i++ ) a1[i] = signal[i];

         //for( ; i < sumLen; i++ ) a1[i] = 0.0;//0: default value

         double Ppower = 0.0;

         for( i = 0; i < patternlen; i++ ) {

                 a2[i] = pattern[i];

                 Ppower = Ppower + pattern[i] * pattern[i];

         }

         //for( ; i < sumLen; i++ ) a2[i] = 0.0;//0: default value

         realFastFourierTransform( a1, power, false );

         realFastFourierTransform( a2, power, false );

         a1[0] = a1[0] * a2[0];

         a1[1] = a1[1] * a2[1];

         for( i = 2; i < sumLen; i += 2 ) {

                 double t1 = a1[i];

                 double t2 = a1[i + 1];

                 a1[i]     = t1 * a2[i] + t2 * a2[i + 1];

                 a1[i + 1] = t2 * a2[i] - t1 * a2[i + 1];

         }

         realFastFourierTransform( a1, power, true );

         a2[0] = 0.0;

         double denominator = power * Ppower;

         if( denominator > 0.0 ) {

                 a2[0] = (a1[0] * a1[0] / denominator);

                 if( a1[0] < 0.0 ) a2[0] = -a2[0];

         }

         i = 1;

         for( int p = 0; i < signallen - patternlen; i++, p++ ) {

                 Spower = Spower - signal[p] * signal[p];

                 Spower = Spower + signal[p + patternlen] * signal[p + patternlen];

                 a2[i] = 0.0;

                 denominator = Spower * Ppower;

                 if( denominator > 0.0 ) {

                         a2[i] = (a1[i] * a1[i] / denominator);

                         if( a1[i] < 0.0 ) a2[i] = -a2[i];

                 }

         }

         return a2;

 }