Tho*_*num 10 c fft wav libsndfile fftw
我正在尝试开发一个简单的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 */
Goz*_*Goz 15
那一切都取决于你所追求的频率范围.FFT通过采用2 ^ n个样本并为您提供2 ^(n-1)个实数和虚数来工作.我不得不承认我对这些价值所代表的东西非常朦胧(我有一位朋友已答应与我一起完成所有这些,而不是在他遇到财务问题时给我的贷款;)除了围绕一个圆的角度.实际上,它们为每个频率区提供正弦和余弦角度参数的arccos,原始的2 ^ n样本可以完美地重建.
无论如何,这有一个巨大的优势,你可以通过获取实部和虚部的欧氏距离来计算幅度(sqrtf((真实*真实)+(图像*图像))).这为您提供了非标准化的距离值.然后,该值可用于为每个频带建立幅度.
所以我们下订单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信号重建,但它会让您遇到频率在某种程度上泄漏到环绕声箱中的问题.为了解决这个问题,一些数学家(比我聪明得多)想出了一个窗口函数的概念.窗口函数在频域中提供了更好的频率隔离,但导致时域中的信息丢失(即,在使用窗口函数AFAIK之后,它不可能完美地重构信号).
现在有各种类型的窗口功能,从矩形窗口(有效地对信号无效)到提供更好的频率隔离的各种功能(尽管有些也可能会杀死你可能感兴趣的周围频率!!).唉,没有一种尺寸适合所有人,但我是blackmann-harris窗函数的忠实粉丝(用于光谱图).我认为它给出了最好看的结果!
但是正如我之前提到的,FFT为您提供了非标准化频谱.要对频谱进行归一化(在欧氏距离计算之后),您需要将所有值除以归一化因子(我在此处详细介绍).
此规范化将为您提供0到1之间的值.因此您可以轻松地将此值乘以100以获得0到100的比例.
然而,这并不是它的结束.你从中获得的光谱相当令人不满意.这是因为您使用线性刻度查看幅度.不幸的是,人耳听到使用对数刻度.这相当导致频谱图/频谱看起来如何.
为了解决这个问题,你需要将这些0到1的值(我称之为'x')转换为分贝比例.标准转换为20.0f*log10f(x).然后,这将为您提供一个值,其中1已转换为0,0已转换为-infinity.你的数量现在处于适当的对数范围.然而,它并不总是那么有用.
此时,您需要查看原始样本位深度.在16位采样时,您将获得介于32767和-32768之间的值.这意味着您的动态范围是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范围.
这比我原本打算的更像是一个怪物帖子,但应该为你如何为输入信号生成一个好的频谱/频谱图提供良好的基础.
呼吸
进一步阅读(对于已经找到它的原始海报以外的人):
编辑:作为一个旁边我发现吻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;
}