Как правильно рассчитать кросс-корреляцию?

Question

Я написал программу, которая (надеюсь) вычисляет взаимную корреляцию между двумя сигналами. Хотя я сделал хороший поиск того, как должен быть выполнен расчет, я не мог понять некоторые важные детали. Меня особенно беспокоит средний расчет. Похоже, что в некоторых алгоритмах используется среднее значение всего набора данных для выполнения расчета корреляции для каждого сдвига (или задержки). Другими словами, они используют постоянное среднее значение. Я даже нашел некоторые алгоритмы, которые вычисляют знаменатель только один раз, используя его как постоянное значение для остальных задержек. Однако я считаю, что как средние, так и знаменатели следует рассчитывать итеративно, учитывая только данные внутри диапазона суперпозиции. Поэтому я написал две версии этой программы. Они, похоже, дают очень похожие результаты при меньших задержках. Я хотел бы знать, какой из них правильный.

Итерационный средний:

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 

FILE *input1, *input2, *output; 
int m = 0, n = 0, t; 
float *VarA, *VarB, *Results, *Results2; 

void open_inputs_and_output(); 

void count_and_check_lines(); 

void allocate_memory(); 

void read_data(); 

void allocate_memory2(); 

void write_output(); 

int main() 
{ 
    float SumAverageA = 0, SumAverageB = 0, AverageA, AverageB, SubA, SubB, SumAB = 0, SumAs = 0, SumBs = 0, Correl; 
    int p = 0, i, delay; 

    open_inputs_and_output(); 

    count_and_check_lines(); 

    rewind(input1); 
    rewind(input2); 

    allocate_memory(); 

    read_data(); 

    fclose(input1); 
    fclose(input2); 

    printf("How many lag steps from the origin do you want to calculate?\nIf you want as many steps as the number of input points, type -1.\n"); 
    scanf("%i", &p); 

    if(p < -1) 
    { 
     printf("Bad number!\n"); 
     exit(1); 
    } 
    else if(p == -1) 
     t = n; 
    else 
     t = p; 

    allocate_memory2(); 

    printf("\nWait...\n\n"); 

    for(delay = 0; delay < t; delay++) 
    { 
     for(i = delay; i < n; i ++) 
     { 
      SumAverageA += VarA[i]; 
      SumAverageB += VarB[(i - delay)]; 
     } 

     AverageA = SumAverageA/(n - delay); 
     AverageB = SumAverageB/(n - delay); 

     for(i = delay; i < n; i++) 
     { 
      SubA = VarA[i] - AverageA; 
      SubB = VarB[(i - delay)] - AverageB; 
      SumAB += (SubA * SubB); 
      SumAs += (SubA * SubA); 
      SumBs += (SubB * SubB); 
     } 

     Correl = SumAB/(sqrt(SumAs * SumBs)); 

     Results[delay] = Correl; 

     SumAverageA = 0; 
     SumAverageB = 0; 
     SumAB = 0; 
     SumAs = 0; 
     SumBs = 0; 

     for(i = delay; i < n; i++) 
     { 
      SubB = VarB[i] - AverageB; 
      SubA = VarA[(i - delay)] - AverageA; 
      SumAB += (SubA * SubB); 
      SumAs += (SubA * SubA); 
      SumBs += (SubB * SubB); 
     } 

     Correl = SumAB/(sqrt(SumAs * SumBs)); 

     Results2[delay] = Correl; 

     SumAverageA = 0; 
     SumAverageB = 0; 
     SumAB = 0; 
     SumAs = 0; 
     SumBs = 0; 
    } 

    printf("Calculations performed.\n"); 

    free(VarA); 
    free(VarB); 

    write_output(); 

    free(Results); 
    free(Results2); 

    fclose(output); 

    return 0; 
} 

void open_inputs_and_output() 
{ 
    input1 = fopen("C:\\...\\test.txt","r"); 

    if (input1 == NULL) 
    { 
     printf("Error! Could not open input 1.\n"); 
     exit(1); 
    } 
    else 
     printf("Input1 opening: OK.\n"); 

    input2 = fopen("C:\\...\\test2.txt","r"); 

    if (input2 == NULL) 
    { 
     printf("Error! Could not open input 2.\n"); 
     exit(1); 
    } 
    else 
     printf("Input2 opening: OK.\n"); 

    output = fopen("C:\\...\\out.txt","w"); 

    if (output == NULL) 
    { 
     printf("Error! Could not open output.\n"); 
     exit(1); 
    } 
    else 
     printf("Output opening: OK.\n"); 
} 

