Analiza plików WAV C (libsndfile, fftw3)

Question

Próbuję stworzyć prostą aplikację C, która może dać wartość od 0-100 w pewnym zakresie częstotliwości przy danym znaczniku czasu w pliku WAV.

Przykład: Mam zakres częstotliwości 44,1 kHz (typowy plik MP3) i chcę podzielić ten zakres na n zakresów (od 0). Następnie trzeba uzyskać amplitudę każdego zakresu, wynosi od 0 do 100.

Co udało mi się do tej pory:

Korzystanie libsndfile jestem teraz w stanie odczytać dane z wav plik.

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

float samples[sfinfo.frames]; 

sf_read_float(infile, samples, 1); 

Jednak moje zrozumienie FFT jest raczej ograniczone. Ale wiem, że to jest wymagane, aby uzyskać amplitudy w zakresie, którego potrzebuję. Ale jak przejść dalej? Znalazłem bibliotekę FFTW-3, która wydaje się być odpowiednia do tego celu.

znalazłem jakąś pomoc tutaj: https://stackoverflow.com/a/4371627/1141483

i spojrzał na tutorialu FFTW tutaj: http://www.fftw.org/fftw2_doc/fftw_2.html

Ale jestem pewien zachowania FFTW, nie wiem, aby przejść stąd .

I jeszcze jedno pytanie, zakładając, że używasz libsndfile: Jeśli wymusisz odczyt jako pojedynczy kanał (z plikiem stereo), a następnie odczytasz próbki. Czy w takim razie będziesz tylko czytał połowę próbek z całego pliku? Czy połowa z nich pochodzi z kanału 1, czy też automatycznie je odfiltrowuje?

Dzięki za pomoc.

EDIT: Mój kod można zobaczyć tutaj:

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

No to wszystko zależy od zakresu częstotliwości jesteś po. FFT działa, pobierając 2^n próbek i dostarczając 2^(n-1) liczb rzeczywistych i urojonych. Muszę przyznać, że jestem dość niepewny, co dokładnie reprezentują te wartości (mam przyjaciela, który obiecał przejść przez to wszystko ze mną zamiast pożyczki, którą mu zrobiłem, gdy miał problemy finansowe;)) innych niż kąt wokół okręgu. Skutecznie dostarczają one arccos parametru kątowego dla sinusa i cosinusa dla każdego przedziału częstotliwości, z którego oryginalne, 2^n próbki mogą być idealnie zrekonstruowane.

W każdym razie ma to ogromną zaletę, że można obliczyć wielkość, przyjmując odległość euklidesową części rzeczywistych i urojonych (sqrtf ((real * real) + (imag * imag))). Zapewni to niezniormalizowaną wartość odległości. Wartość ta może być następnie wykorzystana do zbudowania wielkości dla każdego pasma częstotliwości.

Pozwala więc przyjąć zamówienie 10 FFT (2^10). Wprowadzasz 1024 próbki. Próbkujecie te próbki, a otrzymacie 512 wartości urojonych i rzeczywistych (konkretna kolejność tych wartości zależy od używanego algorytmu FFT). Oznacza to, że w przypadku pliku audio o częstotliwości 44,1Khz każdy odbiornik reprezentuje 44100/512 Hz lub ~ 86 Hz na odbiornik.

Jedną rzeczą, która powinna się wyróżniać, jest to, że jeśli użyjesz więcej próbek (z tego, co nazywa się czasem lub przestrzenną domeną, gdy mamy do czynienia z wielowymiarowymi sygnałami, takimi jak obrazy), uzyskujesz lepszą reprezentację częstotliwości (w tym, co nazywa się domeną częstotliwości). Jednak poświęcasz jedno dla drugiego. Tak po prostu się dzieje i trzeba z tym żyć.

Zasadniczo konieczne będzie dostrojenie skrzynek częstotliwości i czasu/rozdzielczości przestrzennej w celu uzyskania wymaganych danych.

