如何将频率仓与傅里叶变换幅度对齐?

How can I align the frequency bins with the fourier transform magnitude?

我正在尝试在 STM32F411RE Nucleo 开发板上实现具有相关复数幅度函数的快速傅里叶变换。我的目标是将具有多个正弦元素的组合信号分离成具有正确振幅的独立频率分量。

我的问题是我无法正确地将复数幅度函数的频点结果与频率对齐。我也开始质疑这些结果的有效性。

我已经尝试使用人们发布的许多不同的 FFT 算法实现,其中最著名的是 and Blog by LB9MG.

中列出的示例

据我所知,我有一个类似的方法,但不知何故,他们的方法产生了预期的结果,而我的却没有。下面是我修改过的代码,通过 SleuthEye 创建的实现来工作。

int main(void)
{
    fftLen = 32;    // can be 32, 64, 128, 256, 512, 1024, 2048, 4096
    half_fftLen = fftLen/2;
    volatile float32_t sampleFreq = 50 * fftLen;    // Fs = binsize * fft length, desired binsize = 50 hz

    arm_rfft_fast_instance_f32 inst;
    arm_status status;
    status = arm_rfft_fast_init_f32(&inst, fftLen);

    float32_t signalCombined[fftLen] = {0};
    float32_t fftCombined[fftLen] = {0};
    float32_t fftMagnitude[fftLen] = {0};
    volatile float32_t fftFreq[fftLen] = {0};

    float32_t maxAmp;
    uint32_t maxAmpInd;

    while (1)
    {
        for (int i = 0; i< fftLen; i++)
        {
            signalCombined[i] = 40 * arm_sin_f32(450 * i); // 450 frequency at 40 amplitude
        }

        arm_rfft_fast_f32(&inst, signalCombined, fftCombined, 0); // perhaps switch to complex transform to allow for negative frequencies?
        arm_cmplx_mag_f32(fftCombined, fftMagnitude, half_fftLen);
        fftMagnitude[0] = fftCombined[0];
        fftMagnitude[half_fftLen] = fftCombined[1];

        arm_max_f32(fftMagnitude, half_fftLen, &maxAmp, &maxAmpInd); // We need the 3 max values

        for (int k = 0; k < fftLen ; k++)
        {
            fftFreq[k] = ((k*sampleFreq)/fftLen);
        }
}

下面显示的是我从上面列出的代码中得到的结果:虽然我确实从算法中得到了一个幅度(在正确的索引 12 处),但它与频率或幅度不对应输入数组 signalCombined[].

有人知道为什么会这样吗?就像我的许多错误一样,这可能是一件非常微不足道和愚蠢的事情,但我终究无法弄清楚为什么会这样。


编辑:感谢 SleuthEye 的帮助,现在可以找到频率,因为生成 sin() 信号的初始方法不正确。

一些新问题出现了,因为 FFT 似乎只为 32 个样本产生了正确的频率,尽管 bin 大小相应地缩放以适应调整后的样本大小。

我也无法实现幅度固定算法:根据 SleuthEye's Link 示例代码 2*(1/N)*abs(X(k))^2 我自己实现了 2 * powf(fabs(fftMagnitude[j]), 2) / fftLen,如下面的代码所示,但这不会产生接近正确的结果。

while (1)
{

    for (int i = 0; i < fftLen; i++)
    {
        signalCombined[i] =     400 * arm_sin_f32(2 * PI * 450 * i / sampleFreq);       // Sin Alpha, 400 amp at 10 kHz
//                      700 * arm_sin_f32(2 * PI * 33000 * i / sampleFreq) +        // Sin Bravo, 700 amp at 33 kHz
//                          300 * arm_sin_f32(2 * PI * 50000 * i / sampleFreq);             // Sin Charlie, 300 amp at 50 kHz
    }

    arm_rfft_fast_f32(&inst, signalCombined, fftCombined, 0); // calculate the fourier transform of the time domain signal
    arm_cmplx_mag_f32(fftCombined, fftMagnitude, half_fftLen);  // calculate the magnitude of the fourier transform
    fftMagnitude[0] = fftCombined[0];
    fftMagnitude[half_fftLen] = fftCombined[1];

    for (int j = 0; j < sizeof(fftMagnitude); j++)
    {
        fftMagnitude[j] = 2 * powf(fabs(fftMagnitude[j]), 2) / fftLen; // Algorithm to fix the amplitude of each unique frequency
    }

    arm_max_f32(fftMagnitude, half_fftLen, &maxAmp, &maxAmpInd); // We need the 3 max values

    for (int k = 0; k < fftLen ; k++)
    {
        fftFreq[k] = ((k*sampleFreq)/fftLen);
    }
}

您的音调生成没有考虑 1600Hz 的采样频率,因此您实际上是在 450*1600/(2*PI) ~ 114591Hz 的频率下生成音调,混叠到 ~608Hz。当使用 32 的 FFT 大小时,该 608Hz 频率大致对应于 12 左右的频率索引。

在 1600Hz 采样频率下生成 450Hz 音调应按如下方式完成:

for (int i = 0; i< fftLen; i++)
{
    signalCombined[i] = 40 * arm_sin_f32(2 * PI * 450 * i / sampleFreq);
}

就匹配振幅而言,请注意 time-domain 和 frequency-domain 之间的比例因子约为 0.5*fftLen(参见 this other post of mine)。