void count_and_check_lines() 
{ 
    float k; 

    while(fscanf(input1,"%f",&k) == 1) 
     n++; 

    printf("n = %i\n", n); 

    while(fscanf(input2,"%f",&k) == 1) 
     m++; 

    printf("m = %i\n", m); 

    if(m != n) 
    { 
     printf("Error: Number of rows does not match!\n"); 
     exit(1); 
    } 
    else 
     printf("Number of rows matches.\n"); 
} 

void allocate_memory() 
{ 
    VarA = calloc(n, sizeof(float)); 

    if(VarA == NULL) 
    { 
     printf("Could not allocate memory for VarA.\n"); 
     exit(1); 
    } 

    VarB = calloc(m, sizeof(float)); 

    if(VarA == NULL) 
    { 
     printf("Could not allocate memory for VarB.\n"); 
     exit(1); 
    } 
} 

void read_data() 
{ 
    int i; 

    for(i = 0; i < n; i++) 
     fscanf(input1,"%f",&VarA[i]); 

    printf("Data A successfully read.\n"); 

    for(i = 0; i < m; i++) 
     fscanf(input2,"%f",&VarB[i]); 

    printf("Data B successfully read.\n"); 
} 

void allocate_memory2() 
{ 
    Results = calloc(t, sizeof(float)); 

    if(Results == NULL) 
    { 
     printf("Could not allocate memory for Results.\n"); 
     exit(1); 
    } 

    Results2 = calloc(t, sizeof(float)); 

    if(Results2 == NULL) 
    { 
     printf("Could not allocate memory for Results2.\n"); 
     exit(1); 
    } 
} 

void write_output() 
{ 
    int i; 

    for(i = t - 1; i > 0; i--) 
     fprintf(output,"-%i %f\n", i , Results2[i]); 

    for(i = 0; i < t; i++) 
     fprintf(output,"%i %f\n", i , Results[i]); 

    printf("Results written.\n"); 
} 

Constant средний:

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 

FILE *input1, *input2, *output; 
int m = 0, n = 0, t; 
float *VarA, *VarB, *Results, *Results2; 

void open_inputs_and_output(); 

void count_and_check_lines(); 

void allocate_memory(); 

void read_data(); 

void allocate_memory2(); 

void write_output(); 

int main() 
{ 
    float SumAverageA = 0, SumAverageB = 0, AverageA, AverageB, SubA, SubB, SumAB = 0, SumAs = 0, SumBs = 0, Correl; 
    int p = 0, i, delay; 

    open_inputs_and_output(); 

    count_and_check_lines(); 

    rewind(input1); 
    rewind(input2); 

    allocate_memory(); 

    read_data(); 

    fclose(input1); 
    fclose(input2); 

    printf("How many lag steps from the origin do you want to calculate?\nIf you want as many steps as the number of input points, type -1.\n"); 
    scanf("%i", &p); 

    if(p < -1) 
    { 
     printf("Bad number!\n"); 
     exit(1); 
    } 
    else if(p == -1) 
     t = n; 
    else 
     t = p; 

    allocate_memory2(); 

    printf("\nWait...\n\n"); 

    for(i = 0; i < n; i ++) 
    { 
     SumAverageA += VarA[i]; 
     SumAverageB += VarB[i]; 
    } 

    AverageA = SumAverageA/n; 
    AverageB = SumAverageB/n; 

    for(delay = 0; delay < t; delay++) 
    { 
     for(i = delay; i < n; i++) 
     { 
      SubA = VarA[i] - AverageA; 
      SubB = VarB[(i - delay)] - AverageB; 
      SumAB += (SubA * SubB); 
      SumAs += (SubA * SubA); 
      SumBs += (SubB * SubB); 
     } 

     Correl = SumAB/(sqrt(SumAs * SumBs)); 

     Results[delay] = Correl; 

     SumAverageA = 0; 
     SumAverageB = 0; 
     SumAB = 0; 
     SumAs = 0; 
     SumBs = 0; 

     for(i = delay; i < n; i++) 
     { 
      SubB = VarB[i] - AverageB; 
      SubA = VarA[(i - delay)] - AverageA; 
      SumAB += (SubA * SubB); 
      SumAs += (SubA * SubA); 
      SumBs += (SubB * SubB); 
     } 

     Correl = SumAB/(sqrt(SumAs * SumBs)); 

     Results2[delay] = Correl; 

     SumAverageA = 0; 
     SumAverageB = 0; 
     SumAB = 0; 
     SumAs = 0; 
     SumBs = 0; 
    } 

    printf("Calculations performed.\n"); 

    free(VarA); 
    free(VarB); 

    write_output(); 

    free(Results); 
    free(Results2); 

    fclose(output); 

    return 0; 
} 

