我正在尝试开发一个简单的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教程: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 */
发布于 2012-05-17 14:00:13
这完全取决于你所追求的频率范围。FFT的工作原理是获取2^n样本,并为你提供2^(n-1)实数和虚数。我不得不承认,我对这些价值观到底代表了什么(我有一个朋友承诺要和我一起度过这一切,而不是在他有财务问题时我给他的一笔贷款),我很模糊),而不是一个圆圈的角度。有效地,他们提供了一个arccos的角度参数的正弦和余弦的每一个频率宾,其中原始的2^n样本,可以完美地,重建。
无论如何,这有一个巨大的优势,你可以通过取实部和虚部(sqrtf( (实*实)+ (imag * imag ) )的欧几里得距离来计算震级。这为您提供了一个未规范化的距离值。然后,该值可用于为每个频带构建一个幅度。
因此,让我们采用10 FFT (2^10)的顺序。你输入了1024个样本。您可以对这些样本进行FFT处理,得到512个假想值和实数(这些值的特殊顺序取决于您使用的FFT算法)。这意味着,对于一个44.1Khz音频文件,每个bin代表44100/512 Hz或~86 Hz/ bin。
有一点值得注意的是,如果您使用更多的样本(在处理多维信号(如图像)时,使用更多的样本(从所谓的时间域或空间域),则可以获得更好的频率表示(在所谓的频域中)。不管你为了另一个牺牲一个。事情就是这样进行的,你必须接受它。
基本上,您将需要调整频率箱和时间/空间分辨率,以获得您需要的数据。
首先来点命名。我前面提到的1024个时域示例称为“您的窗口”。通常,在执行这类处理时,您需要在窗口上滑动一定数量的窗口,以获得下一个1024个样本(您的FFT )。最明显的做法是取样本0->1023,然后1024->2047,以此类推。不幸的是,这并没有给出最好的结果。理想情况下,您希望在一定程度上重叠窗口,以便随着时间的推移,您的频率变化会更平稳。最常见的情况是,人们把窗户滑动到半个窗口大小。你的第一个窗口是0->1023,第二个512->1535等等。
现在,这就引出了另一个问题。虽然这一信息提供了完美的逆FFT信号重建,它给你留下了一个问题,频率泄漏到环抱箱在一定程度上。为了解决这个问题,一些数学家(比我聪明得多)提出了窗口函数的概念。窗口函数在频域中提供了更好的频率隔离,尽管它会导致时域信息的丢失(即在使用窗口函数AFAIK之后不可能完美地重新构造信号)。
现在有各种类型的窗口函数,从矩形窗口(实际上对信号没有任何作用)到提供更好频率隔离的各种功能(尽管有些函数还可能扼杀您可能感兴趣的周围频率!!)唉,没有一个尺寸适合所有的人,但我是一个狂热的球迷(谱图)的黑曼-哈里斯窗口功能。我认为它能给出最好看的结果!
然而,正如我前面提到的,FFT为你提供了一个未正常化的光谱。要使光谱正常化(在欧几里德距离计算之后),您需要将所有值除以一个归一化因子(我将详细介绍这里)。
这种标准化将为您提供一个介于0到1之间的值。因此,您可以轻松地将此值乘以100,以获得0到100的比例。
然而,这并不是它的终点。你从这得到的光谱是相当不令人满意的。这是因为你用线性尺度来观察震级。不幸的是,人类耳朵听到的是对数音阶。这就导致了光谱图/光谱的外观问题。
为了解决这个问题,您需要将这0到1的值(我称之为'x')转换为分贝标度。标准转换是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范围。
这是一个怪物的帖子,比我原来的打算,但应该给你一个良好的基础如何产生一个良好的光谱/光谱图的输入信号。
呼吸
进一步阅读(供已发现原海报的人参考):
编辑:顺便提一下,我发现kiss使用起来容易得多,我执行前向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;
}
https://stackoverflow.com/questions/10627517
复制