c# - 使用 DirectShow .NET 从 SampleGrabber 计算 FFT

Question

我正在使用 DirectShow .NET 开发一个项目。我正在尝试集成一个名为“WPF Sound Visualization Library”的库，它创建了一个频谱分析仪视觉对象。

为了让视觉效果正常工作，我需要在我的播放器中实现这两种方法:

• GetFFTData(float[] fftDataBuffer) - 将当前 FFT 数据分配给缓冲区。

  备注:缓冲区中的 FFT 数据应仅包含实数强度值。这意味着如果您的 FFT 算法返回复数(很多人都这样做)，您将运行类似于以下的算法:for(int i = 0; i < complexNumbers.Length/2; i++) fftResult[i] = Math.Sqrt (complexNumbers[i].Real * complexNumbers[i].Real + complexNumbers[i].Imaginary * complexNumbers[i].Imaginary);

• GetFFTFrequencyIndex(int frequency) - 获取给定频率的 FFT 数据缓冲区中的索引。

编辑:我已经添加了 SampleGrabber 并将其回调与 GetFFTData 集成(仍未测试)。但是如何集成GetFFTFrequencyIndex方法呢？

    protected int SampleCB(double SampleTime, IMediaSample pSample)
    {
        IntPtr pointer = IntPtr.Zero;
        pSample.GetPointer(out pointer);
        sampleDataBytes = new byte[pSample.GetSize()];
        Marshal.Copy(pointer, sampleDataBytes, 0, sampleDataBytes.Length);
        var sampleTime = SampleTime;
        var actualDataLength = pSample.GetActualDataLength();

        /* Release unmanaged resources */
        Marshal.ReleaseComObject(pSample);
        pSample = null;

        return (int)HResults.S_OK;
    }

    #region ISpectrumPlayer
    byte[] sampleDataBytes = null;

    public bool GetFFTData(float[] fftDataBuffer)
    {
        if (sampleDataBytes != null)
        {
            var sampleData = Utils.GetInt16Array(sampleDataBytes);
            double[] pRealIn = new double[sampleData.Length];

            for (var i = 0; i <= sampleData.Length - 1; i++)
                pRealIn[i] = sampleData[i];

            var pImagIn = new double[sampleDataBytes.Length];
            var pRealOut = new double[sampleDataBytes.Length];
            var pImagOut = new double[sampleDataBytes.Length];

            FFTUtils.Compute((uint) pRealIn.Length, pRealIn, pImagIn, pRealOut, pImagOut, false);

            fftDataBuffer = new float[sampleDataBytes.Length];

            for (int i = 0; i < pRealOut.Length; i++)
                fftDataBuffer[i] = (float) Math.Sqrt(pRealOut[i] * pRealOut[i] + pImagOut[i] * pImagOut[i]);
        }

        return true;
    }

    public int GetFFTFrequencyIndex(int frequency)
    {
        throw new NotImplementedException();
    }
    #endregion

我对这个类有帮助的方法:

public class FFTUtils
{
    public const Double DDC_PI = 3.14159265358979323846;

    /// <summary>
    /// Verifies a number is a power of two
    /// </summary>
    /// <param name="x">Number to check</param>
    /// <returns>true if number is a power two (i.e.:1,2,4,8,16,...)</returns>
    public static Boolean IsPowerOfTwo(UInt32 x)
    {
        return ((x != 0) && (x & (x - 1)) == 0);
    }

    /// <summary>
    /// Get Next power of number.
    /// </summary>
    /// <param name="x">Number to check</param>
    /// <returns>A power of two number</returns>
    public static UInt32 NextPowerOfTwo(UInt32 x)
    {
        x = x - 1;
        x = x | (x >> 1);
        x = x | (x >> 2);
        x = x | (x >> 4);
        x = x | (x >> 8);
        x = x | (x >> 16);
        return x + 1;
    }

    /// <summary>
    /// Get Number of bits needed for a power of two
    /// </summary>
    /// <param name="PowerOfTwo">Power of two number</param>
    /// <returns>Number of bits</returns>
    public static UInt32 NumberOfBitsNeeded(UInt32 PowerOfTwo)
    {
        if (PowerOfTwo > 0)
        {
            for (UInt32 i = 0, mask = 1; ; i++, mask <<= 1)
            {
                if ((PowerOfTwo & mask) != 0)
                    return i;
            }
        }
        return 0; // error
    }