void open_inputs_and_output() 
{ 
    input1 = fopen("C:\\...\\test.txt","r"); 

    if (input1 == NULL) 
    { 
     printf("Error! Could not open input 1.\n"); 
     exit(1); 
    } 
    else 
     printf("Input1 opening: OK.\n"); 

    input2 = fopen("C:\\...\\test2.txt","r"); 

    if (input2 == NULL) 
    { 
     printf("Error! Could not open input 2.\n"); 
     exit(1); 
    } 
    else 
     printf("Input2 opening: OK.\n"); 

    output = fopen("C:\\...\\out.txt","w"); 

    if (output == NULL) 
    { 
     printf("Error! Could not open output.\n"); 
     exit(1); 
    } 
    else 
     printf("Output opening: OK.\n"); 
} 

void count_and_check_lines() 
{ 
    float k; 

    while(fscanf(input1,"%f",&k) == 1) 
     n++; 

    printf("n = %i\n", n); 

    while(fscanf(input2,"%f",&k) == 1) 
     m++; 

    printf("m = %i\n", m); 

    if(m != n) 
    { 
     printf("Error: Number of rows does not match!\n"); 
     exit(1); 
    } 
    else 
     printf("Number of rows matches.\n"); 
} 

void allocate_memory() 
{ 
    VarA = calloc(n, sizeof(float)); 

    if(VarA == NULL) 
    { 
     printf("Could not allocate memory for VarA.\n"); 
     exit(1); 
    } 

    VarB = calloc(m, sizeof(float)); 

    if(VarA == NULL) 
    { 
     printf("Could not allocate memory for VarB.\n"); 
     exit(1); 
    } 
} 

void read_data() 
{ 
    int i; 

    for(i = 0; i < n; i++) 
     fscanf(input1,"%f",&VarA[i]); 

    printf("Data A successfully read.\n"); 

    for(i = 0; i < m; i++) 
     fscanf(input2,"%f",&VarB[i]); 

    printf("Data B successfully read.\n"); 
} 

void allocate_memory2() 
{ 
    Results = calloc(t, sizeof(float)); 

    if(Results == NULL) 
    { 
     printf("Could not allocate memory for Results.\n"); 
     exit(1); 
    } 

    Results2 = calloc(t, sizeof(float)); 

    if(Results2 == NULL) 
    { 
     printf("Could not allocate memory for Results2.\n"); 
     exit(1); 
    } 
} 

void write_output() 
{ 
    int i; 

    for(i = t - 1; i > 0; i--) 
     fprintf(output,"-%i %f\n", i , Results2[i]); 

    for(i = 0; i < t; i++) 
     fprintf(output,"%i %f\n", i , Results[i]); 

    printf("Results written.\n"); 
} 

Ссылки:

http://www.jot.fm/issues/issue_2010_03/column2.pdf

http://paulbourke.net/miscellaneous/correlate/

Answer 1

Если ваши процессы wide sense stationary, то средние значения со временем не будут меняться. Если процессы также являются ergodic, то эти средние значения могут быть получены путем вычисления среднего значения одного временного ряда. В этом случае вы можете использовать все доступные образцы, чтобы получить как можно более точную оценку среднего. Это, естественно, приводит к вашей «постоянной усредненной» реализации.

Если, с другой стороны, ваши процессы не являются широкозначимыми и устойчивыми, то получение хорошей оценки их местных средств может оказаться более сложной задачей. Оценка среднего значения в меньшем временном окне, где процессы являются приблизительно стационарными, может быть разумным подходом (который аналогичен вашей реализации «Итеративно-средний»). Обратите внимание, что существуют другие подходы, и их пригодность будет зависеть от конкретного приложения и свойств конкретных сигналов.

Возможно, вам будет интересно узнать, как вы узнаете, являются ли ваши процессы широкими, неподвижными. К сожалению, если вы не знаете много о том, как эти процессы были сгенерированы, часто лучшее, что вы могли бы сделать, это работать в соответствии с гипотезой о том, что процессы являются широкими, а затем пытаться опровергнуть эту гипотезу (наблюдая результаты, выходящие за пределы ожидаемых диапазонов , см. statistical hypothesis testing).