Audiovisualisierung mit OpenAL

Übersicht BlitzMax, BlitzMax NG Allgemein

Neue Antwort erstellen

 

timmeTheOnly

Betreff: Audiovisualisierung mit OpenAL

 BeitragMi, Okt 19, 2011 21:08
Antworten mit Zitat
Benutzer-Profile anzeigen
Guten Abend,

ich habe eine Frage an euch und dafür ersteinmal eine kleine Geschichte Wink

Als wir neulich gefeiert haben, hab ich meinen Projektor als Lichtshow verwendet und muss sagen, dass das richtig gut kommt. Leider fing das Programm mit dem ich arbeitete (BeamerSoundToLight) nach einer gewissen Zeit an mit Fehlermeldungen um sich zu werfen (wohl inkompatibel mit Windows 7).

Weil mir die Möglichkeiten des Programms zum Teil zu wenig waren (keine Farbwechsel, Flashfiles, ...) bin ich auf die Idee gekommen sowas selbst zu machen.

Jetzt sitze ich hier und muss feststellen, dass ich keine Ahnung habe, wie ich das anstellen soll.

Die Stereomix Daten habe ich bereits in einem TAudioSample
(An dieser Stelle danke an Midimaster für sein Tutorial),
allerdings weiß ich nicht, wie ich diese sinnvoll verwenden kann.

Wie bekomm ich die Daten in eine Form, in der ich sie nicht nur abspielen kann?
Außerdem hab ich gelesen, dass man beim Frequenzen trennen um eine FFT nicht herum kommt und dafür reichen meine Mathekenntnisse (11. Klasse) einfach nicht aus.

Ich hoffe ihr könnt mir helfen Wink

Gruß,
Tim

Vertex

 BeitragDi, Okt 25, 2011 20:52
Antworten mit Zitat
Benutzer-Profile anzeigen
Hallo,

genau, was visualisiert wird sind die Frequenzanteile, die Du mit Hilfe der Fourieranalyse ermitteln kannst. Da die Analyse auf diskreten Werten basiert, musst Du die diskrete Fourieranalyse (DFT) durchführen. FFT (Fast Fourier Transformation) ist ein Alagorithmus, der Symmetrien ausnutzt und zur Berechnung der DFT nur O(n*log(n)) benötigt.

FFT: Du gibst die Messwerte einer beliebigen Funktion (= die aufgenommenen Samples in einem Zeitbereich) rein, und heraus kommen Frequenzanteile. Es wird vesucht die Funktion durch Überlagerung (Aufaddieren) von Sinus- und Kosinus-Schwingungen anzunähern. FFT liefert komplexe Werte mit Imaginär- und Realteil. Sie geben jeweils zu einer Frequenz den Sinus- und Kosinus-Anteil an. Für die Visualisierung interessiert man sich aber für den Betrag des komplexen Wertes also sqrt(re² + im²). Was Du also siehst, sind die Frequenzanteile im aufgenommenen Signal. Tiefe Frequenzen entsprechen dabei Bässen, hohe Frequenzen hellen Tönen. Spielst Du bspw. auf einer Geige die Note A (Kammerton) , sollte im Spektrum dann vor allem 440 Hz dominieren.

Im Übrigen gibt es die Umkehrfunktion iFFT (inverse FFT), die aus den Frequenzanteilen wieder Funktionswerte erzeugt (halt durch Übereinanderlagerung der Schwingungen). Wenn Du Dein Signal erst in den Frequenzbereich transformierst (FFT), die Anteile manipulierst, und dann die geänderten Anteile zurück transformierst (iFFT) hast Du einen Equalizer. Ein Equalizer dämpft (schwächt ab) oder verstärkt (erhöht) die Frequenzanteile.

Falls Du wirklich verstehen möchtest, wie man ein Signal in verschiedene Frequenzen zerlegt, empfehle ich Dir Dich mit der diskreten Kosinustransformation (DCT) beschäftigen. Sie betrachtet nur die Kosinusanteile und ist reellwertig. Eine hübsche Beschreibung findest Du unter http://www.inf.fh-flensburg.de...ft/dkt.htm .

