Analisi file WAV C (libsndfile, fftw3)

Question

Sto cercando di sviluppare una semplice applicazione C in grado di fornire un valore compreso tra 0 e 100 in un determinato intervallo di frequenza a un determinato timestamp in un file WAV.

Esempio: ho un intervallo di frequenza di 44,1 kHz (tipico file MP3) e voglio suddividere tale intervallo in n quantità di intervalli (a partire da 0). Ho quindi bisogno di ottenere l'ampiezza di ciascuna gamma, essendo da 0 a 100.

Quello che sono riuscito finora:

Utilizzando libsndfile sono ora in grado di leggere i dati di un WAV- file.

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

float samples[sfinfo.frames]; 

sf_read_float(infile, samples, 1); 

Tuttavia, la mia comprensione di FFT è piuttosto limitata. Ma so che è necessario inorder per ottenere le ampiezze alle gamme che mi servono. Ma come posso andare avanti da qui? Ho trovato la libreria FFTW-3, che sembra essere adatta allo scopo.

ho trovato un po 'di aiuto qui: https://stackoverflow.com/a/4371627/1141483

e guardò il tutorial FFTW qui: http://www.fftw.org/fftw2_doc/fftw_2.html

Ma, come io sono sicuro circa il comportamento del FFTW, non so a progredire da qui .

E un'altra domanda, supponendo che si usi libsndfile: Se si forza la lettura per essere single channeled (con un file stereo) e quindi leggere i campioni. In questo momento leggerete solo metà dei campioni del file totale? Poiché la metà di essi proviene dal canale 1 o li filtra automaticamente?

Grazie mille per il vostro aiuto.

EDIT: Il mio codice può essere visto qui:

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

Beh, tutto dipende dalla gamma di frequenze che stai cercando. Un FFT funziona prendendo 2^n campioni e fornendoti numeri reali e immaginari 2^(n-1). Devo ammettere che sono piuttosto confuso su cosa rappresentino esattamente questi valori (ho un amico che mi ha promesso di passare tutto con me al posto di un prestito che gli ho fatto quando aveva problemi finanziari;)) diversi da un angolo attorno a un cerchio. Effettivamente forniscono un arccos del parametro angolare per un seno e coseno per ciascun intervallo di frequenze dal quale i campioni originali 2^n possono essere, perfettamente, ricostruiti.

In ogni caso questo ha l'enorme vantaggio di poter calcolare la magnitudine prendendo la distanza euclidea delle parti reali e immaginarie (sqrtf ((real * real) + (imag * imag))). Questo ti fornisce un valore di distanza non normalizzato. Questo valore può quindi essere usato per costruire una grandezza per ogni banda di frequenza.

