c - 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
但由于我不确定飞行队的行为，我不知道从这里开始。
还有一个问题，假设您使用libsndfile：如果您强制读取为单通道（使用立体声文件），然后读取示例。你真的只会读一半的样本吗？因为其中一半来自第一频道，还是会自动过滤掉？
非常感谢你的帮助。
编辑：我的代码可以在这里看到：

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（（实部*实部）+（imag*imag）））来计算震级。这将为您提供一个不规范的距离值。然后，可以使用该值为每个频带建立一个幅值。
所以我们取10阶fft（2^10）。输入1024个样本。对这些样本进行fft运算，得到512个虚值和实值（这些值的特定顺序取决于使用的fft算法）。因此，这意味着对于44.1kHz音频文件，每个bin代表44100/512Hz或~86Hz/bin。
其中一个突出的特点是，如果你使用更多的样本（在处理图像等多维信号时来自所谓的时间或空间域），你会得到更好的频率表示（在所谓的频率域）。不管你如何牺牲一个来换取另一个。这就是事情的发展方向，你必须接受它。
基本上，您需要调整频率箱和时间/空间分辨率，以获得所需的数据。
首先是一些术语。我前面提到的1024个时域示例称为您的窗口。通常，在执行这种类型的处理时，您需要将窗口滑动一定量，以获取下一个1024个fft采样。很明显，要做的事情是采集样本0->1023，然后1024->2047，以此类推。不幸的是，这并不能带来最好的结果。理想情况下，您希望在一定程度上重叠窗口，以便随着时间的推移获得更平滑的频率变化。最常见的是人们把窗户滑动一半。你的第一个窗口是0->1023，第二个是512->1535，以此类推。
现在这又引出了一个问题。虽然这些信息提供了完美的反fft信号重建，但它在一定程度上给您留下了一个问题，即频率会泄漏到包围盒中。为了解决这个问题，一些数学家（比我聪明得多）提出了awindow function的概念。窗口函数在频域提供了更好的频率隔离，但会导致时域信息的丢失（即在使用窗口函数afaik后无法完全重新构造信号）。
现在有各种类型的窗口函数，从矩形窗口（实际上对信号没有任何作用）到提供更好的频率隔离的各种函数（尽管有些函数可能会杀死您可能感兴趣的周围频率！）！）。唉，没有一种尺寸适合所有人，但我是Blackmann-Harris窗口函数的忠实粉丝。我认为它给了最好的结果！
不过，正如我前面提到的，fft为您提供了一个非标准化的频谱。要使光谱正常化（在欧几里德距离计算之后），需要用一个标准化因子来划分所有值（我将更详细地介绍here）。
这种正常化将为您提供一个介于0和1之间的值。所以你可以很容易地将这个值乘以100，得到你的0到100的比例。
然而，这并不是它的终点。你从中得到的光谱相当不令人满意。这是因为你是在用线性比例来观察震级。不幸的是，人类的耳朵使用对数音阶。这会引起光谱图/光谱的外观问题。
要解决这个问题，您需要将这些0到1的值（我称之为“x”）转换为分贝刻度。标准转换是20.0f * log10f( x )。这将为您提供一个值，其中1已转换为0，0已转换为-无穷大。你的震级现在在适当的对数范围内。然而，这并不总是那么有帮助。
此时，您需要查看原始采样位深度。在16位采样时，会得到一个介于32767和-32768之间的值。这意味着您的dynamic range是fabsf（20.0f*log10f（1.0f/65536.0f））或~96.33db。所以现在我们有了这个价值。
取上面db计算得到的值。加上-96.33这个值。显然，最大振幅（0）现在是96.33。现在didivde得到了相同的值，你现在得到了一个从-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;
}