Ciao Olli
vertex.dreamfall.at | GitHub
 

timmeTheOnly

 BeitragMi, Nov 02, 2011 11:12
Antworten mit Zitat
Benutzer-Profile anzeigen
Danke Vertex!

Hab jetzt das Sample in einer Bank abgelegt, komm aber mit der DCT nicht weiter.
Ich hab sie als Funktion implementiert, allerdings weiß ich nicht so recht, was ich damit anfangen soll.

Das liegt zum einen an der Funktion selbst:

Soweit ich es verstanden habe liefert mir die DCT die Summe der Koeeffizienten für verschiedene Kosinusfunktionen, aber was bringt mir das?

Zum anderen weiß ich nicht was die Variable i "macht".

Gruß,
Tim

DCT:
BlitzMax: [AUSKLAPPEN] 

Function DoDCT:Double(Sample:Short[],i:Double)
	Local Sum:Double = 0, s:Double = 0, n:Int = 0, j:Int = 0
	
	n = Sample.Length
	s = Sqr(2)/Sqr(n)
	
	For j = 0 To n-1
		Sum:+ Sample[j] * Cos((2*j*i+i)*Pi/2*n)
	Next
	
	Sum:* s
	
	Return Sum
End Function

Vertex

 BeitragMi, Nov 02, 2011 12:55
Antworten mit Zitat
Benutzer-Profile anzeigen
Ich möchte es nicht beschwören, aber ich vermute mal Dein Code ist falsch.
- Zum einem nimmt Cos unter BMax den Winkel in Grad und nicht im Bogenmaß. Das heißt, wenn Du Dir den Ausdruck aus einer Formel hergeleitet hast , musst Du wissen, dass üblicherweise cos das Bogenmaß entgegen nimmt. Gleiches gilt, wenn Du die Funktion aus einem Java- oder C-Quelltext kopiert hast - auch dort wird der Winkel im Bogenmaß angegeben.
- Dann sehe ich noch ein Problem: s ist sicher der "Vorfaktor". Beachte:
Code: [AUSKLAPPEN]
s0 = 1/sqrt(n)
si = sqrt(2)/sqrt(n) für i > 0

Du nimmst stets i > 0 an.

Und nun zur Interpretation:
Also, Du hast ein Zeitfenster [y_0 ... y_n-1] aufgenommen.
Es wird angenommen, dass sei eine Funktion, die von 0 bis Pi läuft.
Wenn Du mit einer Samplefrequenz von sagen wir 4096 Hz auf nimmst, und Du das Zeitfenster [y_0 ... y_511] betrachtest, bekommst Du genau 4 Zeitfenster in der Sekunde geliefert. In den Zeitfenstern steckt allerdings eine volle Periode, also eine Funktion von 0 bis 2Pi. Das bedeutet, die Frequenzanteile nach der DCT auf [y_0 ... y_511] entsprechen den Frequenzen mal 4 durch 2. Wenn ich mich jetzt nicht irre bedeutet das:
Code: [AUSKLAPPEN]
i entspricht i*4/2 Hz


Die DCT liefert Dir die Anteile der Kosinusschwingungen zurück, um Dein aufgenommenes Fenster zu rekonstruieren. Anteil entspricht den Koeefizienten der Kosinusfunktion und damit den Amplituden der Kosinusschwingungen. Jetzt muss aber beachtet werden, dass die Amplitude mit einem Vorfaktor si skaliert wurde, den Du nach der DCT heraus dividieren musst. Je größer die Amplitude ist, desto stärker ist die Frequenz in Deinem Fenster enthalten.

Zusammengefasst:
i entspricht einer Frequenz
DCT gibt dafür den Anteil der Frequenz in Deinem aufgenommenen Fenster zurück.

Vllt. verwirre ich Dich jetzt noch, aber ich habe hier noch einen alten Java-Code gefunden. Er ist nach http://www.inf.fh-flensburg.de...ft/dkt.htm implementiert:
Code: [AUSKLAPPEN]
public class DCT {
   private final float[][] t; // Transformation-Matrix T
   private final int n; // Größe der Matrix
   private boolean inverted;

