WAV-Dateianalyse C (libsndfile, fftw3)

Question

Ich versuche eine einfache C-Anwendung zu entwickeln, die bei einem bestimmten Zeitstempel in einer WAV-Datei einen Wert von 0-100 in einem bestimmten Frequenzbereich geben kann.

Beispiel: Ich habe einen Frequenzbereich von 44,1 kHz (typische MP3-Datei) und möchte diesen Bereich in n Bereiche aufteilen (beginnend mit 0). Ich muß dann die Amplitude jeden Bereich bekommen, von 0 bis 100 seine

Was ich bisher geschaffen habe:

Mit libsndfile ich nun in der Lage bin, die Daten eines WAV- zu lesen Datei.

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

float samples[sfinfo.frames]; 

sf_read_float(infile, samples, 1); 

Allerdings ist mein Verständnis von FFT eher begrenzt. Aber ich weiß, dass es erforderlich ist, um die Amplituden in den Bereichen zu bekommen, die ich brauche. Aber wie gehe ich von hier weiter? Ich fand die Bibliothek FFTW-3, die für diesen Zweck geeignet scheint.

fand ich etwas Hilfe hier: https://stackoverflow.com/a/4371627/1141483

und sah hier im FFTW Tutorial: http://www.fftw.org/fftw2_doc/fftw_2.html

Aber wie ich über das Verhalten des FFTW nicht sicher bin, weiß ich nicht von hier, um die Fortschritte .

Und eine andere Frage, vorausgesetzt Sie verwenden libsndfile: Wenn Sie erzwingen, dass der Messwert single-channeled (mit einer Stereo-Datei) und lesen Sie dann die Proben. Werden Sie dann tatsächlich nur die Hälfte der Samples der gesamten Datei lesen? Wenn die Hälfte von ihnen aus Kanal 1 kommt oder filtert sie diese automatisch aus?

Vielen Dank für Ihre Hilfe.

EDIT: Mein Code kann hier gesehen werden:

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

Nun, es hängt alles von dem Frequenzbereich, nachdem Sie sind. Eine FFT arbeitet mit 2^n Samples und liefert Ihnen 2^(n-1) reelle und imaginäre Zahlen. Ich muss zugeben, dass ich ziemlich unklar bin, was genau diese Werte darstellen (ich habe einen Freund, der versprochen hat, alles mit mir durchzugehen anstelle eines Kredits, den ich ihm gemacht habe, als er finanzielle Probleme hatte;)) ein Winkel um einen Kreis. Effektiv liefern sie Ihnen einen Arccos des Winkelparameters für einen Sinus und einen Kosinus für jedes Frequenzfach, aus dem die ursprünglichen 2^n Abtastwerte perfekt rekonstruiert werden können.

Wie auch immer, dies hat den großen Vorteil, dass Sie die Größe berechnen können, indem Sie den euklidischen Abstand der Real- und Imaginärteile nehmen (sqrtf ((real * real) + (imag * imag))). Dadurch erhalten Sie einen unnormierten Abstandswert. Dieser Wert kann dann verwendet werden, um eine Größe für jedes Frequenzband zu bilden.

Also nehmen wir eine Bestellung 10 FFT (2^10). Sie geben 1024 Samples ein. Sie fasten diese Samples und Sie erhalten 512 imaginäre und reelle Werte zurück (die bestimmte Reihenfolge dieser Werte hängt von dem verwendeten FFT-Algorithmus ab). Das bedeutet, dass für eine 44.1Khz Audiodatei jeder Bin 44100/512 Hz oder ~ 86Hz pro Bin repräsentiert. Eine Sache, die sich davon abheben sollte ist, dass, wenn Sie mehr Samples verwenden (von dem, was Zeit oder Raumdomäne genannt wird, wenn Sie mit mehrdimensionalen Signalen wie Bildern arbeiten), eine bessere Frequenzrepräsentation erhalten). Wie auch immer du opferst für den anderen. Das ist genau die Art und Weise, wie die Dinge laufen und du musst damit leben.

Grundsätzlich müssen Sie die Frequenzbins und die zeitliche/räumliche Auflösung abstimmen, um die gewünschten Daten zu erhalten.

