为什么理想带通滤波器不能按预期工作？

2018-06-11 19:29:53

这是最近的版本，产生的效果接近所需的

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#6224 - 感谢@SleuthEye以供参考）。
那么你能否建议我基于FFT的使用简单的解决方案，这将允许我至少删除给定范围的frequneces ？

我正试图实现理想的带通滤波器。但它没有像我预期的那样工作 - 只有高频被削减。

这是我的实现描述：

从采样速率为8000hz的PCM（原始）16位格式读取到尺寸为1024的短路缓冲器的放大值

应用FFT从时域到频域

零全部频率<300和> 3700：

逆FFT

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开始）

我做错了什么？

这是之前的声音频谱在这里输入图像描述

在将所需值清零后，这里是声音频率在这里输入图像描述

请帮忙！

更新：

这里是我想出的代码，我应该用Zeroes填充什么？

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;
}
}

你可能想看看这个答案来解释你正在观察的效果。

否则，由于频域中的矩形函数（具有零转换和无限阻带衰减）对应于时间上的无限长度脉冲响应，因此您尝试实现的“理想”滤波器比实际实现更多的是数学工具域。

要获得更实用的滤波器，您必须首先根据您的具体应用需求定义所需的滤波器特性，例如转换宽度和阻带衰减。基于这些规格，滤波器系数可以使用各种滤波器设计方法之一来导出，例如：

窗口方法

频率采样方法

Parks-McClellan方法

双线性变换在模拟模板上的应用

...

也许最接近你正在做的是Window方法。使用这种方法，像三角形窗口这样简单的东西可以帮助增加阻带衰减，但是您可能想要尝试其他窗口选择（许多可从相同链接获得）。增加窗口长度将有助于减小转换宽度。

一旦完成了滤波器设计，就可以使用overlap-add方法或overlap-save方法在频域中应用滤波器。使用这两种方法之一，你可以将输入信号分成长度为L的块，然后填充到一些方便的大小N> = L + M-1，其中M是滤波器系数的数量（例如，如果你有一个滤波器42个系数，你可以选择N = 128，其中L = N-M + 1 = 87）。

使用FFT / IFFT进行快速卷积滤波需要零填充至少为滤波器长度的两倍（出于性能原因，通常为2的下幂），然后使用重叠添加或重叠保存方法来消除循环卷积伪像。

在做完真正的FFT之后，您可以获得两次光谱数据：一次位于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实现自动执行此步骤。

顺便说一句，只是清理频率箱像你做的不是一个好办法做这样的过滤器。期望一个非常讨厌的相位响应和大量的信号振铃。

链接地址: http://www.djcxy.com/p/33841.html

上一篇: Why Ideal band pass filter not working as expected?

下一篇: Is my understanding of FFT and Pitch Detection correct here?