   public DCT(final int n) {
      t = new float[n][n];
      this.n = n;
      buildDCT();
   }

   private void buildDCT() {
      // Vorfaktor
      final float s0 = (float) (1.0 / Math.sqrt(n));
      final float si = (float) (Math.sqrt(2.0) / Math.sqrt(n));
      float s;

      for (int i = 0; i < n; i++) {
         for (int j = 0; j < n; j++) {
            if (i == 0)
               s = s0;
            else
               s = si;

            t[i][j] = (float) (s * Math.cos(i * (j + 0.5) * (Math.PI / n)));

         }
      }

      inverted = true;
   }

   private void transpose() {
      float temp;

      for (int i = 0; i < n; i++) {
         for (int j = 0; j < n; j++) {
            if (i > j) {
               temp = t[i][j];
               t[i][j] = t[j][i];
               t[j][i] = temp;
            }
         }
      }

      inverted = !inverted;
   }

   public void forwardDCT(final float[] y, final float[] a) {
      if (y == null || a == null)
         throw new NullPointerException("y and a must not be null");
      if (!(y.length == n && a.length == n))
         throw new IllegalArgumentException("y and a must have the length "
               + n);

      if (!inverted)
         transpose();

      // a = (T^-1).y
      for (int i = 0; i < n; i++) {
         for (int j = 0; j < n; j++) {
            a[i] += t[i][j] * y[j];
         }
      }
   }


   public void inverseDCT(final float[] y, final float[] a) {
      if (y == null || a == null)
         throw new NullPointerException("y and a must not be null");
      if (y.length == n && a.length == n)
         throw new IllegalArgumentException("y and a must have the length "
               + n);

      if (inverted)
         transpose();

      // y = T.a
      for (int i = 0; i < n; i++) {
         for (int j = 0; j < n; j++) {
            y[i] += t[i][j] * a[j];
         }
      }
   }
}


Ich habe jetzt leider keine Zeit, um Dir das mal an einem Beispiel zu zeigen. Vllt. versuchst Du Dir selber mal ein Fenster zu generieren, in dem drei unterschiedliche Frequenzen mit unterschiedlichen Anteilen vorkommen. Und dann versuchst Du die Frequenzanteile nach der DCT zu bestimmen.

Ciao Olli
vertex.dreamfall.at | GitHub
 

timmeTheOnly

 BeitragMi, Nov 02, 2011 13:05
Antworten mit Zitat
Benutzer-Profile anzeigen
Vielen Dank Vertex!

Das mit dem Bogenmaß hab ich komplett verpennt.

Ich schau mir den Code mal an und versuche ihn auf BlitzMax zu portieren und hoffentlich auch zu verstehen

Gruß,
Tim

Vertex

 BeitragMi, Nov 02, 2011 14:07
Antworten mit Zitat
Benutzer-Profile anzeigen
Habe noch was für Dich...

Es geht darum, ein Fenster zu generieren, womit man mal die DCT testen kann. Das Fenster besteht aus drei übereinander gelagerten Kosinus-Schwingungen mit 440 Hz, 230 Hz und 123 Hz. Dabei kann man die Samplerate und die Fenstergröße angeben. Heraus kommt das Fenster [y_0 ... y_n - 1] (n = WINDOW_SIZE).

Damit man das auch mal anschauen kann, wird eine Raw-Datei erzeugt. Die Raw-Datei kannst Du mit Audacityimportieren. Dazu muss die Projekt-Frequenz auf 4096 Hz (= SAMPLING_RATE) gestellt werden. Die Datei kann durch Datei -> Import > Rohdaten importiert werden, mit:
Codec: Signed 16 bit PCM
Byte-Reihenfolge. Little Endian
0-Byte-Offset: 0
Datenmenge: 100%
Samplefrequenz: 4096 Hz (= SAMPLING_RATE)

Code: [AUSKLAPPEN]
Const SAMPLE_RATE    = 4096 ; Hz
Const WINDOW_SIZE    = 512  ; Samples
Const AMPLITUDE_440# = 0.5  ; Anteil der Cos-Schwingung mit 440 Hz
Const AMPLITUDE_230# = 0.2  ; Anteil der Cos-Schwingung mit 230 Hz
Const AMPLITUDE_123# = 0.1  ; Anteil der Cos-Schwingung mit 123 Hz


