c++ - 为什么理想带通滤波器不能按预期工作?

标签 c++ filtering signal-processing fft frequency-analysis

这是最新版本,效果接近预期

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

enter image description here

但为什么会这样呢???

重要的是我刚开始学习DSP,所以我可能不知道一些重要的想法 (我对此深表歉意,但我有需要解决的任务:我需要减少录音机语音中的背景噪音,我尝试通过切断 <300 && > 3700 范围内的录制语音频率(如人声在 [300;3700] 范围内)我从那个方法开始,因为它很简单,但我发现 out - 它无法应用(请参阅 - https://dsp.stackexchange.com/questions/6220/why-is-it-a-bad-idea-to-filter-by-zeroing-out-fft-bins/6224#6224 - 感谢@SleuthEye 提供引用)。
那么,您能否根据 FFT 的使用建议我一个简单的解决方案,至少允许我删除给定的频率范围

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

这是我的实现描述:

  1. 从采样率为 8000 赫兹的 PCM(原始)16 位格式读取振幅值到大小为 1024 的短裤缓冲区
  2. 应用 FFT 从时域到频域
  3. 将 < 300 和 > 3700 的所有频率归零:
  4. 逆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 开始)

我做错了什么?

这是enter image description here之前的声音频谱

这是我将所需值归零后的声音频率 enter image description here

请帮忙!

更新:

这是我想出的代码,我应该用零填充什么???

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

它产生以下频谱(低频被削减,但高频没有) enter image description here

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

最佳答案

使用 FFT/IFFT 的快速卷积滤波需要零填充到滤波器长度的至少两倍(出于性能原因通常是 2 的下一次幂),然后使用重叠添加或重叠保存方法来消除循环卷积伪像.

关于c++ - 为什么理想带通滤波器不能按预期工作?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/24101814/

相关文章:

c - 定点: out_Q15 = antilog( in_Q25 ) 计算,避免溢出?

c++ - C++ 中的多重声明

c++ - 我想通过windows api制作一个时钟(GUI程序),但是时钟上的文字没有改变?

java - Jtable过滤失败

javascript - 多个键的 JSON 过滤运行缓慢

python - 我可以将使用 librosa 生成的频谱图转换回音频吗?

c++ - 构造函数中抽象类的默认值

c++ - 无法捕获模板特化方法抛出的异常

PHP 在插入 MySQL 时向字符添加斜杠

过滤/归一化不良信号的算法