    /// <summary>
    /// Reverse bits
    /// </summary>
    /// <param name="index">Bits</param>
    /// <param name="NumBits">Number of bits to reverse</param>
    /// <returns>Reverse Bits</returns>
    public static UInt32 ReverseBits(UInt32 index, UInt32 NumBits)
    {
        UInt32 i, rev;

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

        return rev;
    }

    /// <summary>
    /// Return index to frequency based on number of samples
    /// </summary>
    /// <param name="Index">sample index</param>
    /// <param name="NumSamples">number of samples</param>
    /// <returns>Frequency index range</returns>
    public static Double IndexToFrequency(UInt32 Index, UInt32 NumSamples)
    {
        if (Index >= NumSamples)
            return 0.0;
        else if (Index <= NumSamples / 2)
            return (double)Index / (double)NumSamples;

        return -(double)(NumSamples - Index) / (double)NumSamples;
    }

    /// <summary>
    /// Compute FFT
    /// </summary>
    /// <param name="NumSamples">NumSamples Number of samples (must be power two)</param>
    /// <param name="pRealIn">Real samples</param>
    /// <param name="pImagIn">Imaginary (optional, may be null)</param>
    /// <param name="pRealOut">Real coefficient output</param>
    /// <param name="pImagOut">Imaginary coefficient output</param>
    /// <param name="bInverseTransform">bInverseTransform when true, compute Inverse FFT</param>
    public static void Compute(UInt32 NumSamples, Double[] pRealIn, Double[] pImagIn,
                        Double[] pRealOut, Double[] pImagOut, Boolean bInverseTransform)
    {
        UInt32 NumBits;    /* Number of bits needed to store indices */
        UInt32 i, j, k, n;
        UInt32 BlockSize, BlockEnd;

        double angle_numerator = 2.0 * DDC_PI;
        double tr, ti;     /* temp real, temp imaginary */

        if (pRealIn == null || pRealOut == null || pImagOut == null)
        {
            // error
            throw new ArgumentNullException("Null argument");
        }
        if (!IsPowerOfTwo(NumSamples))
        {
            // error
            throw new ArgumentException("Number of samples must be power of 2");
        }
        if (pRealIn.Length < NumSamples || (pImagIn != null && pImagIn.Length < NumSamples) ||
             pRealOut.Length < NumSamples || pImagOut.Length < NumSamples)
        {
            // error
            throw new ArgumentException("Invalid Array argument detected");
        }

        if (bInverseTransform)
            angle_numerator = -angle_numerator;

        NumBits = NumberOfBitsNeeded(NumSamples);

        /*
        **   Do simultaneous data copy and bit-reversal ordering into outputs...
        */
        for (i = 0; i < NumSamples; i++)
        {
            j = ReverseBits(i, NumBits);
            pRealOut[j] = pRealIn[i];
            pImagOut[j] = (double)((pImagIn == null) ? 0.0 : pImagIn[i]);
        }

        /*
        **   Do the FFT itself...
        */
        BlockEnd = 1;
        for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1)
        {
            double delta_angle = angle_numerator / (double)BlockSize;
            double sm2 = Math.Sin(-2 * delta_angle);
            double sm1 = Math.Sin(-delta_angle);
            double cm2 = Math.Cos(-2 * delta_angle);
            double cm1 = Math.Cos(-delta_angle);
            double w = 2 * cm1;
            double ar0, ar1, ar2;
            double ai0, ai1, ai2;

            for (i = 0; i < NumSamples; i += BlockSize)
            {
                ar2 = cm2;
                ar1 = cm1;

                ai2 = sm2;
                ai1 = sm1;

                for (j = i, n = 0; n < BlockEnd; j++, n++)
                {
                    ar0 = w * ar1 - ar2;
                    ar2 = ar1;
                    ar1 = ar0;

                    ai0 = w * ai1 - ai2;
                    ai2 = ai1;
                    ai1 = ai0;

                    k = j + BlockEnd;
                    tr = ar0 * pRealOut[k] - ai0 * pImagOut[k];
                    ti = ar0 * pImagOut[k] + ai0 * pRealOut[k];

                    pRealOut[k] = (pRealOut[j] - tr);
                    pImagOut[k] = (pImagOut[j] - ti);

                    pRealOut[j] += (tr);
                    pImagOut[j] += (ti);
                }
            }

            BlockEnd = BlockSize;
        }

