WAV 文件分析 C (libsndfile, fftw3)

Question

我正在尝试开发一个简单的 C 应用程序，它可以在 WAV 文件中的给定时间戳的特定频率范围内给出 0-100 的值。

示例：我的频率范围为 44.1kHz（典型的 MP3 文件），我想将该范围分成 n 个范围（从 0 开始）。然后我需要得到每个范围的幅度，从 0 到 100。

到目前为止我所管理的：

使用 libsndfile 我现在能够读取 WAV 文件的数据。

infile = sf_open(argv [1], SFM_READ, &sfinfo);

float samples[sfinfo.frames];

sf_read_float(infile, samples, 1);

但是，我对 FFT 的理解相当有限。但我知道为了在我需要的范围内获得振幅是必需的。但我该如何从这里继续前进？我找到了库 FFTW-3，它似乎很适合这个目的。

我在这里找到了一些帮助：https://stackoverflow.com/a/4371627/1141483

并在这里查看了 FFTW 教程：http://www.fftw.org/fftw2_doc/fftw_2.html

但由于我不确定 FFTW 的行为，我不知道从这里开始。

还有一个问题，假设您使用 libsndfile：如果您强制读取为单通道（使用立体声文件），然后读取样本。那么您实际上只会读取整个文件的一半样本吗？其中一半来自通道 1，还是会自动将其过滤掉？

非常感谢您的帮助。

编辑：我的代码可以在这里看到：

double blackman_harris(int n, int N){
double a0, a1, a2, a3, seg1, seg2, seg3, w_n;
a0 = 0.35875;
a1 = 0.48829;
a2 = 0.14128;
a3 = 0.01168;

seg1 = a1 * (double) cos( ((double) 2 * (double) M_PI * (double) n) / ((double) N - (double) 1) );
seg2 = a2 * (double) cos( ((double) 4 * (double) M_PI * (double) n) / ((double) N - (double) 1) );
seg3 = a3 * (double) cos( ((double) 6 * (double) M_PI * (double) n) / ((double) N - (double) 1) );

w_n = a0 - seg1 + seg2 - seg3;
return w_n;
}

int main (int argc, char * argv [])
{   char        *infilename ;
SNDFILE     *infile = NULL ;
FILE        *outfile = NULL ;
SF_INFO     sfinfo ;


infile = sf_open(argv [1], SFM_READ, &sfinfo);

int N = pow(2, 10);

fftw_complex results[N/2 +1];
double samples[N];

sf_read_double(infile, samples, 1);


double normalizer;
int k;
for(k = 0; k < N;k++){
    if(k == 0){

        normalizer = blackman_harris(k, N);

    } else {
        normalizer = blackman_harris(k, N);
    }

}

normalizer = normalizer * (double) N/2;



fftw_plan p = fftw_plan_dft_r2c_1d(N, samples, results, FFTW_ESTIMATE);

fftw_execute(p);


int i;
for(i = 0; i < N/2 +1; i++){
    double value = ((double) sqrtf(creal(results[i])*creal(results[i])+cimag(results[i])*cimag(results[i]))/normalizer);
    printf("%f\n", value);

}



sf_close (infile) ;

return 0 ;
} /* main */

Answer 1

这一切都取决于您所追求的频率范围。 FFT 通过采集 2^n 个样本并为您提供 2^(n-1) 个实数和虚数来工作。我不得不承认，我对这些价值观到底代表什么感到很模糊（我有一个朋友承诺会和我一起经历这一切，而不是在他遇到财务问题时我借给他；））围绕一个圆的一个角。实际上，它们为您提供了每个频率区间的正弦和余弦的角度参数的 arccos，可以从中完美地重建原始 2^n 样本。

无论如何，这具有巨大的优势，您可以通过取实部和虚部的欧几里德距离 (sqrtf( (real * real) + (imag * imag) )) 来计算幅度。这为您提供了一个非标准化的距离值。然后可以使用该值来为每个频带构建幅度。

所以让我们订购 10 FFT (2^10)。您输入了 1024 个样本。您对这些样本进行 FFT，然后返回 512 个虚值和实值（这些值的特定顺序取决于您使用的 FFT 算法）。因此，这意味着对于 44.1Khz 音频文件，每个 bin 代表 44100/512 Hz 或每个 bin 约 86Hz。