Quindi, facciamo un ordine 10 FFT (2^10). Immetti 1024 campioni. FFT quei campioni e ottieni 512 valori reali e immaginari (il particolare ordinamento di questi valori dipende dall'algoritmo FFT che usi). Questo significa che per un file audio da 44.1 Khz ogni scomparto rappresenta 44100/512 Hz o ~ 86Hz per bin.

Una cosa che dovrebbe risaltare è che se si utilizzano più campioni (da quello che viene chiamato il dominio temporale o spaziale quando si gestiscono segnali multidimensionali come le immagini) si ottiene una migliore rappresentazione in frequenza (in che cosa è chiamato dominio della frequenza). Comunque tu sacrifichi uno per l'altro. Questo è solo il modo in cui vanno le cose e dovrai conviverci.

Fondamentalmente è necessario regolare i bin di frequenza e la risoluzione tempo/spazio per ottenere i dati richiesti.

Prima un po 'di nomenclatura. I campioni del dominio del tempo 1024 cui ho fatto riferimento in precedenza sono chiamati finestra. Generalmente quando si esegue questo tipo di processo, si vorrà far scorrere la finestra di qualche ammontare per ottenere i successivi 1024 campioni di FFT. La cosa ovvia da fare sarebbe prendere i campioni 0-> 1023, quindi 1024-> 2047 e così via. Questo purtroppo non dà i risultati migliori. Idealmente, si desidera sovrapporre le finestre in una certa misura in modo da ottenere una variazione di frequenza più uniforme nel tempo. Più comunemente le persone fanno scorrere la finestra di mezza dimensione della finestra. cioè la tua prima finestra sarà 0-> 1023 il secondo 512-> 1535 e così via e così via.

Ora questo fa apparire un ulteriore problema. Mentre questa informazione fornisce una perfetta ricostruzione inversa del segnale FFT, ti lascia un problema che le frequenze perdono in una certa misura nei cassonetti surround. Per risolvere questo problema alcuni matematici (molto più intelligenti di me) hanno ideato il concetto di window function. La funzione finestra fornisce un isolamento della frequenza molto migliore nel dominio della frequenza, ma porta a una perdita di informazioni nel dominio del tempo (cioè è impossibile ricostruire perfettamente il segnale dopo aver usato una funzione finestra, AFAIK).

Ora ci sono vari tipi di funzioni della finestra che vanno dalla finestra rettangolare (che effettivamente non fa nulla al segnale) a varie funzioni che forniscono un isolamento della frequenza molto migliore (anche se alcuni potrebbero anche uccidere le frequenze circostanti che potrebbero interessarti! !). Non c'è, ahimè, nessuna taglia adatta a tutti, ma io sono un grande fan (per gli spettrogrammi) della funzione di finestra blackmann-harris. Penso che dia i risultati migliori!

Tuttavia, come detto in precedenza, la FFT offre uno spettro non normalizzato. Per normalizzare lo spettro (dopo il calcolo della distanza euclidea) è necessario dividere tutti i valori per un fattore di normalizzazione (passo più in dettaglio here).

questa normalizzazione fornirà un valore compreso tra 0 e 1. Quindi è possibile moltiplicare facilmente questo valore per 100 per ottenere la scala da 0 a 100.

Questo, tuttavia, non è dove finisce. Lo spettro che ottieni da questo è piuttosto insoddisfacente. Questo perché stai osservando la magnitudine usando una scala lineare. Sfortunatamente l'orecchio umano si sente usando una scala logaritmica. Questo piuttosto causa problemi con l'aspetto di uno spettrogramma/spettro.

Per aggirare questo è necessario convertire questi valori da 0 a 1 (lo chiamerò 'x') alla scala del decibel. La trasformazione standard è 20.0f * log10f(x). Ciò fornirà quindi un valore per cui 1 è stato convertito in 0 e 0 è stato convertito in -infinity. le tue grandezze sono ora nella scala logaritmica appropriata. Tuttavia non è sempre così utile.

A questo punto è necessario esaminare la profondità di bit del campione originale. A campionamento a 16 bit si ottiene un valore compreso tra 32767 e -32768. Ciò significa che il tuo dynamic range è fabsf (20.0f * log10f (1.0f/65536.0f)) o ~ 96.33dB. Quindi ora abbiamo questo valore.

Prendi i valori che abbiamo ottenuto dal calcolo in dB sopra. Aggiungi questo valore -96.33 ad esso. Ovviamente l'ampiezza massima (0) è ora 96.33. Ora ho fatto con lo stesso valore e ora hai un valore che va da -infinity a 1.0f. Blocca l'estremità inferiore su 0 e ora hai un intervallo da 0 a 1 e moltiplica quello per 100 e hai il tuo intervallo da 0 a 100 finale.

E questo è molto più di un post mostro di quanto avessi inizialmente previsto, ma dovrebbe darvi una buona base su come generare un buon spettro/spettrogramma per un segnale in ingresso.

e respirare

Ulteriori letture (per le persone diverse dal manifesto originale che ha già trovato):

Converting an FFT to a spectogram

Edit: Per inciso ho trovato bacio FFT molto più facile da usare, il mio codice per eseguire un forward fft è il seguente:

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