FFTW3 вычисляет взаимную корреляцию в том же сигнале - PullRequest
Новая
       29

FFTW3 вычисляет взаимную корреляцию в том же сигнале

0 голосов
/ 13 июня 2018

В настоящее время я создаю код C, который принимает в качестве входного файла файл wav (в частности, только один канал исходного файла wav) и выполняет кратковременное преобразование Фурье.Основная часть кода выглядит следующим образом: 

stft_data = (fftw_complex*)(fftw_malloc(sizeof(fftw_complex)*windowSize));

fft_result= (fftw_complex*)(fftw_malloc(sizeof(fftw_complex)*windowSize));

storage = (fftw_complex*)(fftw_malloc(sizeof(fftw_complex)*storage_capacity));

//define the fftw plane
fftw_plan plan_forward;
plan_forward = fftw_plan_dft_1d(windowSize, stft_data, fft_result, FFTW_FORWARD, FFTW_ESTIMATE);

//integer indexes
int i,counter ;
counter = 0 ;
//create a Hamming window
double hamming_result[windowSize];
hamming(windowSize, hamming_result);

//implement the stft position indexes
int chunkPosition = 0; //actual chunk position
int readIndex ; //read the index of the wav file

while (chunkPosition < wav_length ){
    //read the window
    for(i=0; i<windowSize; i++){

      readIndex = chunkPosition + i;

      if (readIndex < wav_length){
        stft_data[i] = wav_data[readIndex]*hamming_result[i]*_Complex_I  + 0.0*I;
      }
      else{
        //if we are beyond the wav_length
        stft_data[i] = 0.0*_Complex_I + 0.0*I;//padding
        break;
      }
    }
    //compute the fft
    fftw_execute(plan_forward);
    //store the stft in a data structure
    for (i=0; i<windowSize;i++)
    {
      //printf("RE: %.2f  IM: %.2f\n", creal(fft_result[i]),cimag(fft_result[i]));
      storage[counter] = creal(fft_result[i]) + cimag(fft_result[i]);
      counter+=1;
    }

    //update indexes
    chunkPosition += hop_size;
    printf("Chunk Position %d\n", chunkPosition);
    printf("Counter position %d\n", counter);
    printf("Fourier transform done\n");

}

Как только БПФ вычислено для выбранного окна, я сохраняю действительную и мнимую части БПФ в переменной storage.

После этого я хотел бы вычислить взаимную корреляцию между точками данных в каждом из N окон, которые у меня есть в конце.В качестве примера я хотел бы вычислить корреляцию между первой точкой данных первого окна (storage[0]) и первым элементом второго окна (storage[windowSize+1]).Однако я сталкиваюсь с некоторыми проблемами, и у меня нет разумных ценностей.Согласно тому, что я изучал, корреляция в пространстве Фурье это просто комплексное умножение между двумя членами Фурье.Таким образом, то, что я делаю, выглядит примерно так: 

correlation = storage[0]*conj(storage[windowSize+1]);

Однако я получил очень большие значения, что заставляет меня задуматься, действительно ли я вычисляю корреляцию.

Где я не прав?Как мне масштабировать результаты корреляции?Как мне вычислить корреляцию со значениями Фурье?а затем, как я должен построить значения Фурье, которые я получаю из расчетов FFTW3?мне сдвинуть все значения или они уже сдвинуты?

Большое спасибо

1 Ответ

0 голосов
/ 13 июня 2018

Линия storage[counter] = creal(fft_result[i]) + cimag(fft_result[i]); делает хранилище чисто реальным.Поскольку вычисление correlation = storage[0]*conj(storage[windowSize+1]); является следующим шагом в вычислении взаимной корреляции, существует проблема.В самом деле, нет смысла спрягать действительное число.

Попытка storage[counter] = fft_result[i]; может частично решить проблему.Кроме того, correlation = storage[0]*conj(storage[windowSize+1]); следует изменить на correlation = storage[0]*conj(storage[windowSize]);

. При выполнении correlation = storage[0]*conj(storage[windowSize]); получается величина индекса [0] ДПФ корреляции.Действительно, storage[0] соответствует среднему значению первого кадра, тогда как storage[windowSize] соответствует среднему значению второго кадра.Он не равен средним, а намного больше, так как он масштабируется по длине кадра windowSize.

Чтобы вычислить корреляцию между двумя сигналами, следующим шагом должно быть: 

for (i=0; i<windowSize;i++)
{
    dftofcorrelation[i]=storage[i]*conj(storage[i+windowSize]
}

Затем обратный ДПФ должен быть применен к массиву dftofcorrelation, чтобы получить корреляцию в виде массива.Следует иметь в виду, что ни прямой, ни обратный ДПФ FFTW не включают в себя никакого масштабирования, см. , что на самом деле FFTW вычисляет : 

fftw_execute(plan_backward);

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

Общее масштабирование, вызванное FFTW, представляет собоймощность windowSize (размер окна ^ 3?).Это можно проверить путем вычисления автокорреляции однородного сигнала (который является равномерным).

...