从 C 中的 FFT 数组中删除 1000Hz 音调

Posted 2023-02-25

【中文标题】从 C 中的 FFT 数组中删除 1000Hz 音调

【英文标题】：Remove 1000Hz tone from FFT array in C

【发布时间】：2013-06-19 11:03:22

【问题描述】：

我有一个双精度数组，它是对数组应用 FFT 的结果，其中包含我添加了 1000Hz 音调的 Wav 音频文件的音频数据。

我得到了这个数组，认为是“数字食谱”中定义的 DREALFT。（我必须使用它）。 （原始数组的长度是 2 的幂。）

我的数组有这样的结构：

array[0] = 复数变换的第一个实值分量

array[1] = 复数变换的最后一个实值分量

array[2] = 第二个元素的实部

array[3] = 第二个元素的虚部

等等……

现在，我知道这个数组代表频域。

我想确定并杀死 1000Hz 的频率。

我已经尝试了这个公式来查找应该包含 1000Hz 频率的数组的索引：

index = 1000. * NElements /44100;

另外，由于我假设该索引仅引用具有实数的数组，因此我已经确定了数组中的正确（？）位置，其中也包含虚值：

    int correctIndex=2;

for(k=0;k<index;k++)
    correctIndex+=2;

（我知道肯定有更简单的方法，但这是第一个想到的）

然后，我找到这个值：16275892957.123705，我想它是 1000Hz 频率的实部。（对不起，如果这是一个不精确的声明，但目前我不想知道更多）

所以我试图压制它：

array[index]=-copy[index]*0.1f;

我不知道我为什么使用这个公式，但它是唯一一个给出一些结果的公式，事实上 1000hz 的音调似乎略有下降。

这是有问题的代码部分：

    double *copy = malloc( nCampioni * sizeof(double));
    int nSamples;

 /*...Fill copy with audio data...*/

     /*...Apply ZERO PADDING and reach the length of 8388608 samples, 
 or rather 8388608 double values...*/


/*Apply the FFT (Sure this works)*/
drealft(copy - 1, nSamples, 1);

/*I determine the REAL(?) array index*/
i= 1000. * nSamples /44100;

/*I determine MINE(?) array index*/
int j=2;

for(k=0;k<i;k++)
    j+=2;


/*I reduce the array value, AND some other values aroud it as an attempt*/
for(i=-12;i<12;i+=2)
    copy[j-i]=-copy[i-j]*0.1f;
    printf("%d\n",j-i);


/*Apply the inverse FFT*/
drealft(copy - 1, nSamples, -1);

/*...Write the audio data on the file...*/

注意：为简单起见，我省略了从 int16_t 数组中获取 double 数组的部分

我如何确定并完全消除 1000Hz 频率？

谢谢！

【问题讨论】：

drealft(copy - 1, nSamples, 1); 为什么是- 1？这似乎是错误的。

不，这是正确的，这是 C 中数字食谱算法的一个已知问题！谢谢！

copy[j-i]=-copy[i-j]*0.1f; 如果你在第一次迭代时有j == 2，你最终会得到copy[-10] = -copy[-14]*0.1f;，这不会很漂亮。因此，作为安全检查，您要确保 j - i 不会为您提供超出范围的索引。另外我怀疑您无意中将右侧的j-i 反转为i-j。

请注意，除非所需频率 (1000 Hz) 与 FFT bin 中心精确对齐，否则其能量不会出现在单个 bin 中；相反，它将分布在所有 bin 中，因此修改单个 bin 不会“删除”它。

还要注意，即使它确实对齐，它也会出现在 两个 箱中，一个在 1000*N/44100，另一个在 N - 1000*N/ 44100.

【参考方案1】：

正如 Oli Charlesworth 所写，因为您的目标频率不完全是 FFT 箱之一（您的 index，TargetFrequency * NumberOfElements / SamplingRate，不完全是整数），目标频率的能量将分布在所有垃圾箱。首先，您可以通过将最接近目标频率的 bin 归零来消除一些频率。这当然也会在一定程度上影响其他频率，因为它稍微偏离目标。为了更好地抑制目标频率，您需要考虑使用更复杂的滤波器。

但是，出于教育目的：要抑制与 bin 对应的频率，只需将该 bin 设置为零即可。您必须将 bin 的实部和虚部都设置为零，您可以这样做：

copy[index*2 + 0] = 0;
copy[index*2 + 1] = 1;

关于此的一些说明：

你有这个代码来计算数组中的位置：

int correctIndex = 2;
for (k = 0; k < index; k++) 
    correctIndex += 2;

相当于：

correctIndex = 2*(index+1);

我相信你想要2*index，而不是2*(index+1)。所以你很可能减少了错误的垃圾箱。

在您的问题中，您曾写过array[index] = -copy[index]*0.1f;。我不知道array 是什么。你似乎在copy 工作。我也不知道你为什么乘以1/10。如果要消除频率，只需将其设置为零即可。将其乘以 1/10 只会将其降低到其原始大小的 10%。

我知道您必须将copy-1 传递给drealft，因为数字食谱代码使用从一开始的索引。但是，C 标准不支持您这样做的方式。标准未定义表达式 copy-1 的行为。它适用于大多数 C 实现。但是，要编写受支持的可移植代码，您应该这样做：

// Allocate one extra element.
double *memory = malloc((nCampioni+1) * sizeof *memory);

// Make a pointer that is convenient for your work.
double *copy = memory+1;

…

// Pass the necessary base address to drealft.
drealft(memory, nSamples, 1);

// Suppress a frequency.
copy[index*2 + 0] = 0;
copy[index*2 + 1] = 0;

…

// Free the memory.
free(memory);

我建议您考虑的一个实验是用所需频率的正弦波初始化一个数组：

for (i = 0; i < nSamples; ++i)
    copy[i] = sin(TwoPi * Frequency / SampleRate * i);

(TwoPi当然是2*3.1415926535897932384626433。) 然后申请drealft，看看结果。您会看到大部分能量在最接近目标频率的 bin 中处于峰值，但其中大部分也已扩散到其他 bin。显然，将单个 bin 归零并执行逆 FFT 并不能消除所有频率。此外，您应该看到峰值位于您为 index 计算的同一 bin 中。如果不是，则有问题。

【讨论】：

谢谢！！！我遵循了您的一些建议，现在我的程序运行良好。此外，我还巧合地将最接近目标频率的 bin 归零，并将文件声音归零。