Local file, i, t#
Local sample
Dim y#(WINDOW_SIZE - 1) ; Fenster mit Samples in [-1, +1]

; Format: Raw, Mono, Little Endian, Signed 16 bit in PCM
; Samplerate: SAMPLE_RATE
file = WriteFile("sample.raw")
For i = 0 To WINDOW_SIZE - 1
   t = (Float(i) / Float(SAMPLE_RATE)) * (2.0 * Pi)
   
   y(i) = AMPLITUDE_440 * CosRad(440.0 * t)
   y(i) = y(i) + AMPLITUDE_230 * CosRad(230.0 * t)
   y(i) = y(i) + AMPLITUDE_123 * CosRad(123.0 * t)
   
   
   If y(i) > 1.0 Then y(i) = 1.0
   If y(i) < -1.0 Then y(i) = -1.0
   
   sample = (y(i) * $7FFF)
   sample = sample And $7FFF ; bereits als 2er Komplement
   If y(i) < 0.0 Then sample = (sample Or $8000) ; Vorzeichenbit

   WriteShort(file, sample)
   
Next
CloseFile(file)

; Kosinus von a im Bogenmaß
Function CosRad#(a#)
   ;  a     a'
   ; --- = ----
   ; 2Pi   360°
   
   ; (a im Bogenmaß, a' im Gradmaß)
   Return Cos((a * 360.0) / (2.0*Pi))
End Function


Das ganze sieht dann so aus:
user posted image

Habe gleich mal die Frequenzanalyse ins Bild rein genommen (hier allerdings FFT). Du siehst die entsprechnenden Peaks bei 440 Hz, 230 Hz und 123 Hz.

Jetzt müsstest Du auf y(i = 0 ... i = WINDOW_SIZE - 1) mal die DCT anwenden, und ebenfalls die Anteile 0.5, 0.2, 0.1 heraus bekommen (wie gesagt, den Vorfaktor nicht vergessen!)

Ciao Olli
vertex.dreamfall.at | GitHub
 

timmeTheOnly

 BeitragMi, Nov 02, 2011 15:45
Antworten mit Zitat
Benutzer-Profile anzeigen
Also,

ich hab das ganze mal versucht und ich bekomme auch die drei Peaks.
Allerdings stimmen meine Werte auf den Peaks nicht.

Muss ich auch wieder durch den Faktor t dividieren?

Also:

Code: [AUSKLAPPEN]
Anteil = a[i] / (s*t)


Gruß,
Tim

Vertex

 BeitragMi, Nov 02, 2011 17:46
Antworten mit Zitat
Benutzer-Profile anzeigen
Nicht vergessen, Audacity gibt die Ampltidue in Dezibel (dB) an. Kann Dir allerdings nicht sagen, wie Du die dB berechnen kannst. Dezibel ist immer ein logarithmisches Verhältnis zweier Größen ("Bel") multipliziert mit 10 ("Dezi"). Was allerdings Audacity als Bezugsamplitude nimmt, weiß ich nicht.

Hier etwas Code für Dich:
BlitzBasic: [AUSKLAPPEN] 
Const SAMPLE_RATE    = 4096 ; Hz
Const WINDOW_SIZE    = 512  ; Samples
Const AMPLITUDE_440# = 0.5  ; Anteil der Cos-Schwingung mit 440 Hz
Const AMPLITUDE_230# = 0.2  ; Anteil der Cos-Schwingung mit 230 Hz
Const AMPLITUDE_123# = 0.1  ; Anteil der Cos-Schwingung mit 123 Hz


Dim y#(WINDOW_SIZE - 1) ; Fenster mit Samples in [-1, +1]
Dim a#(WINDOW_SIZE - 1) ; Stützstellen
Dim t#(WINDOW_SIZE - 1, WINDOW_SIZE - 1) ; Transformationsmatrix T

GenerateWindow()
SaveRaw("sample.raw")
InitDCT()
ForwardDCT()
DebugSpectrum()


; Erzeugt Überlagerung von drei Kosinus-Wellen
; mit 440 Hz, 230 Hz und 123 Hz.
; Ausgabe ist das Fenster y(i = 0 ... i = WINDOW_SIZE - 1)
; Jedes Sample liegt in [-1; +1]
Function GenerateWindow()
	Local i, sample#, t#
	
	For i = 0 To WINDOW_SIZE - 1
		t = (Float(i) / Float(SAMPLE_RATE)) * (2.0 * Pi)
		
		sample = AMPLITUDE_440 * CosRad(440.0 * t)
		sample = sample + AMPLITUDE_230 * CosRad(230.0 * t)
		sample = sample + AMPLITUDE_123 * CosRad(123.0 * t)
		
		; Begrenze Sample auf [-1; +1]
		If sample > 1.0 Then sample = 1.0
		If sample < -1.0 Then sample = -1.0
		
		y(i) = sample
	Next
End Function


; DCT-Transformations-Matrix T initialisieren
; Matrix ist WINDOW_SIZE x WINDOW_SIZE groß
Function InitDCT()
	Local s0#, si#, s# ; Vorfaktor s
	Local i, j
	
	s0 = 1.0 / Sqr(WINDOW_SIZE) ; für i = 0
	si = Sqr(2.0) / Sqr(WINDOW_SIZE) ; für i > 0
	
	For i = 0 To WINDOW_SIZE - 1
		For j = 0 To WINDOW_SIZE  - 1
			If i = 0 Then
				s = s0
			Else
				s = si
			EndIf
			
			t(i, j) = s * CosRad(i * (j + 0.5) * (Pi / WINDOW_SIZE)) 
		Next
	Next
End Function

; Matrix T ausgeben
Function DebugT()
	Local i, j, li$
	
	For i = 0 To WINDOW_SIZE - 1 ; i = Row
		li = ""
		For j = 0 To WINDOW_SIZE  - 1 ; j = Col
			li = li + t(i, j) + "	"
		Next
		li = Left(li, Len(li) - 1) ; Letzter Tab wird entfernt
		DebugLog "| " + li + " |"
	Next
End Function

; Transformiert die Stützstellen-Darstellung y()
; in die Koeeffizientendarstellung a()
; Formal a^T = y^T.T^(-1) mit T^(-1) = T^T
; (^T heißt transponiert)
Function ForwardDCT()
	Local i, j
	
	For i = 0 To WINDOW_SIZE - 1
		a(i) = 0
		For j = 0 To WINDOW_SIZE - 1
			; merke t(i, j) ist die Transponierte von t(j, i) 
			a(i) = a(i) + t(i, j) * y(j)
		Next
	Next
End Function

; Gibt die Amplituden für jede Frequenz aus
Function DebugSpectrum()
	Local s0#, si#, s# ; Vorfaktor s
	Local i
	Local frequency#, amplitude#
	
	s0 = 1.0 / Sqr(WINDOW_SIZE) ; für i = 0
	si = Sqr(2.0) / Sqr(WINDOW_SIZE) ; für i > 0
	
	For i = 0 To WINDOW_SIZE - 1
		If i = 0 Then
			s = s0
		Else
			s = si
		EndIf
		
		frequency = i * (Float(SAMPLE_RATE) / Float(WINDOW_SIZE)) / 2
		amplitude = a(i) * s
		
		; Nur interessante Frequenzen ausgeben
		If Abs(amplitude) > 0.09 Then DebugLog "Amplitude für " + Int(frequency) + " Hz : " + amplitude
	Next
End Function

; Format: Raw, Mono, Little Endian, Signed 16 bit in PCM
; Samplerate: SAMPLE_RATE
Function SaveRaw(file$)
	Local stream
	Local i, sample
	
	stream = WriteFile(file$)
	For i = 0 To WINDOW_SIZE - 1		
		sample = y(i) * $7FFF
		sample = sample And $7FFF ; bereits als 2er Komplement
		If y(i) < 0.0 Then sample = (sample Or $8000) ; Vorzeichenbit
	
		WriteShort(stream, sample)
	Next
	CloseFile(stream)
End Function


; Kosinus von a im Bogenmaß
Function CosRad#(a#)
	;  a     a'
	; --- = ----
	; 2Pi   360°
	
	; (a im Bogenmaß, a' im Gradmaß)
	Return Cos((a * 360.0) / (2.0*Pi))
End Function


Das Fenster wird wie zuvor berechnet, setzt sich also aus 3 Kosinuswellen zusammen.

Dann wird die T-Matrix aufgebaut mit InitDCT. Das Fenster y() wird nun mit der Inversen von T multipliziert, um von der Stützstellendarstellung in die Koeeffizientendarstellung zu kommen. Die Inverse von T ist netterweise die Transponierte von T, weswegen nur eine normale Vektor-Matrix-Multiplikation mit vertauschten Zeilen und Spalten durchgeführt werden muss.

In der Debuglog siehst Du zu den Frequenzen die Amplituden. Die stimmen in etwa überein:
Code: [AUSKLAPPEN]
124 Hz: 0.0915205
228 Hz: 0.151581
440 Hz: 0.472312

vs. die Idealen Amplituden
Code: [AUSKLAPPEN]
124 Hz: 0.1 (eigtl. für 123 Hz!)
228 Hz: 0.2 (eigtl. für 230 Hz)
440 Hz: 0.5


Was Du machen kannst, ist mal alle Frequenzen so auszugeben:
WriteLine(file, Replace(Int(frequency) + " " + amplitude, ".", ","))
Sie in Excel zu laden, und ein X-Y-Plot zu machen. Dann siehst Du das Frequenzspektrum des Fensters, ähnlich wie in Audacity:
user posted image

Ciao Olli
vertex.dreamfall.at | GitHub

Vertex

 BeitragMi, Nov 02, 2011 21:50
Antworten mit Zitat
Benutzer-Profile anzeigen
Habe das ganze mal noch der bass.dll kombiniert. Es wird eine MP3 abgespielt, und davon werden regelmäßig die Samples abgefragt (BASS_ChannelGetData). Dieses Fenster wird DCT transformiert und das Frequenzspektrum in einem Kreis dargestellt (DrawSpectrumRadial).

BlitzBasic: [AUSKLAPPEN] 
Const SAMPLE_RATE    = 4096 ; Hz
Const WINDOW_SIZE    = 512  ; Samples
Const AMPLITUDE_440# = 0.5  ; Anteil der Cos-Schwingung mit 440 Hz
Const AMPLITUDE_230# = 0.2  ; Anteil der Cos-Schwingung mit 230 Hz
Const AMPLITUDE_123# = 0.1  ; Anteil der Cos-Schwingung mit 123 Hz

Include "bass.bb"

Dim y#(WINDOW_SIZE - 1) ; Fenster mit Samples in [-1, +1]
Dim a#(WINDOW_SIZE - 1) ; Stützstellen
Dim t#(WINDOW_SIZE - 1, WINDOW_SIZE - 1) ; Transformationsmatrix T

;GenerateWindow()
;SaveRaw("sample.raw")


InitDCT()

If Not  BASS_Init(-1, SAMPLE_RATE, 0, 0, BASS_NULL) Then RuntimeError "Bass not ready"
BASS_CheckVersion()

Global music = BASS_StreamCreateFile(0, "tiesto.mp3", 0, 0, BASS_MUSIC_FLOAT + BASS_SAMPLE_LOOP)
Global samples = CreateBank(WINDOW_SIZE * 4)

Graphics(640, 480, 0, 2)
SetBuffer(BackBuffer())

BASS_ChannelPlay(music, BASS_TRUE)

Local i
While Not KeyDown(1)
	If BASS_ChannelIsActive(music) Then
		BASS_ChannelGetData(music, samples, (WINDOW_SIZE * 4) Or BASS_DATA_FLOAT)
		For i = 0 To WINDOW_SIZE - 1
			y(i) = PeekFloat(samples, i * 4)
		Next
		ForwardDCT()
	EndIf

	Cls()
	DrawSpectrumRadial()
	DrawSpectrum()
	Flip()
Wend

BASS_ChannelStop(music)
BASS_Free()

End


; Erzeugt Überlagerung von drei Kosinus-Wellen
; mit 440 Hz, 230 Hz und 123 Hz.
; Ausgabe ist das Fenster y(i = 0 ... i = WINDOW_SIZE - 1)
; Jedes Sample liegt in [-1; +1]
Function GenerateWindow()
	Local i, sample#, t#
	
	For i = 0 To WINDOW_SIZE - 1
		t = (Float(i) / Float(SAMPLE_RATE)) * (2.0 * Pi)
		
		sample = AMPLITUDE_440 * CosRad(440.0 * t)
		sample = sample + AMPLITUDE_230 * CosRad(230.0 * t)
		sample = sample + AMPLITUDE_123 * CosRad(123.0 * t)
		
		; Begrenze Sample auf [-1; +1]
		If sample > 1.0 Then sample = 1.0
		If sample < -1.0 Then sample = -1.0
		
		y(i) = sample
	Next
End Function


; DCT-Transformations-Matrix T initialisieren
; Matrix ist WINDOW_SIZE x WINDOW_SIZE groß
Function InitDCT()
	Local s0#, si#, s# ; Vorfaktor s
	Local i, j
	
	s0 = 1.0 / Sqr(WINDOW_SIZE) ; für i = 0
	si = Sqr(2.0) / Sqr(WINDOW_SIZE) ; für i > 0
	
	For i = 0 To WINDOW_SIZE - 1
		For j = 0 To WINDOW_SIZE  - 1
			If i = 0 Then
				s = s0
			Else
				s = si
			EndIf
			
			t(i, j) = s * CosRad(i * (j + 0.5) * (Pi / WINDOW_SIZE)) 
		Next
	Next
End Function

; Matrix T ausgeben
Function DebugT()
	Local i, j, li$
	
	For i = 0 To WINDOW_SIZE - 1 ; i = Row
		li = ""
		For j = 0 To WINDOW_SIZE  - 1 ; j = Col
			li = li + t(i, j) + "	"
		Next
		li = Left(li, Len(li) - 1) ; Letzter Tab wird entfernt
		DebugLog "| " + li + " |"
	Next
End Function

; Transformiert die Stützstellen-Darstellung y()
; in die Koeeffizientendarstellung a()
; Formal a^T = y^T.T^(-1) mit T^(-1) = T^T
; (^T heißt transponiert)
Function ForwardDCT()
	Local i, j
	
	For i = 0 To WINDOW_SIZE - 1
		a(i) = 0
		For j = 0 To WINDOW_SIZE - 1
			; merke t(i, j) ist die Transponierte von t(j, i) 
			a(i) = a(i) + t(i, j) * y(j)
		Next
	Next
End Function

Function DrawSpectrum()
	Local i
	Local x, y
	
	Local s0#, si#, s# ; Vorfaktor s
	s0 = 1.0 / Sqr(WINDOW_SIZE) ; für i = 0
	si = Sqr(2.0) / Sqr(WINDOW_SIZE) ; für i > 0
	
	For i = 0 To WINDOW_SIZE - 1
		If i = 0 Then
			s = s0
		Else
			s = si
		EndIf
	
		x = i + 64
		y = -Abs(a(i) * s) * 200 + 240
		
		Line(x, 240, x, y)
	Next
End Function

; Zeichnet die Punkte eines Kreises
; Der Radius jedes Punktes ist abhängig von der Amplitude
Function DrawSpectrumRadial()
	Local alpha#, t
	Local x, y, r#
	
	Local s0#, si#, s# ; Vorfaktor s
	s0 = 1.0 / Sqr(WINDOW_SIZE) ; für i = 0
	si = Sqr(2.0) / Sqr(WINDOW_SIZE) ; für i > 0
	
	LockBuffer()
	For alpha = 0.0 To 2.0 * Pi Step 0.04
		; alpha = 0 ... 2*PI wird abgebildet auf t = 0 ... WINDOW_SIZE / 8
		; also nur das erste Achtel des Spektrums wird ausgegeben
		t = Float(alpha / (2.0 * Pi)) * (WINDOW_SIZE / 8)
		If t < 0 Then t = 0
		If t >= WINDOW_SIZE Then t = WINDOW_SIZE - 1
		
		If t = 0 Then
			s = s0
		Else
			s = si
		EndIf
		
		r = 100 + 100 * Abs(a(t) * s) ; Radius = 100 px + Amplitude für Frequenz an Stelle t
		x = r * CosRad(alpha) + 320 ; Mittelpunkt ist (320, 240)
		y = r * SinRad(alpha) + 240
		
		; Cyan = RGB(0, 255, 255)
		WritePixelFast(x, y, $0000FFFF)
	Next
	UnlockBuffer()
End Function

; Gibt die Amplituden für jede Frequenz aus
Function DebugSpectrum()
	Local s0#, si#, s# ; Vorfaktor s
	Local i
	Local frequency#, amplitude#
	
	s0 = 1.0 / Sqr(WINDOW_SIZE) ; für i = 0
	si = Sqr(2.0) / Sqr(WINDOW_SIZE) ; für i > 0
	

	For i = 0 To WINDOW_SIZE - 1
		If i = 0 Then
			s = s0
		Else
			s = si
		EndIf
		
		frequency = i * (Float(SAMPLE_RATE) / Float(WINDOW_SIZE)) / 2
		amplitude = a(i) * s
		
		; Nur interessante Frequenzen ausgeben
		If Abs(amplitude) > 0.09 Then DebugLog "Amplitude für " + Int(frequency) + " Hz : " + amplitude
	Next
End Function

; Format: Raw, Mono, Little Endian, Signed 16 bit in PCM
; Samplerate: SAMPLE_RATE
Function SaveRaw(file$)
	Local stream
	Local i, sample
	
	stream = WriteFile(file$)
	For i = 0 To WINDOW_SIZE - 1		
		sample = y(i) * $7FFF
		sample = sample And $7FFF ; bereits als 2er Komplement
		If y(i) < 0.0 Then sample = (sample Or $8000) ; Vorzeichenbit
	
		WriteShort(stream, sample)
	Next
	CloseFile(stream)
End Function


; Kosinus von a im Bogenmaß
Function CosRad#(a#)
	;  a     a'
	; --- = ----
	; 2Pi   360°
	
	; (a im Bogenmaß, a' im Gradmaß)
	Return Cos((a * 360.0) / (2.0*Pi))
End Function

; Sinus von a im Bogenmaß
Function SinRad#(a#)
	Return Sin((a * 360.0) / (2.0*Pi))
End Function


Das sieht dann so aus:
user posted image

Die BASS.dll findest Du unter: http://www.blitzbasic.com/tool...p?tool=207

Da BASS.dll praktischerweise eine FFT-Analyse gleich mitliefert (statt BASS_DATA_FLOAT einfach BASS_DATA_FFT512 an BASS_GetChannelData übergeben), konnte ich die handgeschriebene DCT und FFT vergleichen. Siehe da, die tiefen Frequenzen sind gleich, die hohen nicht. Das liegt aber sicher daran, dass keine Fenster-Funktion, wie das Hamming-Fenster, vorher auf die Samples angewendet werden. Es wird bei der DCT angenommen, dass das Signal periodisch ist. Deswegen gibt es Sprünge am Ende des Fensters zum Anfang des Fensters, was sich durch hohe Frequenzen auszeichnet. Siehe Fensterfunktion

Ciao Olli

Edit: Meine Vermutung wird bestätigt, wenn man in der Doku von BASS nachliest:
Zitat:
[...] A Hanning window is applied to the sample data to reduce leakage, unless the BASS_DATA_FFT_NOWINDOW flag is used.
vertex.dreamfall.at | GitHub
 

timmeTheOnly

 BeitragDo, Nov 03, 2011 7:16
Antworten mit Zitat
Benutzer-Profile anzeigen
Ok vielen Dank!

Ich glaub jetzt bin ich soweit, dass ich das umgesetzt bekommen.
Ich melde mich, wenn es geklappt hat Wink

Gruß,
Tim

Neue Antwort erstellen


Übersicht BlitzMax, BlitzMax NG Allgemein

Nach Oben

Gehe zu:

Powered by phpBB © 2001 - 2006, phpBB Group