WAV-файл анализа C (libsndfile, fftw3)

Question

Я пытаюсь разработать простое приложение на C, которое может дать значение от 0 до 100 в определенном диапазоне частот при заданной отметке времени в WAV-файле.

Пример: у меня частотный диапазон 44,1 кГц (типичный файл 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, которая, кажется, подходит для этой цели.

Я нашел некоторую помощь здесь: /questions/22432504/kak-poluchit-chastotyi-kazhdogo-znacheniya-v-bpf/22432519#22432519

и посмотрел учебник 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

Ну, все зависит от того, какой диапазон частот вам нужен. БПФ работает, беря 2^n выборок и предоставляя вам 2^(n-1) действительных и мнимых чисел. Я должен признать, что я довольно туманный в том, что именно представляют эти ценности (у меня есть друг, который пообещал пройти через все это со мной вместо кредита, который я сделал ему, когда у него были финансовые проблемы;)) кроме угол вокруг круга. По сути, они предоставляют вам arccos параметра угла для синуса и косинуса для каждого частотного бина, из которого можно идеально восстановить исходные 2^n выборки.

В любом случае это имеет огромное преимущество: вы можете рассчитать величину, взяв евклидово расстояние от реальной и мнимой частей (sqrtf( (real * real) + (imag * imag))). Это дает вам ненормированное значение расстояния. Это значение затем можно использовать для построения величины для каждой полосы частот.

Итак, давайте возьмем порядка 10 БПФ (2^10). Вы вводите 1024 образца. Вы БПФ эти образцы, и вы получаете 512 мнимых и реальных значений обратно (конкретный порядок этих значений зависит от используемого вами алгоритма БПФ). Таким образом, это означает, что для аудиофайла с частотой 44,1 кГц каждая ячейка представляет 44100/512 Гц или ~86 Гц на ячейку.

Из этого следует выделить одну вещь: если вы используете больше выборок (из того, что называется временной или пространственной областью при работе с многомерными сигналами, такими как изображения), вы получите лучшее представление частоты (в том, что называется частотной областью). Однако вы жертвуете одним ради другого. Так обстоят дела, и вам придется с этим жить.

По сути, вам нужно настроить частотные интервалы и временное / пространственное разрешение, чтобы получить необходимые данные.

Сначала немного номенклатуры. 1024 примера временной области, на которые я ссылался ранее, называются вашим окном. Как правило, при выполнении такого рода процессов вам нужно на некоторое время сдвинуть окно, чтобы получить следующие 1024 семпла, которые вы БПФ. Очевидное, что нужно сделать, это взять образцы 0->1023, затем 1024->2047 и так далее. Это, к сожалению, не дает лучших результатов. В идеале вы хотите до некоторой степени перекрывать окна, чтобы с течением времени получить более плавное изменение частоты. Чаще всего люди надевают окно на половину размера окна. т.е. ваше первое окно будет 0->1023, второе 512->1535 и так далее, и так далее.

Теперь это поднимает еще одну проблему. Хотя эта информация обеспечивает идеальное восстановление сигнала обратного БПФ, у вас возникает проблема, связанная с тем, что частоты до некоторой степени просачиваются в приемники объемного звучания. Чтобы решить эту проблему, некоторые математики (гораздо более умные, чем я) придумали концепцию оконной функции. Оконная функция обеспечивает гораздо лучшую частотную изоляцию в частотной области, хотя и приводит к потере информации во временной области (то есть невозможно полностью восстановить сигнал после использования оконной функции AFAIK).

Теперь существуют различные типы оконных функций, начиная от прямоугольного окна (фактически ничего не делая с сигналом) до различных функций, которые обеспечивают гораздо лучшую частотную изоляцию (хотя некоторые могут также убивать окружающие частоты, которые могут вас заинтересовать!!). Увы, не один размер подходит всем, но я большой поклонник (для спектрограмм) оконной функции Блэкмана-Харриса. Я думаю, что это дает лучшие результаты!

Однако, как я упоминал ранее, БПФ предоставляет вам ненормализованный спектр. Чтобы нормализовать спектр (после расчета евклидова расстояния), вам нужно разделить все значения на коэффициент нормализации (я здесь более подробно).

эта нормализация предоставит вам значение от 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.33 дБ. Так что теперь у нас есть это значение.

Возьмите значения, которые мы получили из расчета дБ выше. Добавьте к этому значение -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;
}