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

В настоящее время я создаю код 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

Решение

Линия 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 (WindowSize ^ 3?). Это можно проверить, рассчитав автокорреляцию однородного сигнала (который является равномерным).

1

Другие решения

Других решений пока нет …