应该从中脱颖而出的一件事是，如果您使用更多样本（在处理图像等多维信号时来自所谓的时域或空间域），您将获得更好的频率表示（在所谓的频域中） .然而，你为另一个牺牲了一个。这就是事情的发展方式，你将不得不忍受它。

基本上，您需要调整频率箱和时间/空间分辨率以获得所需的数据。

首先是一些命名法。我前面提到的 1024 个时域样本称为你的窗口。通常，在执行此类过程时，您需要将窗口滑动一些量以获得您 FFT 的下一个 1024 个样本。显而易见的做法是抽取样本 0->1023，然后是 1024->2047，依此类推。不幸的是，这并没有给出最好的结果。理想情况下，您希望在一定程度上重叠窗口，以便随着时间的推移获得更平滑的频率变化。最常见的是人们将窗口滑动一半的窗口大小。即你的第一个窗口将是 0->1023 第二个 512->1535 等等。

现在这又带来了另一个问题。虽然此信息提供了完美的逆 FFT 信号重建，但它给您留下了一个问题，即频率在一定程度上泄漏到环绕声箱中。为了解决这个问题，一些数学家（比我聪明得多）提出了window function 的概念。窗函数在频域中提供了更好的频率隔离，但会导致时域中的信息丢失（即，在您使用窗函数之后，不可能完美地重建信号，AFAIK）。

现在有各种类型的窗口函数，从矩形窗口（实际上对信号没有任何作用）到提供更好频率隔离的各种函数（尽管有些函数也可能会消除您可能感兴趣的周围频率！！ ）。唉，没有一种尺寸适合所有人，但我是 blackmann-harris 窗口函数的忠实粉丝（对于频谱图）。我认为它给出了最好看的结果！

但是，正如我之前提到的，FFT 为您提供了一个非归一化的频谱。要对光谱进行归一化（在计算欧几里得距离之后），您需要将所有值除以归一化因子（我会更详细地介绍 here）。

这种标准化将为您提供一个介于 0 和 1 之间的值。因此您可以轻松地将该值乘以 100 以获得 0 到 100 的比例。

然而，这并不是它的结束。您从中获得的频谱相当不令人满意。这是因为您正在使用线性比例查看幅度。不幸的是，人耳使用对数刻度听到。这反而会导致频谱图/频谱的外观出现问题。

要解决这个问题，您需要将这些 0 到 1 的值（我称之为“x”）转换为分贝刻度。标准转换为20.0f * log10f( x )。然后，这将为您提供一个值，其中 1 已转换为 0，0 已转换为 -infinity。您的大小现在处于适当的对数刻度中。然而，它并不总是那么有用。

此时您需要查看原始样本位深度。在 16 位采样时，您会得到一个介于 32767 和 -32768 之间的值。这意味着您的 dynamic range 是 fabsf( 20.0f * log10f( 1.0f / 65536.0f ) ) 或 ~96.33dB。所以现在我们有了这个值。

取我们从上面的 dB 计算中得到的值。将这个 -96.33 值添加到它。显然，最大幅度 (0) 现在是 96.33。现在按相同的值进行除法运算，您现在的值范围从 -infinity 到 1.0f。将下端限制为 0，您现在有一个从 0 到 1 的范围，然后将其乘以 100，您就有了最终的 0 到 100 范围。

这比我最初的意图更像一个怪物帖子，但应该为您提供如何为输入信号生成良好频谱/频谱图的良好基础。

呼吸

进一步阅读（对于已经找到它的原始海报以外的人）：

Converting an FFT to a spectogram

编辑：顺便说一句，我发现 Kiss FFT 更容易使用，我执行前向 fft 的代码如下：

CFFT::CFFT( unsigned int fftOrder ) :
    BaseFFT( fftOrder )
{
    mFFTSetupFwd    = kiss_fftr_alloc( 1 << fftOrder, 0, NULL, NULL );
}

bool CFFT::ForwardFFT( std::complex< float >* pOut, const float* pIn, unsigned int num )
{
    kiss_fftr( mFFTSetupFwd, pIn, (kiss_fft_cpx*)pOut );
    return true;
}