Почему идеальный полосовой фильтр не работает должным образом?

Вот последняя версия, которая производит эффект, близкий к желаемому

void DeleteFrequencies(short *audioDataBuffer, const int bufferSize, int lowestFrequency, int highestFrequency, int sampleRate )
{
    int frequencyInHzPerSample = sampleRate / bufferSize;
    /*             __________________________
    /* ___________                           __________________________  filter kernel   */
    int nOfPointsInFilterKernel = (lowestFrequency / frequencyInHzPerSample) + ( bufferSize - highestFrequency / frequencyInHzPerSample);
    U u;
    double *RealX = new  double[bufferSize];
    double *ImmX = new  double[bufferSize];
    ShortArrayToDoubleArray(audioDataBuffer, RealX, bufferSize);

    // padd with zeroes, so that inputSignalSamplesNumber + kernelLength - 1 = bufferSize

    // convert to frequency domain
    ForwardRealFFT(RealX, ImmX, bufferSize);
    // cut frequences < 300 && > 3400
    int Multiplyer = 1;
    for (int i = 0; i < 512; ++i)
    {
        if (i * 8000 / 1024 > 3400 || i * 8000 / bufferSize < 300 )
        {
            RealX[i] = 0;
            ImmX[i] = 0;
        }
        if (i < lowestFrequency / frequencyInHzPerSample || i > highestFrequency / frequencyInHzPerSample )
            Multiplyer = 0;
        else 
            Multiplyer = 1;
        RealX[i] = RealX[i] * Multiplyer /*ReH[f]*/  - ImmX[i] * Multiplyer;
        ImmX[i] = ImmX[i] * Multiplyer + RealX[i] * Multiplyer;

    }
    ReverseRealFFT(RealX, ImmX, bufferSize);
    DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);
    delete [] RealX;
    delete [] ImmX;
}

а почему так работает???