        /*
        **   Need to normalize if inverse transform...
        */
        if (bInverseTransform)
        {
            double denom = (double)(NumSamples);

            for (i = 0; i < NumSamples; i++)
            {
                pRealOut[i] /= denom;
                pImagOut[i] /= denom;
            }
        }
    }

    /// <summary>
    /// Calculate normal (power spectrum)
    /// </summary>
    /// <param name="NumSamples">Number of sample</param>
    /// <param name="pReal">Real coefficient buffer</param>
    /// <param name="pImag">Imaginary coefficient buffer</param>
    /// <param name="pAmpl">Working buffer to hold amplitude Xps(m) = | X(m)^2 | = Xreal(m)^2  + Ximag(m)^2</param>
    public static void Norm(UInt32 NumSamples, Double[] pReal, Double[] pImag, Double[] pAmpl)
    {
        if (pReal == null || pImag == null || pAmpl == null)
        {
            // error
            throw new ArgumentNullException("pReal,pImag,pAmpl");
        }
        if (pReal.Length < NumSamples || pImag.Length < NumSamples || pAmpl.Length < NumSamples)
        {
            // error
            throw new ArgumentException("Invalid Array argument detected");
        }

        // Calculate amplitude values in the buffer provided
        for (UInt32 i = 0; i < NumSamples; i++)
        {
            pAmpl[i] = pReal[i] * pReal[i] + pImag[i] * pImag[i];
        }
    }

    /// <summary>
    /// Find Peak frequency in Hz
    /// </summary>
    /// <param name="NumSamples">Number of samples</param>
    /// <param name="pAmpl">Current amplitude</param>
    /// <param name="samplingRate">Sampling rate in samples/second (Hz)</param>
    /// <param name="index">Frequency index</param>
    /// <returns>Peak frequency in Hz</returns>
    public static Double PeakFrequency(UInt32 NumSamples, Double[] pAmpl, Double samplingRate, ref UInt32 index)
    {
        UInt32 N = NumSamples >> 1;   // number of positive frequencies. (numSamples/2)

        if (pAmpl == null)
        {
            // error
            throw new ArgumentNullException("pAmpl");
        }
        if (pAmpl.Length < NumSamples)
        {
            // error
            throw new ArgumentException("Invalid Array argument detected");
        }

        double maxAmpl = -1.0;
        double peakFreq = -1.0;
        index = 0;

        for (UInt32 i = 0; i < N; i++)
        {
            if (pAmpl[i] > maxAmpl)
            {
                maxAmpl = (double)pAmpl[i];
                index = i;
                peakFreq = (double)(i);
            }
        }

        return samplingRate * peakFreq / (double)(NumSamples);
    }
}

非常感谢!

Answer 1

如果我没记错的话，该算法采用实数(例如 int[n])信号或复数信号(例如 int[n][2])和返回复数 FFT 结果。

所以，这看起来很简单:

您获取输入值(您可以在图表中绘制为时间值，例如左扬声器音频值)并将它们输入 pRealIn 参数。在 pImagIn 中放置零(与在 pRealIn 中一样多)。在 bInverseTransform 中你输入 false(当然)。

然后您将把结果取回pRealOut & pImagOut。结果缓冲区在逻辑上应与输入缓冲区大小相同。您必须采用这两个输出缓冲区并将它们像这样成对组合(对于 OUT 数组的每个元素):

fftDataBuffer[k] = Math.Sqrt(pRealOut[k] * pRealOut[k] + pImagOut[k] * pImagOut[k]); // Do this from 1 to n

FFT 结果是一组复数值(x = 实部，y = 虚部 - 您可以在笛卡尔系统中将其描述为向量)。你想要向量的大小/振幅，这就是你执行上述操作的原因。

那是为了 GetFFTData。

---

我看到您有一个名为 IndexToFrequency 的函数。所以这可能有效。您所要做的就是为缓冲区的每个索引调用此方法。即:

for(int i=0; i<n; i++) freq[i] = IndexToFrequency(i, n);

保存这些值，然后在您的 GetFFTFrequencyIndex(int frequency) 中找到输入参数 (frequency) 与 freq 中的元素最接近的匹配项[n] 并返回其索引。

我认为这就足够了。

重要提示:确保您的缓冲区大小为二次方(NextPowerOfTwo 似乎旨在帮助您做到这一点)。

如果您的数据有时碰巧较小，您可以在末尾用零填充值(也就是将零附加到输入缓冲区)。另外:为了获得更好的分辨率，您可以再次用零填充数据。这将增加您可能需要的结果的“平滑度”。

如果我可能会问(只是好奇 :))，您在哪里找到这段代码？

所以，就是这样! :)