Najpierw trochę nomenklatury. Próbki z 1024 domenami czasowymi, o których wspomniałem wcześniej, nazywa się twoim oknem. Generalnie podczas wykonywania tego rodzaju procesu będziesz chciał przesuwać okno o pewną kwotę, aby uzyskać następne 1024 próbki, które FFT. Oczywistym rozwiązaniem byłoby pobranie próbek 0-> 1023, a następnie 1024-> 2047 i tak dalej. To niestety nie daje najlepszych rezultatów. Idealnie chciałbyś pokryć okna w pewnym stopniu, aby z czasem uzyskać płynniejszą zmianę częstotliwości. Najczęściej ludzie przesuwają okno o połowę rozmiaru okna. tzn. twoje pierwsze okno będzie 0-> 1023 drugie 512-> 1535 i tak dalej i tak dalej.

Teraz pojawia się kolejny problem. Podczas gdy te informacje zapewniają doskonałą odwrotną rekonstrukcję sygnału FFT, pozostawia to problem w pewnym stopniu, że częstotliwości przeciekają do pojemników otaczających. Aby rozwiązać ten problem, niektórzy matematycy (o wiele bardziej inteligentni ode mnie) wymyślili koncepcję window function. Funkcja okna zapewnia znacznie lepszą izolację częstotliwościową w domenie częstotliwości, ale prowadzi do utraty informacji w dziedzinie czasu (tj. Niemożliwe jest jej pełne odtworzenie po użyciu funkcji okna AFAIK).

Teraz istnieją różne rodzaje funkcji okna, począwszy od prostokątnego okna (skutecznie nie robiąc nic do sygnału) do różnych funkcji, które zapewniają znacznie lepszą izolację częstotliwości (chociaż niektóre mogą również zabijać otaczające częstotliwości, które mogą Cię zainteresować! !). Istnieje, niestety, nie jeden rozmiar pasuje do wszystkich, ale jestem wielkim fanem (dla spektrogramów) funkcji okna blackman-harris. Myślę, że daje to najlepiej wyglądające wyniki!

Jednak jak już wspomniałem wcześniej, FFT oferuje nieznehormalizowane widmo. Aby znormalizować widmo (po obliczeniu odległości euklidesowej), musisz podzielić wszystkie wartości przez współczynnik normalizacji (bardziej szczegółowo omówię here).

Ta normalizacja zapewni Ci wartość z zakresu od 0 do 1. Możesz więc łatwo zwiększyć tę wartość o 100, aby uzyskać skalę od 0 do 100.

To jednak nie jest miejsce, w którym się kończy. Widmo, które można uzyskać z tego, jest raczej niesatysfakcjonujące. To dlatego, że patrzysz na wielkość używając skali liniowej. Niestety ludzkie ucho słyszy używając skali logarytmicznej. To raczej powoduje problemy z wyglądem spektrogramu/widma.

Aby to obejść, należy przekonwertować te wartości z 0 na 1 (nazam to "x") na skalę decybeli. Standardowa transformacja to 20.0f * log10f(x). W ten sposób otrzymasz wartość, w której 1 został przekonwertowany na 0, a 0 na konwertowany na -infinity. wasze wielkości są teraz w odpowiedniej skali logarytmicznej. Jednak nie zawsze jest to pomocne.

W tym momencie należy przyjrzeć się oryginalnej głębokości bitowej próbki. Przy 16-bitowym próbkowaniu otrzymujesz wartość od 32767 do -32768. Oznacza to, że Twój dynamic range to fabsf (20.0f * log10f (1.0f/65536.0f)) lub ~ 96,33dB. Teraz mamy tę wartość.

Przyjmij wartości, które otrzymaliśmy z powyższego obliczenia dB. Dodaj do niego wartość -96.33. Oczywiście maksymalna amplituda (0) wynosi teraz 96,33. Teraz zrobiłem to samo, a teraz masz wartość od -infinity do 1.0f. Przymocuj dolny koniec do 0, a teraz masz zakres od 0 do 1 i pomnóż go przez 100, a otrzymasz końcowy zakres od 0 do 100.

I to jest znacznie bardziej potworny słupek niż pierwotnie zamierzałem, ale powinien dać dobre podstawy do generowania dobrego spektrum/spektrogramu dla sygnału wejściowego.

i oddychać

Dalsza lektura (dla osób innych niż oryginalny plakat, który już znalazł IT):

Converting an FFT to a spectogram

Edycja: Tak na marginesie znalazłem pocałunek FFT jest znacznie łatwiejszy w użyciu, mój kod do wykonania funkcji FFT jest następujący:

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; 
}