Zuerst ein bisschen Nomenklatur. Die 1024 Zeitdomänen-Beispiele, auf die ich früher Bezug nahm, werden als Ihr Fenster bezeichnet. Im Allgemeinen, wenn Sie diese Art von Prozess durchführen, möchten Sie das Fenster um einen gewissen Betrag schieben, um die nächsten 1024 Samples zu erhalten, die Sie FFT haben. Die naheliegende Sache wäre, die Beispiele 0-> 1023, dann 1024-> 2047 usw. zu nehmen. Dies liefert leider nicht die besten Ergebnisse. Idealerweise sollten Sie die Fenster zu einem gewissen Grad überlappen, damit Sie im Laufe der Zeit eine gleichmäßigere Frequenzänderung erhalten. Am häufigsten wird das Fenster um eine halbe Fenstergröße verschoben. dh Ihr erstes Fenster wird 0-> 1023 sein, das zweite 512-> 1535 und so weiter und so fort.

Das bringt nun ein weiteres Problem auf. Während diese Information für eine perfekte inverse FFT-Signalrekonstruktion sorgt, ergibt sich das Problem, dass Frequenzen gewissermaßen in Surround-Bins münden. Um dieses Problem zu lösen, kamen einige Mathematiker (viel intelligenter als ich) auf das Konzept einer window function. Die Fensterfunktion sorgt für eine weit bessere Frequenzisolierung im Frequenzbereich, führt jedoch zu einem Informationsverlust im Zeitbereich (dh es ist unmöglich, das Signal nach Verwendung einer Fensterfunktion, AFAIK, perfekt wieder aufzubauen).

Jetzt gibt es verschiedene Arten von Fensterfunktionen, angefangen vom rechteckigen Fenster (effektiv nichts mit dem Signal zu tun) bis zu verschiedenen Funktionen, die eine weit bessere Frequenzisolierung bieten (obwohl einige auch die umliegenden Frequenzen zerstören könnten, die Sie interessieren könnten! !). Es gibt leider keine Einheitsgröße für alle, aber ich bin ein großer Fan (für Spektrogramme) der Blackmann-Harris-Fensterfunktion. Ich denke, es gibt die besten Ergebnisse!

Wie bereits erwähnt, bietet die FFT jedoch ein unnormalisiertes Spektrum. Um das Spektrum zu normalisieren (nach der euklidischen Abstandsberechnung) müssen Sie alle Werte durch einen Normierungsfaktor dividieren (ich gehe genauer auf here).

Diese Normalisierung liefert Ihnen einen Wert zwischen 0 und 1. Sie können also einfach diesen Wert um 100 multiplizieren, um Ihren Maßstab von 0 bis 100 zu erhalten.

Dies ist jedoch nicht, wo es endet. Das Spektrum, das Sie daraus erhalten, ist eher unbefriedigend. Dies liegt daran, dass Sie die Größe mit einer linearen Skala betrachten. Leider hört das menschliche Ohr eine logarithmische Skala. Dies verursacht eher Probleme mit dem Aussehen eines Spektrogramms/Spektrums.

Um das zu umgehen, müssen Sie diese 0 zu 1 Werte (ich nenne es "x") auf die Dezibelskala umrechnen. Die Standardumwandlung ist 20.0f * log10f(x). Dies liefert Ihnen dann einen Wert, wobei 1 zu 0 konvertiert wurde und 0 zu -infinity konvertiert wurde. Eure Größen sind jetzt in der entsprechenden logarithmischen Skala. Aber es ist nicht immer so hilfreich.

An dieser Stelle müssen Sie in die ursprüngliche Bittiefe des Samples schauen. Bei 16-Bit-Sampling erhalten Sie einen Wert zwischen 32767 und -32768. Das heißt, Ihr dynamic range ist fabsf (20.0f * log10f (1.0f/65536.0f)) oder ~ 96.33dB. Jetzt haben wir diesen Wert.

Nehmen Sie die Werte, die wir aus der obigen dB-Berechnung erhalten haben. Fügen Sie diesen Wert -96,33 hinzu. Offensichtlich ist die maximale Amplitude (0) jetzt 96,33. Jetzt didivde durch den gleichen Wert und Sie haben jetzt einen Wert von -unendlich bis 1,0f. Klemmen Sie das untere Ende auf 0 und Sie haben jetzt einen Bereich von 0 bis 1 und multiplizieren Sie das mit 100 und Sie haben Ihren letzten 0 bis 100 Bereich.

Und das ist viel mehr ein Monster Beitrag als ich ursprünglich beabsichtigt hatte, aber sollte Ihnen eine gute Grundlage in wie ein gutes Spektrum/Spektrogramm für ein Eingangssignal zu erzeugen.

und atmen

Weiterführende Literatur (für Menschen, andere als die Original-Poster, die es bereits gefunden hat):

Converting an FFT to a spectogram

bearbeiten: Als Nebenwirkung fand ich Kuss FFT viel einfacher zu verwenden, ist mein Code, um ein Vorwärtsfft durchzuführen, wie folgt:

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