Schnelle Fouriertransformation
Übersicht
BlitzMax, BlitzMax NG Codearchiv & Module
cathodeBetreff: Schnelle Fouriertransformation |
 Mi, Feb 20, 2008 18:26Antworten mit Zitat 
|
|
|---|---|---|
|
Hallo Blitz*-Coder,
hier ist eine Implementierung des Cooley-Tukey-FFT-Algorithmus. Es gibt noch ein Problem mit dem Baum: wenn er außerhalb der fft-Funktion initialisiert wird, dann bleiben die Werte in den FT-Arrays der Knoten, so dass bei wiederholter Anwendung der FFT die Werte aufsummieren. Ich habe leider keine Ahnung, woran das liegen könnte. Code: [AUSKLAPPEN] Type node Field parent:Int 'index of parent node Field length:Int 'length of corresponding subsequence Field offset: Int 'offset within main sequence Field ft_re:Float[] 'node's fourier transform Field ft_im:Float[] End Type Type tree Field nodes:node[] Field lastentry:Int End Type Type spectrum Field a:Float[] 'real part Field b:Float[] 'imaginary part End Type Function GenerateTree:tree(l:Int) 'computational tree setup Local numnodes = 2*l-1 Local t:tree = New tree t.nodes = New node[numnodes] If (l > 0) 'init root node t.nodes[0] = New node t.nodes[0].ft_re = New Float[l] t.nodes[0].ft_im = New Float[l] t.nodes[0].parent = -1 t.nodes[0].length = l t.nodes[0].offset = 0 Local q_even:Int = 0 'length of even subsequence Local q_odd:Int = 0 'length of odd subsequence Local j:Int = 1 Local i:Int = 0 Local _i:Int = 0 Local dist:Int = 1 Repeat q = t.nodes[i].length If (q > 1) q_odd = q shr 1 q_even = q-q_odd t.nodes[j] = New node t.nodes[j+1] = New node t.nodes[j].ft_re = New Float[q_even] t.nodes[j].ft_im = New Float[q_even] t.nodes[j+1].ft_re = New Float[q_odd] t.nodes[j+1].ft_im = New Float[q_odd] t.nodes[j].parent = i t.nodes[j].length = q_even t.nodes[j].offset = t.nodes[i].offset If (_i=dist) dist = dist + dist _i=0 End If _i=_i+1 t.nodes[j+1].parent = i t.nodes[j+1].length = q_odd t.nodes[j+1].offset = t.nodes[i].offset + dist j = j + 2 End If i = i + 1 Until (q = 1) t.lastentry = j-1 End If Return t End Function Function fft:spectrum(s:Float[])', t:tree) Local spec:spectrum = New spectrum Local t:tree = GenerateTree(s.length) spec.a = New Float[s.length] spec.b = New Float[s.length] Local par:Int Local even_re:Float, odd_re:Float, even_im:Float, odd_im:Float Local f_re:Float, f_im:Float Local of_re:Float, of_im:Float Local ooff_c:Int, eoff_c:Int Local arg:Float Local _N:Float Local h:Float Local lastnode:Int = t.lastentry For i = lastnode To 2 Step -2 par = t.nodes[i].parent 'index of parent node N = t.nodes[i].length 'length of subsequence _N = N h = 180/_N ooff_c = t.nodes[i].offset 'offsets of subsequences eoff_c = t.nodes[i-1].offset If (N = 1) t.nodes[i-1].ft_re[0] = s[eoff_c] t.nodes[i].ft_re[0] = s[ooff_c] End If For j = 0 To N-1 even_re = t.nodes[i-1].ft_re[j] 'real and imaginary parts even_im = t.nodes[i-1].ft_im[j] 'of subsequence fourier transforms odd_re = t.nodes[i].ft_re[j] odd_im = t.nodes[i].ft_im[j] _j = j + N arg = h * j f_re = Cos(arg) 'phase factor f_im = -Sin(arg) 'of odd subseq. of_re = f_re*odd_re - f_im*odd_im of_im = f_re*odd_im + f_im*odd_re 'if i=2 then write values into final-spectrum array If (i>2) t.nodes[par].ft_re[j] = t.nodes[par].ft_re[j] + even_re + of_re t.nodes[par].ft_im[j] = t.nodes[par].ft_im[j] + even_im + of_im t.nodes[par].ft_re[_j] = t.nodes[par].ft_re[_j] + even_re - of_re t.nodes[par].ft_im[_j] = t.nodes[par].ft_im[_j] + even_im - of_im Else spec.a[j] = spec.a[j] + even_re + of_re spec.b[j] = spec.b[j] + even_im + of_im spec.a[_j] = spec.a[_j] + even_re - of_re spec.b[_j] = spec.b[_j] + even_im - of_im EndIf Next Next Return spec End Function 'Beispiel Graphics 512,384 Local s:Float[512] s[10] = 100 Local spec:spectrum = fft(s) For k = 0 To s.length-1 Plot k, 192-spec.a[k] Next Flip Input gruß, cathode |
||
|
Markus2 |
Mi, Feb 20, 2008 20:42 Antworten mit Zitat
|
|---|---|---|
|
Kannste bitte auch mal kurz erkären was es macht und
was man damit alles tolles machen bzw wofür das gut ist ? Gehört habe ich mal das man es für Spracherkennung einsetzen kann aber genau wie hab ich nicht so kapiert . |
||
cathode |
Mi, Feb 20, 2008 21:33 Antworten mit Zitat
|
|
|---|---|---|
|
hallo markus.
es geht im prinzip darum, eine zahlenfolge in frequenzkomponenten zu zerlegen. mathematisch gesehen lässt sich jede funktion als unendliche summe von "schwingungen" darstellen. eine fouriertransformation berechnet sozusagen für jede frequenz, wie stark diese in einer funktion enthalten ist. zahlenfolgen sind nichts anderes als diskrete funktionen und können ebenfalls in frequenzen zerlegt werden. mit der fouriertransformation kann man zum beispiel das frequenzspektrum von tönen gewinnen. wer winamp benutzt kennt sicherlich das kleine fenster mit den balken drin. anwendungen der FT gibt es viele, als beispiel seien nur frequenzfilter genannt. jetzt zur schnellen FT: ein "naiver" FT-algorithmus benötigt zur verarbeitung einer folge mit N elementen N^2 rechenoperationen, so dass bei großen fouriertransformationen die rechenzeit nicht mehr akzeptabel ist. ein "schneller" algorithmus nutzt einige mathematische eigenschaften der FT aus, um den rechenaufwand auf NlogN operationen zu senken. gruß, cathode [edit] zur benutzung: der algorithmus schluckt grundsätzlich nur samples der länge 2^n. vor der transformation muss GenerateTree(l) aufgerufen werden (wird in fft() automatisch gemacht), l ist dabei die länge der sequenz. die zu transformierende sequenz ist ein eindimensionales array des typs float. der rückgabetyp ist spectrum und enthält für den real- und imaginärteil jeweils ein array der länge l. |
||
|
Vertex |
Fr, Feb 22, 2008 0:57 Antworten mit Zitat
|
|---|---|---|
|
Die Balken die man aus der Spektrumanalyse von Winamp bspw. kennt, müssten noch in reale Werte "umgewandelt" werden mit Sqrt(rl²+im²).
Habe selber noch einen alten Code dafür herausgekramt: Code: [AUSKLAPPEN] Function IntFFT(In:Int[], Out:Float[], Length:Int)
Local Temp : Float[,], .. Index : Int, .. Radix : Float[2] Temp = New Float[Length, 2] For Index = 0 Until Length Temp[Index, 0] = In[Index] Temp[Index, 1] = 0.0 Next Radix[0] = Cos(360.0/Length) Radix[1] = Sin(360.0/Length) DoFFT(Temp, Length, 0, Radix) For Index = 0 Until Length Out[Index] = Sqr(Temp[Index, 0]*Temp[Index, 0] + .. Temp[Index, 1]*Temp[Index, 1]) Next End Function Function DoFFT(Data:Float[,], Length:Int, Start:Int, Radix:Float[]) Local Index : Int, .. Middle : Int, .. Z : Float[2], .. V : Float[2], .. H : Float[2], .. TempA : Float, .. TempB : Float If Length =< 1 Then Return Middle = Length Shr 1 Z[0] = 1.0 Z[1] = 0.0 For Index = Start To Start + Middle - 1 H[0] = Data[Index, 0] - Data[Index + Middle, 0] H[1] = Data[Index, 1] - Data[Index + Middle, 1] Data[Index, 0] :+ Data[Index + Middle, 0] Data[Index, 1] :+ Data[Index + Middle, 1] Data[Index + Middle, 0] = H[0]*Z[0] - H[1]*Z[1] Data[Index + Middle, 1] = H[0]*Z[1] + H[1]*Z[0] TempA = Z[0] TempB = Z[1] Z[0] = TempA*Radix[0] - TempB*Radix[1] Z[1] = TempA*Radix[1] + TempB*Radix[0] Next ' v=w*w V[0] = Radix[0]*Radix[0] - Radix[1]*Radix[1] V[1] = Radix[0]*Radix[1] + Radix[1]*Radix[0] DoFFT(Data, Middle, Start, v) DoFFT(Data, Middle, Start + Middle, v) Shuffle(Data, Length, Start) End Function Function Shuffle(Data:Float[,], Length:Int, Start:Int) Local Index : Int, .. Middle : Int, .. Temp : Float[,] Middle = Length Shr 1 Temp = New Float[Middle, 2] MemCopy(Temp, Byte Ptr(Data) + Start*8, Middle*4) For Index = 0 Until Middle Data[Start + Index + Index + 1, 0] = Data[Start + Index + Middle, 0] Data[Start + Index + Index + 1, 1] = Data[Start + Index + Middle, 1] Next For Index = 0 Until Middle Data[Start + Index + Index, 0] = Temp[Index, 0] Data[Start + Index + Index, 1] = Temp[Index, 1] Next End Function Damit habe ich früher Sampledaten(sind Integerwerte) vom Mikrofon analysiert. Beim Klavier konnte man bspw. dann im Graphen genau ablesen, welche Frequenz gerade ertönt und daraus die Note ermitteln. Die Theorie ist etwas lange her, aber wie ich es verstanden habe, basiert FFT auf komplexen Zahlen. Um als Eingabe Integer- bzw. Realwerte zu benutzen, müssen diese in komplexe Werte gewandelt werden. Dazu wird der Realanteil der Komplexen Zahl mit dem Integer- resp. Realwert belegt. Als Ausgabe hat der Algorithmus wiederum komplexe Werte die dann über den Betrag von Real- und Imaginäranteil in so ein Balkendiagramm gebracht werden kann. Der erste Wert ist die Gleichspannung. Weiterhin spiegeln sich die Werte ab Länge/2. Es lassen sich nur Frequenzen mit einer Frequenz von Samplerate/2 (Nyquest Frequenz?!) ermitteln. Ein Equalizer ist übrigens die Umkehrfunktion. mfg olli |
||
| vertex.dreamfall.at | GitHub | ||
Übersicht
BlitzMax, BlitzMax NG Codearchiv & Module
Powered by phpBB © 2001 - 2006, phpBB Group