Fourier-Transformation



  • Hallo Leute,
    ich benötige hilfe in Fragen Fourier-Transformation!
    Ich habe einen Funktionsverlauf durch Wertepaare gegeben, z.B: "zeit = 0.305s", "Strom = 125A". nächstes Wertepaar "zeit = 0.306s", Strom=127A" usw. Das alles ergibt dann einen nicht ganz Sinusförmigen Verlauf! Ich benötige jetzt eine Möglichkeit diesen Verlauf mittels Fouriertransformation als Frequenzspektrum darzustellen.
    Kann mir jemand eine Möglichkeit nennen, wie das umzusetzen geht und ob schon Quelltexte dafür existieren?

    Gruß
    Ronald



  • Hallo,

    ich habe da auch schon danach gesucht.

    Schau mal bei diesem Link nach, da solltest Du finden, was Du suchst:
    http://www.library.cornell.edu/nr/cbookcpdf.html (Numerical Recipes in C)

    Dort findest Du eine normale FFT(12.3) und eine Leistungs-FFT(13.4)



  • Dieser Thread wurde von Moderator/in HumeSikkins aus dem Forum C++ in das Forum Rund um die Programmierung verschoben.

    Im Zweifelsfall bitte auch folgende Hinweise beachten:
    C/C++ Forum :: FAQ - Sonstiges :: Wohin mit meiner Frage?

    Dieses Posting wurde automatisch erzeugt.



  • Mmoment, hatten wir nicht eins zu eins diese Frage gerade schonmal?



  • ja, hatten wir, ich habe nicht gewußt das die beiträge verschoben werden und die frage in zwei verschiedenen foren gestellt!

    mein problem ist jetzt folgendes, ich habe einen quellcode, und weiß nun nicht an welcher stelle er mir die frequenz und den dazugehörigen betrag ausgibt!

    kann mir bitte jemand helfen!



  • Er gibt dir ger keine Frequenzen. Die kannst du selber errechnen. Es geht von DC (0Hz) bis Nyquist (Samplerate/2).
    Die Amplitude erhälst du durch sqrt (cos_anteil*cos_anteil+sin_anteil*sin_anteil).



  • Jetzt muss ich aber auch noch mal etwas fragen...

    Ich bin gerade dabei einen Analyser per FFT zu programmieren.

    Momentan verwende ich hierfür die Power Spektrum Analyse von Numerical Recipes.
    Diese Funktion erledigt das Windowing (& Korrektur) und errechnet mir die Leistung aus meiner Frequenz.

    Nur ist es so, dass ich für die Funktion ca. zwei Perioden übergeben muss. Meine FFT-Analyse soll jedoch bis 20 Hz gehen. Dies bedeutet ich benötige einen Arrey mit einer Größe von 2048.

    Von Echtzeit kann man an dieser Stelle nun niht mehr reden. Auch gibt mir die Funktion feste Frequenzwerte zurück. Diese ergeben sich aus:

    Frequenzschritt=1/Sampledauer(Zeit) ==> bei 2048 sind dies Schritte von 10,8 Hz.

    Gibt es da eine bessere und schnellere Methode?
    Ich sollte von 20-20kHz das Spektrum wiedergeben können. Habe jedoch kein periodisches Signal, was bedeutet, dass ich eine Fensterung benötige.

    danke

    Hier mein Code:

    void spctrm(short int fp[], float p[], int m, int k, int ovrlap)
    //Reads data from input stream specified by file pointer fp and returns as p[j]
    //the data's power (mean square amplitude) at frequency (j-1)/(2*m) cycles per
    //gridpoint, for j=1,2,...,m, based on (2*k+1)*m data points (if ovrlap is set
    //true (1)) or 4*k*m data points (if ovrlap is set false (0)).
    //The number of segments of the data is 2*k in both cases: The routine calls
    //four1 k times, each call with 2 partitions each of 2*m real data points.
    {
            int mm,m44,m43,m4,kk,joffn,joff,j2,j,i;
            float w,facp,facm,*w1,*w2,sumw=0.0,den=0.0;
            mm=m+m; //Useful factors.
            m43=(m4=mm+mm)+3;
            m44=m43+1;
            w1=vector(1,m4);
            w2=vector(1,m);
            facm=m;
            facp=1.0/m;
            for (j=1;j<=mm;j++) sumw += SQR(WINDOW(j,facm,facp));
    //        Accumulate the squared sum of the weights.
            for (j=1;j<=m;j++) p[j]=0.0; //Initialize the spectrum to zero.
            if (ovrlap){ i=0;//Initialize the "save" half-buffer.
               for (j=1;j<=m;j++){
                  w2[j]=fp[i]; i++;
               }
            }
            for (kk=1;kk<=k;kk++) {
            //Loop over data set segments in groups of two.
               for (joff = -1;joff<=0;joff++) {
                //Get two complete segments into workspace.
                 if (ovrlap) {
                    for (j=1;j<=m;j++) w1[joff+j+j]=w2[j];
                    i=0;
                    for (j=1;j<=m;j++) {
                        w2[j]=fp[i];
                        i++;
                        }
                    joffn=joff+mm;
                    for (j=1;j<=m;j++) w1[joffn+j+j]=w2[j];
                 }else {i=0;
                    for (j=joff+2;j<=m4;j+=2){
                       w1[j]=fp[i];
                       i++;
                       }
                 }
               }
               for (j=1;j<=mm;j++) { //Apply the window to the data.
                  j2=j+j;
                  w=WINDOW(j,facm,facp);
                  w1[j2] *= w;
                  w1[j2-1] *= w;
               }
               four1(w1,mm,1); //Fourier transform the windowed data.
               p[1] += (SQR(w1[1])+SQR(w1[2])); //Sum results into previous segments.
               for (j=2;j<=m;j++) {
                  j2=j+j;
                  p[j] += (SQR(w1[j2])+SQR(w1[j2-1])
                       +SQR(w1[m44-j2])+SQR(w1[m43-j2]));
                  }
                  den += sumw;
            }
            den *= m4; //Correct normalization.
            for (j=1;j<=m;j++) p[j] /= den; //Normalize the output.
            free_vector(w2,1,m);
            free_vector(w1,1,m4);
    }
    

    ...und der von der FFT

    void four1(float data[], unsigned long nn, int isign)
    {
    	unsigned long n, mmax, m, j, istep, i;
    	double wtemp, wr, wpr, wpi, wi, theta;
    	float tempr, tempi;
    
    	n = nn << 1;			
    	j = 1;
    	for (i = 1; i < n; i += 2) {
    		if (j > i) {
    			SWAP(data[j], data[i]);
    			SWAP(data[j+1], data[i+1]);
    		}
    		m=n >> 1;
    		while (m >= 2 && j > m) {
    			j -= m;
    			m >>= 1;
    		}/* while */
    		j += m;
    	}/* for */
    	mmax = 2;
    
    	while (n > mmax) {
    		istep = mmax << 1;
    		theta = isign * (6.28318530717959 / mmax);
    		wtemp = sin(0.5 * theta);
    		wpr = -2.0 * wtemp * wtemp;
    		wpi = sin(theta);
    		wr = 1.0;
    		wi = 0.0;
    		for (m = 1; m < mmax; m += 2) {
    			for (i = m; i <= n; i += istep) {
    				j = i + mmax;
    			   tempr = (float)(wr * data[j] - wi * data[j+1]);	
    			   tempi = (float)(wr * data[j+1] + wi * data[j]);	
    				data[j] = data[i] - tempr;
    				data[j+1] = data[i+1] - tempi;
    				data[i] += tempr;
    				data[i+1] += tempi;
    			}
    			wr = (wtemp = wr) * wpr - wi * wpi + wr;
    			wi = wi * wpr + wtemp * wpi + wi;
    		}
    		mmax = istep;
    	}
    }
    
    /* (C) Copr. 1986-92 Numerical Recipes Software ) */
    


  • Nur ist es so, dass ich für die Funktion ca. zwei Perioden übergeben muss. Meine FFT-Analyse soll jedoch bis 20 Hz gehen. Dies bedeutet ich benötige einen Arrey mit einer Größe von 2048.

    Ohne nähere Angaben kann ich nichts dazu sagen.

    Von Echtzeit kann man an dieser Stelle nun niht mehr reden. Auch gibt mir die Funktion feste Frequenzwerte zurück.

    Ach.

    Diese ergeben sich aus:
    Frequenzschritt=1/Sampledauer(Zeit) ==> bei 2048 sind dies Schritte von 10,8 Hz.

    Aha, du scheinst eine Samplerate von 44100Hz zu haben, korrekt? 10,8Hz * 2048=22118,4, käme also ungefähr hin.

    Gibt es da eine bessere und schnellere Methode?
    Ich sollte von 20-20kHz das Spektrum wiedergeben können. Habe jedoch kein periodisches Signal, was bedeutet, dass ich eine Fensterung benötige.

    1024er FFT? Das zweite Bin ligt bei 21,5Hz.



  • OK, werde mal kurz meine Daten zusammenfassen:

    Samplerate von 44.1 kHz mit 16 Bit Auflösung
    Frequenzgang 20Hz-20kHz
    Aufnahmezeit: 0,09288 Sekunden ==> 2048 Werte

    1024er FFT? Das zweite Bin ligt bei 21,5Hz.

    Wenn ich es mit einem Arrey mit 1024 mache, dann ist die Funktion schon deutlich schneller, jedoch würde dies bedeuten, dass ich mit einer einzigen Periode arbeiten würde. Liefert diese mir ein zuverlässiges Ergebnis?

    Ein weiteres Problem ist, dass die FFT mir bei jedem Zyklus einen anderen Wert liefert, selbst bei gleich bleibendem Signal.

    Gibt es noch weitere Routinen, die schneller und genauer arbeiten?



  • Er gibt dir ger keine Frequenzen. Die kannst du selber errechnen. Es geht von DC (0Hz) bis Nyquist (Samplerate/2).
    Die Amplitude erhälst du durch sqrt (cos_anteil*cos_anteil+sin_anteil*sin_anteil).

    @Helium:
    So wie ich Deine Antwort verstehe, kann ich hier jeweils die Frequenz, die ich haben möchte berechnen. z.B. 20 Hz, 50 Hz, 100 Hz, 1kHz,...
    Also jede X-beliebige und bin nicht auf die Schritte meiner Funktion angewiesen?

    Wie sieht die Funktion aus? Hat die auch eine Fensterung usw.?



  • So wie ich Deine Antwort verstehe, kann ich hier jeweils die Frequenz, die ich haben möchte berechnen. z.B. 20 Hz, 50 Hz, 100 Hz, 1kHz,...
    Also jede X-beliebige und bin nicht auf die Schritte meiner Funktion angewiesen?

    Ups, nein, so war das nicht gemeint. RonaldoS sagte er wisse nicht, wo die Frequenz ausgegeben werde. Die ist aber allein durch die Samplerate festgelegt.



  • ...dann stimmt das zumindest bei mir bisher... danke 🙂

    D.h. ich sollte schauen, dass ich die Funktion oben so zum laufen bekomme, wenn ich sowohl die Fensterung als auch die Leistung benötige.

    Nun nochmals zu 1024-FFT. Meinst Du es genügt der Routine wenn ich nur eine Periode übergebe?

    Macht es Sinn, die FFT-Routine als Thread zu programmieren?



  • hallo "Helium"

    Ich habe noch ein kleines Verständnissproblem! Ich übergebe der FFT folgendes:

    rfft(tau, y[], synt);
    

    wobei tau: Zweierlogarithmus der Anzahl der Funktionswerte oder, anders ausgedrueckt: Die Anzahl der Funktionswerte ist N = 2^tau.
    wobei y[]: die Funktionswerte enthält.

    tau ist in meinem Fall 16, das bedeuted, ich habe 65536 Werte also auch mein Feld y[] beinhaltet 65536 Werte.

    Mein Problem jetzt:
    diese Werte stammen von einer oberschwingungsbehafteten Sinusfunktion, genauer einem Fahrzeugstrom eines elektrischen Triebfahrzeuges, mit entweder 50 Hz oder 16 2/3 Hz. wobei genau nur eine Periode betrachtet wird!
    ich benötige nun alle in dieser einen Periode vorkommenden Frequenzen und dazugehörigen Amplituden.

    ich weiß nicht an welcher stelle welches ergebniss geliefert wird! kannst du mir das sagen?
    Hier ist mein quellcode:

    FILE *testaus;
      int  N,           // Anzahl der Funktionswerte
           Nd2,         // N / 2
           Nd4,         // N / 4
           sigma,       // die Spiegelung der tau-1 Ziffern der Binaer-
                        // darstellung von j an ihrer Mitte
           min_n,       // 2 ^ (tau - 1 - n)
           n_min_0,     // 2 ^ n
           n_min_1,     // 2 ^ (n - 1)
           ind1,        // Zwischenspeicher fuer haeufig verwendete
           ind2,        // komplizierte Indexausdruecke
           k, j,        // Laufvariablen
           n, l;        // Laufvariablen
    double faktor,      // Normierungsfaktor 2/N fuer die Fourier-
                        // koeffizienten im Fall der Analyse, sonst 1
                        // fuer die Synthese der Funktionswerte
           arg,         // Winkel von (wr,wi)
           arg_m,       // Winkel der N. Einheitswurzel
           arg_md2,     // Winkel der (N/2). Einheitwurzel
           ew_r,        // Real- und Imaginaerteil
           ew_i,        // der N. Einheitswurzel
           eps_r,       // Real- und Imaginaerteil von
           eps_i,       // (N. Einheitswurzel) ^ k   bzw. von
                        // ((N/2). Einheitswurzel) ^ (l * 2^min_n)
           ur, ui,      // Zwischenspeicher fuer eine komplexe Zahl
           wr, wi,      // Real- bzw. Imaginaerteil von
                        // ((N/2). Einheitswurzel) ^ (2^min_n)+7
           rett,        // Zwischenspeicher fuer den alten Wert von eps_r
           yhilf,       // Hilfsvariable
           hilf1,       // Hilfsvariable
           hilf2,       // Hilfsvariable
           hilf3,       // Hilfsvariable
           hilf4;       // Hilfsvariable
    //---------------------------------------------------------------------------
      testaus= fopen("f:\\tmp\\fft.txt", "wt");
      if (tau < 1)                        /* zu kleiner Wert fuer tau?    */
        return 1;
    
      if (tau > 8 * sizeof(int) - 2)      /* zu grosser Wert fuer tau?    */
        return 2;                         /* (Ueberlaufgefahr bei 2^tau!) */
    
      N       = 1 << tau;    // N  =  2 hoch tau
      Nd2     = N / 2;
      Nd4     = Nd2 / 2;
      faktor  = ONE / Nd2;
      arg_md2 = TWO * M_PI * faktor;
      arg_m   = HALF * arg_md2;
    
      if (synthese == 0)
        faktor = ONE;
    
      if (synthese == 0)       // die reellen Daten zusammenfassen zur
      {                        // Ausfuehrung einer FFT halber Laenge
        yhilf =  y[1];
        y[1]  =  y[0] - y[N - 1];
        y[0]  += y[N - 1];
    
        ew_r = COS(arg_m);       // (ew_r,ew_i) = N. Einheitswurzel
        ew_i = SIN(arg_m);
        eps_r = ONE;             // (eps_r,eps_i) = (N. Einheitswurzel)^k
        eps_i = ZERO;
    
        for (k = 1; k < Nd4; k++)
        {
          ind1 = 2 * k;
          ind2 = N - ind1;
          rett = eps_r;
          eps_r = rett * ew_r - eps_i * ew_i;
          eps_i = rett * ew_i + eps_i * ew_r;
          hilf1 = HALF * (eps_r * (yhilf   - y[ind2 - 1]) +
                         eps_i * (y[ind1] + y[ind2]));
          hilf2 = HALF * (eps_i * (yhilf   - y[ind2 - 1]) -
                          eps_r * (y[ind1] + y[ind2]));
          hilf3 = HALF * (yhilf   + y[ind2 - 1]);
          hilf4 = HALF * (y[ind1] - y[ind2]);
          yhilf = y[ind1 + 1];
          y[ind1]     = hilf3 - hilf2;
          y[ind1 + 1] = hilf1 - hilf4;
          y[ind2]     = hilf2 + hilf3;
          y[ind2 + 1] = hilf1 + hilf4;
        }
        y[Nd2 + 1] = y[Nd2];
        y[Nd2]     = yhilf;
      }
    
      for (j = 0; j < Nd2; j++)                       // mit der Bit-
      {                                               // Umkehrfunktion
        for (k = j, n = 1, sigma = 0; n < tau; n++)   // umspeichern, im
          sigma <<= 1,                                // Fall der Analyse
          sigma |=  k & 1,                            // gleichzeitig
          k     >>= 1;                                // normieren
        if (j <= sigma)
          ind1 = 2 * j,
          ind2 = 2 * sigma,
          ur = y[ind1],
          ui = y[ind1 + 1],
          y[ind1]     = y[ind2]     * faktor,
          y[ind1 + 1] = y[ind2 + 1] * faktor,
          y[ind2]     = ur          * faktor,
          y[ind2 + 1] = ui          * faktor;
      }
    
      for (min_n = Nd2, n_min_1 = 1, n = 1; n < tau; n++)   // die FFT
      {                                                     // halber
        min_n   /= 2;                                       // Laenge
        n_min_0 =  2 * n_min_1;                             // ausfuehren
        arg =  arg_md2 * min_n;
        wr  =  COS(arg);
        wi  =  synthese ? ONE : -ONE;
        wi  *= SIN(arg);
        eps_r = ONE; //(eps_r,eps_i)=((N/2). Einheitswurzel) ^(l * 2^min_n)
        eps_i = ZERO;
    
        for (l = 0; l < n_min_1; l++)
        {
          for (j = 0; j <= Nd2 - n_min_0; j += n_min_0)
          {
            ind1 = (j + l) * 2;
            ind2 = ind1 + n_min_0;
            ur = y[ind2] * eps_r - y[ind2 + 1] * eps_i;
            ui = y[ind2] * eps_i + y[ind2 + 1] * eps_r;
            y[ind2]     =  y[ind1]     - ur;
            y[ind2 + 1] =  y[ind1 + 1] - ui;
            y[ind1]     += ur;
            y[ind1 + 1] += ui;
          }
          rett  = eps_r;
          eps_r = rett * wr - eps_i * wi;
          eps_i = rett * wi + eps_i * wr;
    		fprintf(testaus,"%12.6f\n", y[l]);
        }
        n_min_1 = n_min_0;
      }
    


  • Nochmal was zum thema echtzeit:

    Also soweit von echtzeit ist ganze nicht weg, bzw wenn man genung dafür bereit stellt geht schon in echtzeit, ich habs mal ausprobiert, bei einer 2048 FFT dürfte man bei 44100khz max 47 ms benötigen (bei mono). So, bei einem von mir selbst verfassten algorithmus (der jetzt endlich funktioniert :-), zumindest der rekursive), brauche ich dafür 31 ms, also noch im echtzeitfenster. Jetzt muss man ja auch noch sehen, das der algo von mir is und optimierungstechnisch unter aller sau ist.
    Aber man wird wohl einen rechner brauchen der nicht zu langsam ist. (bei mir Notebook mit 2,67 Ghz).

    Aber prinzipiell kann man das unter umständen schon noch in echtzeit lösen, denk ich .



  • Bei mir braucht die Funktion ungefähr eine halbe Sekunde bei einer 2048 FFT. Habe einen Athlon 1400+ Mobil Prozessor.

    Mir scheint Deine Funktion arbeitet etwas schneller... 😞



  • Flow_cplus schrieb:

    Nochmal was zum thema echtzeit:

    Also soweit von echtzeit ist ganze nicht weg,

    Und wie schaut es mit der echtzeit-Darstellung aus ? 😉



  • Logger schrieb:

    Flow_cplus schrieb:

    Nochmal was zum thema echtzeit:

    Also soweit von echtzeit ist ganze nicht weg,

    Und wie schaut es mit der echtzeit-Darstellung aus ? 😉

    naja ich denke schon das ich das auf meinem rechner noch so ungefähr hinbekomme, wenn ich einen ordentlichen alogrithmus nehme, zur Darstellung opengl. und hier und da noch etwas Abstriche mache ...
    im moment nutzte man ich das ganze mit double, da ich aber nur 16 bit brauch, und bei echtzeit muss die ausgabe ja auch nicht so genau sein, weil es ja nur um optische 'Anschauungsgenauigkeit' geht. (dh die Ausgangswerte dürfen in einem bereich von 0 - max ~1000 pixel liegen)
    Messtechnische Erfassungen, wird man hier wohl nicht in echtzeit machen.

    Was anders wäre es wenn man z.B Filter oder solche sachen mit FFT und FFT^-1 lösen wollte, da ist dann wohl in echtzeit nix zu machen.



  • Richtig, wöfur brauchst du überhaupt 2048 Werte, du bekommst ja eh maximal 1000+x auf deinen Bildschirm (1024,1280...)



  • Ich denke mal, dies liegt dann wohl an den Frequenzschritten...



  • naja

    das ist das Problem des FFt algorithmus, wenn ich es soweit verstanden habe. die Anzahl der Eingangswerte = die anzahl der Ausgangswerte, dh. bei 1024 samples bekommt man 1024 spektrallinien im definierten freqenzabstand von fs/N
    Mit orginal dft summenformel ist das nicht so, allerdings kommt man da auf N*N rechenschritte, während die fft nur N*ld N.
    Ich überleg bloß grad wie es ist wenn ich nur ca. 15 spektrallienen haben will, dann sollte dft und fft ungefähr ähnlich schnell sein, nur kann ich mir dft aussuchen welche spektrallienien ich ausrechne, und wenn ich dann die auflösung soweit reduziere, das es langt, zb von den 16 bit audiosignal wenn nur die ersten(höherwertigen) 8 bit nehme, dann käme ich auf eine Anzeige immer noch genung Pixel Höhe, das würde mir langen, ich will ja nur ne relative kleine Anzeige.
    Und wieweit es wichtig is bei optischen anzeige alles samples durch die dft zu jagen, weiß ich auch nicht vielleicht lang es auch nur einen teil zu nehmen ohne dass man die sache zu arg verfälscht.. ?


Log in to reply