3. Bestimmung der Fundamentalfrequenz eines Audiosignals

Bestimmung der Fundamentalfrequenz eines Audiosignals

  • P

    Liebe C++-Community,

    Ich hab ein kleines C++-Programm geschrieben, dass ein Audiosignal aufnimmt und davon die Fundamentalfrequenz berechnen soll.

    Das ganze läuft so ab:
    Die Aufnahme erfolgt mit 8 KHz (mono, 16 bit) und alle 10 ms wird mit den bis dahin aufgenommenen 800 Samples eine FFT-Transformation ausgeführt.
    Das Signal so wie das jeweilige Ergebnis der FFT werden in eine .csv-Datei geschrieben, auf diese Weise kann ich die Werte in Mathematica plotten.

    Leider schaut das Fourrier-Signal nicht so aus, wie ich es erwartet hätte (ich hab dann die FFT kurzfristig durch eine DFT ersetzt, mit dem selben Ergebnis).
    Es gibt in der Regel 3 bis 4 Frequenzen mit fast der selben Amplitude.
    Wenn ich den Verlauf des Audiosignals anschaue, dann kann ich da klar eine Fundamentalperiode von in etwa 100 Hz erkennen.

    Das Audiosignal sieht so aus: http://dl.dropbox.com/u/16528181/audiosignal.gif
    und die dazugehörige FFT so: http://dl.dropbox.com/u/16528181/fft.gif

    Den Code für FFT und DFT hab ich übrigens von dieser Seite:
    http://paulbourke.net/miscellaneous/dft/

    Was muss ich ändern, damit ich die Fundamentalfrequenz berechnen kann?

    LG
    Michael

    0
  • R

    Du kannst ja einfach überprüfen, ob du oder die FFT recht hat. Filter einfach die Frequenzen raus, die du für unpassend hältst und schau, wie sich das im Vergleich zum Originalsignal auswirkt. Wenn die Komponenten wirklich so dominant sind, dann sollte sich das ja auch drastisch auswirken. Aber ich finde, dass die FFT schon hinkommen kann. Man sieht ja eindeutig, dass es da mehrere dominante Schwingungen gibt.

    btw. Für die FFT würde ich dir FFTW empfehlen. Ist schnell, gut und wird zB auch in Matlab verwendet.

    0
  • G

    Zeigst Du da von der FFT den Real- oder den Imaginaerteil? Oder den Absolutbetrag?

    0
  • P

    Gregor schrieb:

    Zeigst Du da von der FFT den Real- oder den Imaginaerteil? Oder den Absolutbetrag?

    Absolutbetrag.

    0

  • Erstmal solltest du das Signal mit einer Fensterfunktion multiplizieren bevor du die (F)FT machst.
    Hamming oder Blackman oder sowas.
    http://en.wikipedia.org/wiki/Window_function

    Dann könntest du probieren mehr als 800 Werte zu nehmen. Idealerweise auch eine Zweierpotenz, da so ne FFT damit am besten (schnellsten) funktioniert. Ansonsten, wenn du ganz bestimmte Frequenzwerte im Frequenzraum brauchst nimm 1600, 2400, 3200 statt 800.

    Und dann kann es natürlich sein, dass das, was im Zeit-Raum für das freie Auge wie eine dominate Frequenz aussieht, in Wirklichkeit die Überlagerung einer Grundschwindung mit ein paar Oberwellen ist. Wobei die Oberwellen durchaus gleich stark oder sogar stärker sein können als die Grundschwindung.

    Mach mal z.B. ein Array von 800 Werten...

    *25 = 1
    *75 = -1
    Alle anderen = 0

    ... und lass das durch deine FFT laufen. Dann siehst du was ich meine.

    0
  • P

    @hustbaer:

    Vielen Dank für deine Tipps!

    Ich hab mal versucht, mehr Werte zu nehmen. Ich hab die Samplerate auf 44100 Hz erhöht und nehm jetzt immer 2^14 (=16384) Werte.

    Leider ist es noch immer nicht ganz das gewünschte Ergebnis:
    Audiosignal->http://dl.dropbox.com/u/16528181/audio2.gif
    FFT->http://dl.dropbox.com/u/16528181/fft2.gif

    Jetzt werd ich mir noch das mit der Fensterfunktion anschauen...

    0
  • G

    1. Zeig mal den kompletten relevanten Code inklusive der Nutzung der FFT.

    2. Was für ein Audiosignal ist das? Ist das etwas spezielles?

    EDIT: Eine Sache, die Du wohl auf jeden Fall falsch machst, ist folgende: Du sagst, dass Du 800 Samples hast. Die FFT von der Seite da oben ist aber sehr einfach gehalten: Deshalb wird dort ausschließlich der Fall von Arrays mit 2^m Elementen behandelt. Da es kein ganzzahliges m mit 2^m=800 gibt, stimmt da schonmal etwas nicht. Aber es könnte auch mehr kaputt sein. Deswegen: Zeig Code!

    0
  • ?

    verdoppel das Fenster doch mal von 10 auf 20 ms.

    das entkoppelt die Fensterbreite von den 100 Hz.

    0
  • P

    Gregor schrieb:

    1. Zeig mal den kompletten relevanten Code inklusive der Nutzung der FFT.

    Hier der Code nach dem letzten Stand (mit jeweils 2^14 Samples):

    void ProcessData (signed short* AudioData, unsigned long SampleCount, fstream &Output)
{
    for (unsigned long l = 0; l < SampleCount; l++)
        Output << AudioData[l] << ","; // Audiodaten werden in die CSV-Datei geschrieben.

    Output << endl;

    double* RealVals = new double[16384];
    double* ImagVals = new double[16384];

    for (unsigned long l = 0; l < 16384; l++) {
        RealVals[l] = 0.0;
        ImagVals[l] = 0.0;
    }
    for (unsigned long l = 0; l < SampleCount; l++)
        RealVals[l] = (double)AudioData[l];

    // Die reellen Werte sind die Samples, der Rest wird mit 0 aufgefuellt,
    // die imaginaeren Werte sind alle 0.

    FFT(1, 14, RealVals, ImagVals);

    for (unsigned long l = 0; l < 8192; l++)
        Output << sqrt(RealVals[l] * RealVals[l] + ImagVals[l] * ImagVals[l]) << ",";
        // Der Absolutbetrag wird in die CSV-Datei geschrieben.

    Output << endl;

    delete [] RealVals;
    delete [] ImagVals;
}