Важно то, что я только начал изучать DSP, поэтому я могу не знать о некоторых важных идеях (я извиняюсь за это, но у меня есть задача, которую мне нужно решить: мне нужно уменьшить фоновый шум в речи диктофона, я пытаюсь подойти к этому с помощью отрезав от записанных речевых частот в диапазонах <300 && > 3700 (как человеческий голос в диапазоне [300;3700]), я начал с этого метода, поскольку он прост, но я обнаружил - его нельзя применять (см. - https://dsp.stackexchange.com/questions/6220/why-is-it-a-bad-idea-to-filter-by-zeroing-out-fft-bins/6224 - спасибо @SleuthEye за справку).
Не могли бы вы предложить мне простое решение, основанное на использовании БПФ, которое позволит мне хотя бы удалить заданные диапазоны частот?

Я пытаюсь реализовать идеальный полосовой фильтр. Но это не работает, как я ожидаю - только высокие частоты обрезаются.

Вот мое описание реализации:

  1. Считать значения амплитуды из 16-битного формата PCM (необработанного) с частотой дискретизации 8000 Гц в буфер шорт размером 1024
  2. Применить БПФ для перехода от временной области к частотной области
  3. Обнулить все частоты <300 и> 3700:
  4. Обратное БПФ
union U
{
    char ch[2];
    short sh;
};
std::fstream in;
std::fstream out;
short audioDataBuffer[1024];
in.open ("mySound.pcm", std::ios::in | std::ios::binary);
out.open("mySoundFilteres.pcm", std::ios::out | std::ios::binary);
int i = 0;
bool isDataInBuffer = true;
U u;
while (in.good())
{
    int j = 0;
    for (int i = 0; i < 1024 * 2; i+=2)
    {
        if (false == in.good() && j < 1024) // padd with zeroes
        {
            audioDataBuffer[j] = 0;
        }
        in.read((char*)&audioDataBuffer[j], 2);
        cout << audioDataBuffer[j];
        ++j;
    }
    // Algorithm
    double RealX [1024] = {0};
    double ImmX [1024] = {0};
    ShortArrayToDoubleArray(audioDataBuffer, RealX, 1024);

    // convert to frequency domain
    ForwardRealFFT(RealX, ImmX, 1024);
    // cut frequences < 300 && > 3400
    for (int i = 0; i < 512; ++i)
    {
        if (i * 8000 / 1024 > 3400 || i * 8000 / 1024 < 300 )
        {
            RealX[i] = 0;
            ImmX[i] = 0;
        }
    }
    ReverseRealFFT(RealX, ImmX, 1024);
    DoubleArrayToShortArray(RealX, audioDataBuffer, 1024);
    for (int i = 0; i < 1024; ++i) // 7 6 5 4 3 2 1 0 - byte order hence we write ch[1]  then ch[0]
    {
        u.sh = audioDataBuffer[i];
        out.write(&u.ch[1], 1);
        out.write(&u.ch[0], 1);
    }
}
in.close();
out.close();

когда я записываю результат в файл, открываю его на храбрость и проверяю анализ спектра, и вижу, что высокие частоты обрезаются, но низкие остаются (они начинаются с 0)

Что я делаю не так?

Вот звуковая частота спектра перед

Вот частота звука после того, как я обнулил необходимые значения

Пожалуйста помоги!

Обновить:

Вот код, который я придумал, что я должен дополнить нулями???

void DeleteFrequencies(short *audioDataBuffer, const int bufferSize, int lowestFrequency, int highestFrequency, int sampleRate )
{
    // FFT must be the same length as output segment - to prevent circular convultion
    // 
    int frequencyInHzPerSample = sampleRate / bufferSize;
    /*             __________________________
    /* ___________                           __________________________  filter kernel   */
    int nOfPointsInFilterKernel = (lowestFrequency / frequencyInHzPerSample) + ( bufferSize - highestFrequency / frequencyInHzPerSample);
    U u;
    double *RealX = new  double[bufferSize];
    double *ImmX = new  double[bufferSize];
    ShortArrayToDoubleArray(audioDataBuffer, RealX, bufferSize);

    // padd with zeroes, so that inputSignalSamplesNumber + kernelLength - 1 = bufferSize

    // convert to frequency domain
    ForwardRealFFT(RealX, ImmX, bufferSize);
    // cut frequences < 300 && > 3400
    int Multiplyer = 1;
    for (int i = 0; i < 512; ++i)
    {
        /*if (i * 8000 / 1024 > 3400 || i * 8000 / bufferSize < 300 )
        {
            RealX[i] = 0;
            ImmX[i] = 0;
        }*/
        if (i < lowestFrequency / frequencyInHzPerSample || i > highestFrequency / frequencyInHzPerSample )
            Multiplyer = 0;
        else 
            Multiplyer = 1;
        RealX[i] = RealX[i] * Multiplyer /*ReH[f]*/  - ImmX[i] * Multiplyer;
        ImmX[i] = ImmX[i] * Multiplyer + RealX[i] * Multiplyer;

    }
    ReverseRealFFT(RealX, ImmX, bufferSize);
    DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);
    delete [] RealX;
    delete [] ImmX;
}

это производит следующий спектр (низкие частоты сокращены, но высокие нет)

void ForwardRealFFT(double* RealX, double* ImmX, int nOfSamples)
{

short nh, i, j, nMinus1, nDiv2, nDiv4Minus1, im, ip, ip2, ipm, nOfCompositionSteps, LE, LE2, jm1;
double ur, ui, sr, si, tr, ti;

// Step 1 : separate even from odd points
nh = nOfSamples / 2 - 1; 
for (i = 0; i <= nh; ++i)
{
    RealX[i] = RealX[2*i];
    ImmX[i] = RealX[2*i + 1];
}
// Step 2: calculate nOfSamples/2 points using complex FFT
// advantage in efficiency, as nOfSamples/2 requires 1/2 of the time as nOfSamples point FFT
nOfSamples /= 2;
ForwardDiscreteFT(RealX, ImmX, nOfSamples );
nOfSamples *= 2;

// Step 3: even/odd frequency domain decomposition
nMinus1 = nOfSamples - 1; 
nDiv2 = nOfSamples / 2;
nDiv4Minus1 = nOfSamples / 4 - 1;
for (i = 1; i <= nDiv4Minus1; ++i)
{
    im = nDiv2 - i;
    ip2 = i + nDiv2;
    ipm = im + nDiv2;
    RealX[ip2] = (ImmX[i] + ImmX[im]) / 2;
    RealX[ipm] = RealX[ip2];
    ImmX[ip2] = -(RealX[i] - RealX[im]) / 2;
    ImmX[ipm] = - ImmX[ip2];
    RealX[i] = (RealX[i] + RealX[im]) / 2;
    RealX[im] = RealX[i];
    ImmX[i] = (ImmX[i] - ImmX[im]) / 2;
    ImmX[im] = - ImmX[i];
}
RealX[nOfSamples * 3 / 4] = ImmX[nOfSamples / 4];
RealX[nDiv2] = ImmX[0];
ImmX[nOfSamples * 3 / 4] = 0;
ImmX[nDiv2] = 0;
ImmX[nOfSamples / 4] = 0;
ImmX[0] = 0;
// 3-rd step: combine the nOfSamples frequency spectra in the exact reverse order
// that the time domain decomposition took place
nOfCompositionSteps = log((double)nOfSamples) / log(2.0);
LE = pow(2.0,nOfCompositionSteps);
LE2 = LE / 2;
ur = 1;
ui = 0;
sr = cos(M_PI/LE2);
si = -sin(M_PI/LE2);
for (j = 1; j <= LE2; ++j)
{
    jm1 = j - 1;
    for (i = jm1; i <= nMinus1; i += LE)
    {
        ip = i + LE2;
        tr = RealX[ip] * ur - ImmX[ip] * ui;
        ti = RealX[ip] * ui + ImmX[ip] * ur;
        RealX[ip] = RealX[i] - tr;
        ImmX[ip] = ImmX[i] - ti;
        RealX[i] = RealX[i] + tr;
        ImmX[i] = ImmX[i] + ti;
    }
    tr = ur;
    ur = tr * sr - ui * si;
    ui = tr * si + ui * sr;
}
}
Источник

3 ответа

Решение

Возможно, вы захотите взглянуть на этот ответ для объяснения эффектов, которые вы наблюдаете.

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

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

Возможно, самым близким к тому, что вы делаете, является метод Window. Используя этот метод, простое, например, треугольное окно может помочь увеличить затухание в полосе задержания, но вы можете поэкспериментировать с другими вариантами выбора окна (многие из которых доступны по той же ссылке). Увеличение длины окна поможет уменьшить ширину перехода.

После того, как вы завершили проектирование фильтра, вы можете применить фильтр в частотной области, используя метод overlap-add или метод overlap-save. Используя любой из этих методов, вы бы разбили свой входной сигнал на куски длины L и добавили к некоторому удобному размеру N >= L+M-1, где M - число коэффициентов фильтра (например, если у вас есть фильтр с 42 коэффициента, вы можете выбрать N = 128, из которых L = N-M+1 = 87).

Источник

Быстрая сверточная фильтрация с использованием FFT/IFFT требует заполнения нулями, по меньшей мере, вдвое большей длины фильтра (и обычно до следующей степени 2 по соображениям производительности), а затем с использованием методов добавления с наложением или с перекрытием для удаления артефактов циклической свертки.

Источник

После выполнения реального БПФ вы получаете ваши спектральные данные дважды: один раз в ячейках с 0 по 512 и зеркальный спектр в ячейках с 513 по 1024. Однако ваш код очищает только нижний спектр.

Попробуй это:

for (int i = 0; i < 512; ++i)
{
    if (i * 8000 / 1024 > 3400 || i * 8000 / 1024 < 300 )
    {
        RealX[i] = 0;
        ImmX[i] = 0;

        // clear mirror spectrum as well:
        RealX[1023-i] = 0;
        ImmX[1023-i]  = 0;
    }
}

Это может помочь, если ваша реализация FFT не сделает этот шаг автоматически.

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

Источник
Другие вопросы по тегам