Реализация свертки в C++ с использованием fftw 3

Question

ОБНОВИТЬ

Смотрите мой основополагающий вопрос о DSP stackexchange здесь

ОБНОВИТЬ

Я все еще испытываю потрескивание на выходе. Эти трещины теперь менее выражены и слышны только при увеличении громкости

ОБНОВИТЬ

Следуя приведенному здесь совету, я удалил потрескивающий звук из моего вывода. Я протестирую с другими доступными HRIR, чтобы убедиться, что свёртка действительно работает правильно, и отвечу на этот вопрос, как только я проверил, что мой код теперь работает

ОБНОВИТЬ

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

#define HRIR_LENGTH 512
#define WAV_SAMPLE_SIZE 256

    while (signal_input_wav.read(&signal_input_buffer[0], WAV_SAMPLE_SIZE) >= WAV_SAMPLE_SIZE)
    {
#ifdef SKIP_CONVOLUTION
        // Copy the input buffer over
        std::copy(signal_input_buffer.begin(),
                  signal_input_buffer.begin() + WAV_SAMPLE_SIZE,
                  signal_output_buffer.begin());

        signal_output_wav.write(&signal_output_buffer[0], WAV_SAMPLE_SIZE);
#else
        // Copy the first segment into the buffer
        // with zero padding
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            if (i < WAV_SAMPLE_SIZE)
            {
                signal_buffer_fft_in[i] = signal_input_buffer[i];
            }
            else
            {
                signal_buffer_fft_in[i] = 0; // zero pad
            }
        }

        // Dft of the signal segment
        fftw_execute(signal_fft);

        // Convolve in the frequency domain by multiplying filter kernel with dft signal
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            signal_buffer_ifft_in[i] = signal_buffer_fft_out[i] * left_hrir_fft_out[i]
                - signal_buffer_fft_out[HRIR_LENGTH - i] * left_hrir_fft_out[HRIR_LENGTH - i];

            signal_buffer_ifft_in[HRIR_LENGTH - i] = signal_buffer_fft_out[i] * left_hrir_fft_out[HRIR_LENGTH - i]
                + signal_buffer_fft_out[HRIR_LENGTH - i] * left_hrir_fft_out[i];

            //double re = signal_buffer_out[i];
            //double im = signal_buffer_out[BLOCK_OUTPUT_SIZE - i];
        }

        // inverse dft back to time domain
        fftw_execute(signal_ifft);

        // Normalize the data
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            signal_buffer_ifft_out[i] = signal_buffer_ifft_out[i] / HRIR_LENGTH;
        }

        // Overlap-add method
        for (int i = 0; i < HRIR_LENGTH; ++i)
        {
            if (i < WAV_SAMPLE_SIZE)
            {
                signal_output_buffer[i] = signal_overlap_buffer[i] + signal_buffer_ifft_out[i];
            }
            else
            {
                signal_output_buffer[i] = signal_buffer_ifft_out[i];
                signal_overlap_buffer[i] = signal_output_buffer[i]; // record into the overlap buffer
            }
        }

        // Write the block to the output file
        signal_output_wav.write(&signal_output_buffer[0], HRIR_LENGTH);

#endif
    }

Результирующий звуковой файл содержит треск; предположительно, артефакты остались от ошибочной реализации fftw. Кроме того, запись блоков 512 (HRIR_LENGTH), по-видимому, приводит к некоторому наложению псевдонимов, при этом звуковой файл при воспроизведении звучит как виниловая пластинка, замедленная. Запись блоков размером WAV_SAMPLE_SIZE (256, половина fft-вывода) воспроизводится с нормальной скоростью. Тем не менее, независимо от этого треск остается.

ОРИГИНАЛ

Я пытаюсь реализовать свертку с использованием библиотеки fftw в C++. Я могу прекрасно загрузить свой фильтр и заполнить нулями как фильтр (длиной 512), так и входной сигнал (длиной 513), чтобы получить блок вывода сигнала 1024 и использовать его в качестве размера БПФ.

Вот мой код:

#define BLOCK_OUTPUT_SIZE 1024
#define HRIR_LENGTH 512

#define WAV_SAMPLE_SIZE 513
#define INPUT_SHIFT 511

while (signal_input_wav.read(&signal_input_buffer[0], WAV_SAMPLE_SIZE) >= WAV_SAMPLE_SIZE)
{
#ifdef SKIP_CONVOLUTION
    // Copy the input buffer over
    std::copy(signal_input_buffer.begin(),
              signal_input_buffer.begin() + WAV_SAMPLE_SIZE,
              signal_output_buffer.begin());

    signal_output_wav.write(&signal_output_buffer[0], WAV_SAMPLE_SIZE);
#else
    // Zero pad input
    for (int i = 0; i < INPUT_SHIFT; ++i)
        signal_input_buffer[WAV_SAMPLE_SIZE + i] = 0;

    // Copy to the signal convolve buffer
    for (int i = 0; i < BLOCK_OUTPUT_SIZE; ++i)
    {
        signal_buffer_in[i] = signal_input_buffer[i];
    }

    // Dft of the signal segment
    fftw_execute(signal_fft);

    // Convolve in the frequency domain by multiplying filter kernel with dft signal
    for (int i = 1; i < BLOCK_OUTPUT_SIZE; ++i)
    {
        signal_buffer_out[i] = signal_buffer_in[i] * left_hrir_fft_in[i]
            - signal_buffer_in[BLOCK_OUTPUT_SIZE - i] * left_hrir_fft_in[BLOCK_OUTPUT_SIZE - i];

        signal_buffer_out[BLOCK_OUTPUT_SIZE - i]
            = signal_buffer_in[BLOCK_OUTPUT_SIZE - i] * left_hrir_fft_in[i]
                    + signal_buffer_in[i] * left_hrir_fft_in[BLOCK_OUTPUT_SIZE - i];

        double re = signal_buffer_out[i];
        double im = signal_buffer_out[BLOCK_OUTPUT_SIZE - i];
    }

    // inverse dft back to time domain
    fftw_execute(signal_ifft);

    // Normalize the data
    for (int i = 0; i < BLOCK_OUTPUT_SIZE; ++i)
    {
        signal_buffer_out[i] = signal_buffer_out[i] / i;
    }

    // Overlap and add with the previous block
    if (first_block)
    {
        first_block = !first_block;
        for (int i = 0; i < BLOCK_OUTPUT_SIZE; ++i)
        {
            signal_output_buffer[i] = signal_buffer_out[i];
        }
    }
    else
    {
        for (int i = WAV_SAMPLE_SIZE; i < BLOCK_OUTPUT_SIZE; ++i)
        {
            signal_output_buffer[i] = signal_output_buffer[i] + signal_buffer_out[i];
        }
    }

    // Write the block to the output file
    signal_output_wav.write(&signal_output_buffer[0], BLOCK_OUTPUT_SIZE);
#endif
}

В конце концов, полученный выходной файл содержит мусор, но не все нули. Вещи, которые я пробовал:

1) Использование стандартного комплексного интерфейса fftw_plan_dft_1d с соответствующим типом fftw_complex. Те же проблемы возникают.

2) Использование меньшего размера входной выборки и итерация по заполненным нулями блокам (overlap-add).

Я также отмечаю, что это не ошибка libsndfile; переключение SKIP_CONVOLUTION успешно приводит к копированию входного файла в выходной файл.