从FFT阵列中移除1000Hz的音调

Question

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

我通过"数值食谱"中定义的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("%dn",j-i);
}
/*Apply the inverse FFT*/
drealft(copy - 1, nSamples, -1);
/*...Write the audio data on the file...*/

注意:为了简单起见，我省略了从int16_t

数组中获得double类型数组的部分。

我怎样才能确定并完全消除1000Hz的频率?

谢谢!

Answer 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，因为Numerical Recipes代码使用基于1的索引。然而，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归零并执行逆FFT不能消除所有频率。此外，您应该看到峰值位于您为index计算的同一个bin中。