waveIn 音频信号的 FFT

• 下载演示项目 - 33.3 KB

引言

快速傅里叶变换 (FFT) 允许用户查看音频信号的频谱内容。 这里展示的 FFT 代码由 Don Cross 编写，他的主页似乎已被撤下。 我将尝试解释 FFT 的实用性，因为它与音频信号相关，而不是解释 FFT 的数学理论。

FFT 允许用户获得音频信号的频谱构成，获得其各种频率的分贝，或者获得其各种频率的强度。 频谱查看器（如上图所示）、均衡器或 VU 表都可以使用 FFT 来显示其结果。 那么它们之间的区别取决于几个等式中的一个，这些等式采用 FFT 的实部和虚部，并返回要在图形结果中使用的强度或分贝级别。 以下代码采用 FFT 结果的实部和虚部，并返回强度和分贝。

inline double GetFrequencyIntensity(double re, double im)
{
    return sqrt((re*re)+(im*im));
}
#define mag_sqrd(re,im) (re*re+im*im)
#define Decibels(re,im) ((re == 0 && im == 0) ? (0) : 
        10.0 * log10(double(mag_sqrd(re,im))))
#define Amplitude(re,im,len) (GetFrequencyIntensity(re,im)/(len))
#define AmplitudeScaled(re,im,len,scale) ((int)Amplitude(re,im,len)%scale)

FFT 使用音频信号作为它的实部，并使用一个 NULL 指针作为它的虚部，表明虚部数据不存在。 在返回时，FFT 将返回基于作为实部给出的数据的实部和虚部数据分量。 这与返回的样本镜像，因此 0-FFT_LEN/2 包含数据，而 FFT_LEN/2 到 FFT_LEN 包含数据的反转。 这个错误在我的代码中得到了纠正。 执行 FFT 的代码如下

    DWORD nCount = 0;
    for (DWORD dw = 0; dw < FFT_LEN; dw++)
    {
        {
            //copy audio signal to fft real component for left channel
            finleft[nCount] = (double)((short*)pwh->lpData)[dw++];
            //copy audio signal to fft real component for right channel
            finright[nCount++] = (double)((short*)pwh->lpData)[dw];
        }
    }
    //Perform FFT on left channel
    fft_double(FFT_LEN/2,0,finleft,NULL,fout,foutimg);
    float re,im,fmax=-99999.9f,fmin=99999.9f;
    for(int i=1;i < FFT_LEN/4-1;i++)
    //Use FFT_LEN/4 since the data is mirrored within the array.
    {
        re = fout[i];
        im = foutimg[i];
        //get amplitude and scale to 0..256 range
        fdraw[i]=AmplitudeScaled(re,im,FFT_LEN/2,256);
        if (fdraw[i] > fmax)
        {
            fmax = fdraw[i];
        }
        if (fdraw[i] < fmin)
        {
            fmin = fdraw[i];
        }
    }
    //Use this to send the average band amplitude to something
    int nAvg, nBars=16, nCur = 0;
    for(int i=1;i < FFT_LEN/4;i++)
    {
        nAvg = 0;
        for (int n=0; n < nBars; n++)
        {
            nAvg += (int)fdraw[i];
        }
        nAvg /= nBars;
        //Send data here to something,
        //nothing to send it to so we print it.
        TRACE("Average for Bar#%d is %d\n",nCur++,nAvg);
        i+=nBars-1;
    }
    DataHolder* pDataHolder = (DataHolder*)lpData;
    // Draw left channel
    CFrequencyGraph* pPeak = (CFrequencyGraph*)pDataHolder->pData;
    if (::IsWindow(pPeak->GetSafeHwnd()))
    {
        pPeak->SetYRange(0,256);
        pPeak->Update(FFT_LEN/4,fdraw);
    }
    
    // Perform FFT on right channel
    fmax=-99999.9f,fmin=99999.9f;
    fft_double(FFT_LEN/2,0,finright,NULL,fout,foutimg);
    fdraw[0] = fdraw[FFT_LEN/4] = 0;
    for(i=1;i < FFT_LEN/4-1;i++)
    //Use FFT_LEN/4 since the data is mirrored within the array.
    {
        re = fout[i];
        im = foutimg[i];
        //get Decibels in 0-110 range
        fdraw[i] = Decibels(re,im);
        if (fdraw[i] > fmax)
        {
            fmax = fdraw[i];
        }
        if (fdraw[i] < fmin)
        {
            fmin = fdraw[i];
        }
    }
    //Draw right channel
    CFrequencyGraph* pPeak2 = (CFrequencyGraph*)pDataHolder->pData2;
    if (::IsWindow(pPeak2->GetSafeHwnd()))
    {
        pPeak2->SetNumberOfSteps(50);
        //Use updated dynamic range for scaling
        pPeak2->SetYRange((int)fmin,(int)fmax);
        pPeak2->Update(FFT_LEN/4,fdraw);
    }

此代码包含在每次 waveIn 函数返回更新的音频信号数据时调用的回调函数中。 实际执行 FFT 的代码如下所示

