c++ - 以与英特尔性能原语相同的方式构建 MFCC 滤波器组

Question

我正在尝试构建用于生成 MFCC 的三角滤波器。我有基于 IPP 6 的现有代码，但随着 IPP 8 即将推出，我真的很想获得一个可以正常工作并且不依赖于现在不受支持的旧库的实现。

我已经生成了相关的梅尔缩放中心频率(加上两端的 2)。

然后我尝试按如下方式构建过滤器:

std::vector< std::vector< float > > ret;
int numFilters  = freqPositions.size() - 2;

for( int f = 1; f < numFilters + 1; f++ )
{
    float freqLow   = freqPositions[f - 1];
    float freqMid   = freqPositions[f];
    float freqHigh  = freqPositions[f + 1];

    float binLow    = (freqLow  / (sampleRate / 2)) * (numSamples + 1);
    float binMid    = (freqMid  / (sampleRate / 2)) * (numSamples + 1);
    float binHigh   = (freqHigh / (sampleRate / 2)) * (numSamples + 1);

    std::vector< float > fbank;
    for( int s = 0; s < (numSamples + 1); s++ )
    {
        if      ( s >= binLow && s < binMid )
        {
            const float fAmpl   = (s - binLow) / (float)(binMid - binLow);
            fbank.push_back( fAmpl );
        }
        else if ( s >= binMid && s <= binHigh )
        {
            const float fAmpl   = 1.0f - ((s - binMid) / (float)(binHigh - binMid));
            fbank.push_back( fAmpl );
        }
        else
        {
            fbank.push_back( 0.0f );
        }

    }

    ret.push_back( fbank );
}

然后我将上述 vector 与 FFT 结果分段相乘(其中 bin 0 是 0Hz 或 DC Offset bin)并将它们相加(本质上是一个点积)。

这似乎工作得相当好，但我得到的结果与 IPP 相比有很大不同，足以让我有点担心。

我做错了什么吗？

整个过程包括进行 FFT、计算返回的复 vector (std::abs) 的大小，然后应用按上述方法计算的滤波器组。代码如下:

std::vector< float > ApplyFilterBanks( std::vector< std::vector< float > >& filterBanks, std::vector< float >& fftMags )
{
    std::vector< float > ret;
    for( int fb = 0; fb < (int)filterBanks.size(); fb++ )
    {
        float res = 0.0f;
        Vec::Dot( res, &filterBanks[fb].front(), &fftMags.front(), filterBanks[fb].size() );
        ret.push_back( res );
    }
    return ret;
}

{
    const int kFFTSize      = 1 << mFFT.GetFFTOrder();
    const int kFFTSizeDiv2  = kFFTSize >> 1;
    std::vector< float > audioToFFT;
    audioToFFT.reserve( kFFTSize );
    std::copy( pAudio, pAudio + numSamples, std::back_inserter( audioToFFT ) );
    audioToFFT.resize( kFFTSize );

    std::vector< float > hammingWindow( numSamples );
    Vec::BuildHammingWindow( hammingWindow );
    Vec::Multiply( &audioToFFT.front(), &audioToFFT.front(), &hammingWindow.front(), numSamples );

    std::vector< std::complex< float > > fftResult( kFFTSize + 1 );

    // FFT the incoming audio.
    mFFT.ForwardFFT( &fftResult.front(), &audioToFFT.front(), kFFTSize );

    // Calculate the magnitudes of the resulting FFT.
    Vec::Magnitude( &audioToFFT.front(), &fftResult.front(), kFFTSizeDiv2 + 1 );
    //Vec::Multiply( &audioToFFT.front(), &audioToFFT.front(), &audioToFFT.front(), kFFTSizeDiv2 + 1 );

    // Apply the MFCC filter banks.
    std::vector< float > filtered   = ApplyFilterBanks( mFilterBanks, audioToFFT );
}

这是一个图，其中系列 1 是我的 MFCC，系列 2 是 IPP:

在日志和举重阶段(我已经确认其工作方式与 IPP 相同)之后，结果更加错误。

如有任何想法和建议，我们将不胜感激!

编辑:我应该指出这里有一些关于 IPP 函数的文档:

http://software.intel.com/sites/products/documentation/hpc/ipp/ipps/ipps_ch8/functn_MelFBankInitAlloc.html

这似乎显示了数学。但是，我不确定 yk 和 ck 到底是什么......

Answer 1

好的，我现在在这个问题上做得更好了。

我发现了两个问题，首先是:

float binLow    = (freqLow  / (sampleRate / 2)) * (numSamples + 1);
float binMid    = (freqMid  / (sampleRate / 2)) * (numSamples + 1);
float binHigh   = (freqHigh / (sampleRate / 2)) * (numSamples + 1);

应该是:

float binLow    = (freqLow  / (sampleRate / 2)) * (numSamples);
float binMid    = (freqMid  / (sampleRate / 2)) * (numSamples);
float binHigh   = (freqHigh / (sampleRate / 2)) * (numSamples);

其次，我错误地计算了通过梅尔空间的步数。我正在做以下事情:

const float melStep     = melDiff / (numFilterBanks + 2);

当我应该做的时候:

const float melStep     = melDiff / (numFilterBanks + 1);

现在我的结果虽然不完全相同，但现在显示出更好的对应关系:

以及最终的 MFCC: