Forum: Digitale Signalverarbeitung / DSP / Machine Learning dft in java. Frage wie zurück in Zeitdomäne rechnen?


von jonas biensack (Gast)


Lesenswert?

Hallo,

aus Spass an der Freude wollte ich mich mal mit was anderem 
beschäftigen, heute mal dft. Ich habe dazu eine kleine java Klasse 
gebaut, da ich gerade erfolgreich eine asio Schnittstelle für Java zum 
Laufen gebracht hatte.
Hier mal der Code:
1
package com.synthbot.jasiohost;
2
3
public class myDFT {
4
  private float sampleRate;
5
  
6
  private static float sinTime[];
7
  private static float cosTime[];
8
  private int N;
9
  private float freq;
10
  
11
  public myDFT(float sampleRa,float freq, int N)
12
  {
13
    this.sampleRate = 44100;
14
    this.N = 512;
15
    this.freq = 440;
16
    sinTime = new float[N];
17
    cosTime = new float[N];
18
    sinTime = createSinTimeSeries(this.N, this.freq);
19
    cosTime = createCosTimeSeries(this.N, this.freq);
20
  }
21
  public float findFreq (float[] x, int minfreq, int maxfreq, int stepfreq)
22
  {
23
    float real = 0;
24
    float imag = 0;
25
    float toreturn=0;
26
    float powerAtFreq = 0;
27
    float MaxpowerAtFreq = 0;
28
    float[] samples = new float[this.N];
29
    for (int a = minfreq; a < maxfreq; a += stepfreq)
30
    {
31
      sinTime = createSinTimeSeries(this.N, a);
32
      cosTime = createCosTimeSeries(this.N, a);
33
      for(int bb = 0; bb < this.N; bb++)
34
      {
35
        real += x[bb] * sinTime[bb];
36
        imag += x[bb] * cosTime[bb];
37
        
38
      }
39
      powerAtFreq = (real * real) + (imag * imag);
40
      
41
      if (powerAtFreq > MaxpowerAtFreq)
42
      {
43
        MaxpowerAtFreq = powerAtFreq;
44
        toreturn = a;
45
      }
46
    }
47
    return toreturn;
48
    
49
  }
50
  public float[] inverseDFt(float[] real, float[] imag)
51
  {
52
    float[] samples = new float[this.N];
53
    
54
    for (int a = 0; a < this.N; a++)
55
    {
56
      
57
    }
58
    return samples;
59
  }
60
  public float calculateMyDFT(float[] x)
61
  {
62
    float real = 0;
63
    float imag = 0;
64
    
65
    float powerAtFreq = 0;
66
    
67
    for(int a = 0; a < this.N; a++)
68
    {
69
      real += x[a] * sinTime[a];
70
      imag += x[a] * cosTime[a];
71
    }
72
    powerAtFreq = (real * real) + (imag * imag);
73
    return powerAtFreq;
74
    //System.out.println("Power at freq " +  freq + " is " + powerAtFreq);
75
  }
76
  
77
  public float[] createSinTimeSeries(int N, float freq)
78
  {
79
    float output[] = new float[N];
80
    for (int i = 0; i < N; i++) {
81
          output[i] = (float) (Math.sin(2 * Math.PI * i * freq / this.sampleRate));
82
        }
83
    return output;
84
    
85
  }
86
  public float[] createCosTimeSeries(int N, float freq)
87
  {
88
    float output[] = new float[N];
89
    for (int i = 0; i < N; i++) {
90
          output[i] = (float) (Math.cos(2 * Math.PI * i * freq / this.sampleRate));
91
        }
92
    return output;
93
  }
94
}

Funktioniert auch soweit. In die eine Richtung also von Timedomaine nach 
Frequenzdomaine, aber wie zurück. Kann man mal jemand auf die Sprünge 
helfen, wie bekomme ich aus einem Array von real und imaginär Anteil 
wieder ein Sample zurück? Hoffe ihr was ich meine ;)

Gruß Jonas

von Antworter (Gast)


Lesenswert?

Ohne deinen Code jetzt größer angeschaut zu haben:

Vom Frequenzbereich in den Zeitbereich kommst du mit der inversen DFT:

http://de.wikipedia.org/wiki/IDFT#DFT_und_iDFT_f.C3.BCr_einen_komplexen_Vektor


Ich hoffe das hilft.

von Christopher B. (chrimbo) Benutzerseite


Lesenswert?

>public myDFT(float sampleRa,float freq, int N)
>  {
>    this.sampleRate = 44100;
>    this.N = 512;
>    this.freq = 440;
>    ...
>  }

welchen Sinn machen die Übergabeparameter?

von jonas biensack (Gast)


Lesenswert?

>welchen Sinn machen die Übergabeparameter?

Gar keinen, wollte nur ein paar Werte hier vorgeben. Sonst hätte ich den 
Aufruf noch mit posten müssen.

>http://de.wikipedia.org/wiki/IDFT#DFT_und_iDFT_f.C...
>Ich hoffe das hilft.

Das weiss ich schon, aber wie implemtiere ich das? Und noch schlimmer, 
ich will ja wissen, was ich da programmiere. Einfach die java.fft libs 
einbinden hat keinen Lerneffekt für mich.

gruß jonas

von jonas biensack (Gast)


Lesenswert?

Nochmal zusammenfassen was ich bis dato gemacht habe und warum.
Laut dieser Seite:
http://www.developer.com/java/other/article.php/3374611
soll ich folgendes tun um eine dst auszuführen:

Gegeben sei ein Array (x) mit den zu analysierenden Signal. Die Analyse 
selbst besteht dann darin, für eine zu untersuchende Frequenz (freq) ein 
Sinus-Zeitserie (frei übersetzt "sine time series") sowie eine 
Cosinus-Zeitserie zu erstellen. Alle Elemente aus der Sinus-Zeitserie 
werden nun iterativ mit den Elementen in x multipliziert und das 
Ergebnis kontinueirlich aufaddiert (real). Das selbe für die 
Cosinus-Zeitserie, diese ergibt dann den imaginären Anteil dieser 
komplexen Zahl. Addiere ich nun zum Schluss das Quadrat aus dem reelen 
Teil (real) plus dem Quadrat aus dem imaginären Teil (imag) bekomme ich 
einen Wert "power".
Wäre das der TRUE-RMS wert??
Also die Fläche unter dem Sinusteil + die Fläche über dem Cosinusteil?

Nun habe ich also zum Beispiel einen Sinus mit 440 hz als zu 
analysierendes Signal. Ich weiss also welches Frequenzband ich mir näher 
anschauen sollte (zur Vereinfachung) Also mach ich zum Beispiel zehn 
Analysen, beginnend bei 435 hz bis zu 445 hz in ein herz Schritten. Als 
Ergbenis habe ich 10 komplexe Zahlen. Wie kann ich nun mit diesen meinen 
440 hz Sinuston vom Anfang rekunstruieren?

gruß jonas

von chris (Gast)


Lesenswert?

Du bildest als erstes die Beträge aus den sin- und cos-Koeffizienten.
Dann suchst Du den maximalen Betrag und weist, welche Frequenz vorhanden 
ist.
Die ursprüngliche Schwingung erhälst Du einfach so

x(n)=a*cos(2*pi*fsig/fab*n)+b*sin(2*pi*fsig/fab*n)

Stichwort: Fourierreihe
http://de.wikipedia.org/wiki/Fourierreihe

von chris (Gast)


Lesenswert?

>Du bildest als erstes die Beträge aus den sin- und cos-Koeffizienten.
Das war etwas missverständlich. Mit Betrag meine ich den geometrischen 
Betrag aus den Koeffizienten:

Betrag=sqrt(a²+b²)

von jonas biensack (Gast)


Lesenswert?

Danke erstmal für die Antwort.