void fft_double (unsigned int p_nSamples, bool p_bInverseTransform, 
    double *p_lpRealIn, double *p_lpImagIn, 
    double *p_lpRealOut, double *p_lpImagOut)
{

    if(!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;


    unsigned int NumBits;
    unsigned int i, j, k, n;
    unsigned int BlockSize, BlockEnd;

    double angle_numerator = 2.0 * PI;
    double tr, ti;

    if( !IsPowerOfTwo(p_nSamples) )
    {
        return;
    }

    if( p_bInverseTransform ) angle_numerator = -angle_numerator;

    NumBits = NumberOfBitsNeeded ( p_nSamples );


    for( i=0; i < p_nSamples; i++ )
    {
        j = ReverseBits ( i, NumBits );
        p_lpRealOut[j] = p_lpRealIn[i];
        p_lpImagOut[j] = (p_lpImagIn == NULL) ? 0.0 : p_lpImagIn[i];
    }


    BlockEnd = 1;
    for( BlockSize = 2; BlockSize <= p_nSamples; BlockSize <<= 1 )
    {
        double delta_angle = angle_numerator / (double)BlockSize;
        double sm2 = sin ( -2 * delta_angle );
        double sm1 = sin ( -delta_angle );
        double cm2 = cos ( -2 * delta_angle );
        double cm1 = cos ( -delta_angle );
        double w = 2 * cm1;
        double ar[3], ai[3];

        for( i=0; i < p_nSamples; i += BlockSize )
        {

            ar[2] = cm2;
            ar[1] = cm1;

            ai[2] = sm2;
            ai[1] = sm1;

            for ( j=i, n=0; n < BlockEnd; j++, n++ )
            {

                ar[0] = w*ar[1] - ar[2];
                ar[2] = ar[1];
                ar[1] = ar[0];

                ai[0] = w*ai[1] - ai[2];
                ai[2] = ai[1];
                ai[1] = ai[0];

                k = j + BlockEnd;
                tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
                ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];

                p_lpRealOut[k] = p_lpRealOut[j] - tr;
                p_lpImagOut[k] = p_lpImagOut[j] - ti;

                p_lpRealOut[j] += tr;
                p_lpImagOut[j] += ti;

            }
        }

        BlockEnd = BlockSize;

    }


    if( p_bInverseTransform )
    {
        double denom = (double)p_nSamples;

        for ( i=0; i < p_nSamples; i++ )
        {
            p_lpRealOut[i] /= denom;
            p_lpImagOut[i] /= denom;
        }
    }

}

它需要以下支持函数

///////////////////////////////////////////////////////////
// check is a number is a power of 2
///////////////////////////////////////////////////////////

bool IsPowerOfTwo( unsigned int p_nX )
{

    if( p_nX < 2 ) return false;

    if( p_nX & (p_nX-1) ) return false;

    return true;

}


///////////////////////////////////////////////////////////
// return needed bits for fft
///////////////////////////////////////////////////////////

unsigned int NumberOfBitsNeeded( unsigned int p_nSamples )
{

    int i;

    if( p_nSamples < 2 )
    {
        return 0;
    }

    for ( i=0; ; i++ )
    {
        if( p_nSamples & (1 << i) ) return i;
    }

}



///////////////////////////////////////////////////////////
// ?
///////////////////////////////////////////////////////////

unsigned int ReverseBits(unsigned int p_nIndex, unsigned int p_nBits)
{

    unsigned int i, rev;

    for(i=rev=0; i < p_nBits; i++)
    {
        rev = (rev << 1) | (p_nIndex & 1);
        p_nIndex >>= 1;
    }

    return rev;
}



///////////////////////////////////////////////////////////
// return a frequency from the basefreq and num of samples
///////////////////////////////////////////////////////////

double Index_to_frequency(unsigned int p_nBaseFreq, 
    unsigned int p_nSamples, unsigned int p_nIndex)
{

    if(p_nIndex >= p_nSamples)
    {
        return 0.0;
    }
    else if(p_nIndex <= p_nSamples/2)
    {
        return ( (double)p_nIndex / 
                 (double)p_nSamples * p_nBaseFreq );
    }
    else
    {
        return ( -(double)(p_nSamples-p_nIndex) / 
                  (double)p_nSamples * p_nBaseFreq );
    }

}

包含的示例类 CFrequencyGraph 将使用频率强度绘制一个示例均衡器、峰值表和频谱图。 希望这可以作为对音频信号的快速傅里叶变换的用途和基础知识的体面介绍。 FFT 的其他功能包括将其与 节拍检测算法 结合使用，以检测音频信号中的节拍。 另一个包含有用的 FFT 信息的页面位于 FFT 频谱分析仪。

本文使用 waveIn* 函数从播放源检索声卡的数据。 因此，您必须手动打开声卡属性，并更改录音源以使用单声道/立体声混音或 waveIn 作为选定的录音线路。 此外，如果您将播放 wave 选项设置为 false，这将覆盖录音选项，并且不会出现声音。 有关执行反向 FFT 以将频率转换为音频信号的信息，请在 Google 上搜索“反向 FFT”或“IFFT”。

历史

• 版本 1.3：修复了频谱绘制，更改了 Pixelgram 颜色，更改了 Process 例程以匹配新的录音参数。
• 版本 1.2：修复了许多绘图错误，更改为使用幅度缩放和分贝。