/*
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
short FFT(short int dir,long m,double *x,double *y)
{
   long n,i,i1,j,k,i2,l,l1,l2;
   double c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for (i=0;i<m;i++) 
      n *= 2;

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for (i=0;i<n-1;i++) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   /* Compute the FFT */
   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }

   return(TRUE);
}

    Der Code für die FFT ist copy&paste von der Seite http://paulbourke.net/miscellaneous/dft/.

    0
  • G

    OK, die Nutzung sieht aus meiner Sicht auf den ersten Blick in Ordnung aus. Aber vielleicht solltest Du das alles mal explizit testen...

    Generier Dir ein Testsignal, bei dem das Resultat absolut klar ist. Einen Sinus oder so, bei dem Du die Frequenz kennst. Und dann guck mal, was rauskommt. Die Frequenz kann einmal mit der Länge Deines Samples synchronisiert sein und einmal nicht, nimm vielleicht einmal eine einzige Frequenz und dann nochmal 2. Wenn die Wellenlänge mit Deiner Samplelänge übereinstimmt, kannst Du statt dem Sinus auch mal den Cosinus nehmen und Dir Real- und Imaginärteil getrennt angucken.

    Ich meine, letztendlich ist es momentan so, dass Deine Vorstellung nicht mit dem Resultat übereinstimmt. Das kann 2 Gründe haben: Zum einen kann die Vorstellung falsch sein, zum anderen kann ein Fehler im Code sein. Deshalb sollte man versuchen, an Daten zu testen, bei denen man noch sehr genau nachvollziehen kann, was passiert. ...und dann geht man schrittweise zu komplizierteren Daten und schaut, ob man die Vorstellung dann mit der Zeit besser mit den Resultaten in Übereinstimmung bringen kann.

    EDIT: Du kannst auch mal ganz generelle Konsistenzchecks machen. Vergleich zum Beispiel mal die Energie im Zeitbereich mit der im Frequenzbereich...

    http://de.wikipedia.org/wiki/Parsevalsches_Theorem

    (unten unter "Anwendungen" gucken)

    0
  • B

    Hallo

    Das Feld der Tonhöhenbestimmung ist nicht ganz trivial. Ich hab mal ein Stimmgerät für den PC geschrieben (HarmonicTune) und da mit verschiedenen Algorithmen herum experimentiert.
    Hier ist mal eine erste Anlaufstelle:
    http://en.wikipedia.org/wiki/Pitch_detection_algorithm

    Aber einfach nur die FFT ausrechnen reicht meistens nicht. Das Problem ist dass das was du mit dem Auge aus dem Plot der FFT rauslesen kannst auch maschinell hinbekommen musst. Die wichtigsten Stichworte hierzu sind erstmal "Peakdetection" und "Noise Threshold". Aber vielleicht passt ja bei deinem Signal eher ein Ansatz aus dem Zeitbereich (wie bei mir z.B.).

    Ich würde dir empfehlen das Audiosignal erst mal ausgiebig mit Matlab (Scilab, Octave) und den da vorhandenen Plotmöglichkeiten genauestens zu untersuchen und dort den Algorithmus verfeinern bis du ihn schließlich in der Sprache deiner Wahl implementierst.

    Gruß

    0
  • K

    Ein paar Kommentare:

    DFT vs FFT: Ist dasselbe, zumindest rein mathematisch. Soll heißen, dass Du mit dem FFT-Algorithmus eine DFT berechnest. Die naive Implementierung einer DFT ist nur recht langsam. Der FFT-Algorithmus nutzt aus, dass sich die lineare Abbildung mit O(n^2) nicht-verschwindenden Elementen der dazugehörigen Matrix in O(log n) Matritzen mit jeweils O(n) nicht-verschwindenden Elementen faktorisieren lässt. Das heißt, man kommt auch mit O(n * log n) Rechenoperationen aus, wenn man weiß, welche das sind. 🙂

    In der Sprachkodierung/kompression (das, was jedes Handy macht zB), was ja auch typischerweise bei niedrigen Abtastraten wie 8000 Hz oder 16000 Hz stattfindet, ermitteln die Kodierer den Pitch soweit ich weiß über eine Autokorrelation. Zunächst wird das Signal durch einen Filter "entfärbt". Danach hört sich jeder Vokal fast gleich an wie ein helles Brummen. Danach kann man die Autokorrelation dieses Signals bestimmen (das geht auch per FFT/iFFT ganz fix). Das Ergebnis ist wieder ein Zeitsignal, welches bei 0 einen großen Peak hat. War das ursprüngliche Signal "voiced" (ein Vokal), so sieht die Autokorrelation etwa so aus:

    +    +    +
   +    |    |    |    +
   |    |    |    |    |