>Du bildest als erstes die Beträge aus den sin- und cos-Koeffizienten.
>Dann suchst Du den maximalen Betrag und weist, welche Frequenz vorhanden >ist.

Ist das noch Teil der dft oder gehört das zur inversiven dft, also ist 
die Antwort auf meine Frage?

>Die ursprüngliche Schwingung erhälst Du einfach so

>x(n)=a*cos(2*pi*fsig/fab*n)+b*sin(2*pi*fsig/fab*n)

Was ist in dem Fall:

a?
b?
fsig?
fab = betrag von cosinus und sinus?

Sorry aber jetzt bin ich noch verwirrter...
Noch ein paar Gedanken von mir. Angenommen ich habe im Signal nur eine 
Sinusschwingung (wie oben 440 hz). Mit dem von mir verwendeten 
Alogrithmus findfreq kann ich die freq mit der höchsten Amplitude aus 
dem Signal bestimmen. Könnte ich nicht zur Rekonstruktion wieder die 
gleichen Sinus und Cosinus-Zeitserien benutzen und die gewichtet (nach 
der ermittelten Amplitude) miteinander mischen?

gruß Jonas

von jonas biensack (Gast)


Lesenswert?

>>Du bildest als erstes die Beträge aus den sin- und cos-Koeffizienten.
>>Das war etwas missverständlich. Mit Betrag meine ich den geometrischen
>>Betrag aus den Koeffizienten:

>Betrag=sqrt(a²+b²)

Mach ich nicht genau das hier, ok die Wurzel hab ich vergessen.
1
 public float findFreq (float[] x, int minfreq, int maxfreq, int stepfreq)
2
  {
3
    float real = 0;
4
    float imag = 0;
5
    float toreturn=0;
6
    float powerAtFreq = 0;
7
    float MaxpowerAtFreq = 0;
8
    float[] samples = new float[this.N];
9
    for (int a = minfreq; a < maxfreq; a += stepfreq)
10
    {
11
      sinTime = createSinTimeSeries(this.N, a);
12
      cosTime = createCosTimeSeries(this.N, a);
13
      for(int bb = 0; bb < this.N; bb++)
14
      {
15
        real += x[bb] * sinTime[bb];
16
        imag += x[bb] * cosTime[bb];
17
        
18
      }
19
      //vorher
20
      //powerAtFreq = (real * real) + (imag * imag);
21
      //nachher nach deinem Hinweis:
22
      //powerAtFreq = (float)Math.sqrt((real * real) + (imag * imag));
23
24
      if (powerAtFreq > MaxpowerAtFreq)
25
      {
26
        MaxpowerAtFreq = powerAtFreq;
27
        toreturn = a;
28
      }
29
    }
30
    return toreturn;
31
    
32
  }

Gruß JOnas

von chris (Gast)


Lesenswert?

>Mach ich nicht genau das hier, ok die Wurzel hab ich vergessen.
Genau, mit der Wurzel rechnest Du die Amplitude der Schwingung aus.
Da die Schwingung zeitlich verschoben sein kann ( Phasenverschiebung ), 
erhälst Du aus der DFT einen Cosinus- und Sinusanteil, die in Deinem 
Programm als Real und Imaginärteil benannt sind

        real += x[bb] * sinTime[bb];
        imag += x[bb] * cosTime[bb];

Diese Bezeichnungen kommen allerdings aus der Darstellung einer 
Schwingung als komplexer Drehzeiger

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

Bei der in Deinem Programm Sinus- und Cosinus vertauscht ist. Was nicht 
schlimm ist, wenn man das bei der Wiederherstellung der Schwingung 
beibehält.

Du müsstest Dir den Imaginär- und Realteil bei der maximalen Amplituden 
merken, dann kannst Du die Schwingung mit

sig(n)= real * sin ... + imag*cos ...

wiederherstellen.

Siehe Überlagerung von Schwingungen gleicher Frequenz:
http://schulphysikwiki.de/index.php/%C3%83%C5%93berlagerung_von_harmonischen_Schwingungen

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.