---+----+----+----+----+---    (Ein Kamm, eine Reihe von Peaks)
             0

              \__/
               |
          Pitch Period

    Hast Du dagegen einen Zischlaut von Dir gegeben, entsteht nach dem Entfärben kein Brummen sondern ein aperiodisches weißes Rauschen, dessen Autokorrelation so aussieht:

    +     
             |          
             |          
-------------+-------------    (Nur ein Peak und etwas "Dreck" drum herum)
             0

    Die Autokorrelation bekommst Du so:

    • Block fenstern. Die Blöcke können ruhig überlappen, sofern Du Dir die Verzögerung erlauben kannst. Bzgl Blockgröße: Es kommt hier auch darauf an, wie groß die größte Pitch-Period ist, die Du erfassen willst.
    • FFT berechnen
    • Betrag-Quadrat ausrechnen (je Koeffizient)
      neu = real(alt)*real(alt)+imag(alt)*imag(alt)
    • inverse FFT
    • An Offset 0 ist auch die 0 (siehe oben). Da das Ergebnis reell und symmetrisch ist, brauchst Du dich nur um die rechte Hälfte zu kümmern (Offsets 0...blockgroesse/2)
    • Die Peakhöhe bei Offset 0 ist proportional zur Leistung des Signals

    Das funktioniert natürlich nur, wenn es nur eine einzige Stimme bzw Note eines Instruments war. Legst Du mehrere Stimmen übereinander, wird die Sache relativ kompliziert. Hier weiß ich auch nicht, ob dann die Korrelationsmethode immer noch brauchbar ist. Du kannst ja mal die Firma Melodyne fragen, wie sie ihr "direct note access"-Feature implementiert haben. Dieses Geschäftsgeheimnis werden sie Dir aber wahrscheinlich nicht verraten.

    Wenn ich Code wie diesen hier sehe:

    double* p = new double[16384]; 
  double* q = new double[16384]; 
  delete[] q;
  delete[] p;

    kommen mir folgende Gedanken:

    • Der Autor weiß RAII nicht zu schätzen. Sonst hätte er vector<double> benutzt.
    • Magic Numbers sind doooof 😉
    0
Anmelden zum Antworten