diskrete Fourier Transformation (DFT) - exakte Frequenz entfernen

Nicht wirklich.

Hirn einschalten. Höchste Frequenz ist 22050 Hz.

Korrekt.

Du wirst um ein digitales Notch-Filter nicht herumkommen. Das kann man schon über FFT implementieren, aber dann über einen overlap-shift bzw. overlap-add-Algorithmus. Das Filterdesign hat dann nichts mit FFT zu tun, sondern macht einfach ein FIR-Filter der gewünschten Länge, das dann mit Hilfe der FFT implementiert wird. 50Hz ist allerdings im Verhältnis zu 44100Hz eine verdammt tiefe Frequenz (Frequenzverhältnis 882), da ist ein FIR-Filter eine miserabel ineffiziente Lösung, wenn man auch nur geringe Genauigkeit erzielen will: die Filterantwort müßte wahnsinnig lang werden.

Viel besser ist hier ein IIR-Filter. Das kann man mit wenigen Elementen designen, und wenn man es einmal vorwärts und einmal rückwärts über die Daten laufen läßt, dann führt es nicht einmal zu Phasenfehlern.

Typische Transferfunktion wäre

H(z) = (1+a z^{-1} + z^{-2}) / (1 + a r z^{-1} + r² z^{-2})

Mit a=-2 cos w (wobei w = 2 pi * 50Hz / (44100 Hz))

Und r die Güte ist (zwischen 0 und 1, wobei 1 die höchste Güte ist, bei der das ausgeschnittene Frequenzband unendlich schmal ist, das Filter also nichts tut. Knapp darunter ist die Trennschärfe am größten, aber auch die Einschwingvorgänge am längsten).

Implementation:

a = ...; r = ...; y1 = y2 = x1 = x2 = 0.0; for (i=0; i

"Markus" schrieb im Newsbeitrag news:452211a2$0$17405$ snipped-for-privacy@newsspool2.arcor-online.net...

Hallo Markus,

Dein Ansatz ist schon mal schlecht, wenn du daran denkst das Signal später wieder in den Zeitbereich zurückzutransformieren. Da du von Audio sprichst, gehe ich davon aus, dass das der Fall ist. Dauert dein Signal nur 0.1sec? Wohl eher nicht. Mit den DFTs läufst du schon mal gleich in Probleme bei der Rücktransformation. An den Übergangsstellen zwichen den Zeitabschnitten werden Artifakte auftreten.

Die Lösung: Programmiere ein Notch-Filter für 50Hz. Da läßt du dein Signal einfach komplett durchlaufen und fertig ist die Filterung.

Gruß Helmut

"David Kastrup"

s.o.

Kannst du das bitte nochmal mit Einheiten ergänzen? Ich verstehe so nicht, was bei dir die 10.7666 sind.

Wenn man die Einheiten korrekt nimmt teilst du nur die Anzahl der Samples durch eine einheitenlose Zahl. Also soll die Zahl wohl eine Wellenlänge repräsentieren, wenn man sie mit 2 multipliziert.

Ich habe gerade nochmal nachgeschaut. Mein Algorithmus gibt mir bei einem 4096 langen Sampleblock

Aber welche Werte der 2048 Zahlen repräsentieren jetzt welche Frequenzen? 0 = 0Hz und 2048 = 22050Hz ist klar. Der Rest ist mir unklar.

Danke, werde mich mal drüber schlaumachen.

MfG,

Markus

Scherzkeks. Du hast rausredigiert, daß Du dem Koeefizienten 4095 eine Frequenz von 44100Hz zugewiesen hast, die weit über Deiner höchsten deklarierten Frequenz von 22050Hz lag.

Herrschaftszeiten. 44100Samples/s / 4096 Samples = 10.7666015625 Hz.

Das ist ja nun wirklich nicht so schwer, das zu ergänzen.

Eine Wellenlänge?!? Wo kommen denn hier bitteschön Längen ins Spiel?

Falsch. Wenn man einen Algorithmus verwendet, der die komplexen Redundanzen wegläßt, bekommt man 2050 Realwerte und 2046 Imaginärwerte, weil, wie bereits erwähnt, die Werte bei 0 und 2048 keine Phaseninformation haben.

Nö.

Dann solltest Du vielleicht meine Antwort von vorhin nochmal lesen. 0 und 2048 stehen alleine. Ansonsten bildet 1 und 4095 ein Quadraturpaar für 1*(44100/4096)Hz, 2 und 4094 für 2*(44100/4096)Hz,

Du mußt Dir aber auch klar machen, daß die Werte einer FFT nur dann sinnvoll interpretierbar sind, wenn es sich um ein 4096-periodisches Signal handelt, nicht aber um einen Ausschnitt aus Signalen mit anderer Periodenlänge.

Probier einfach den Code aus.

Scherzkeks. Du hast rausredigiert, daß Du dem Koeefizienten 4095 eine Frequenz von 44100Hz zugewiesen hast, die weit über Deiner höchsten deklarierten Frequenz von 22050Hz lag.

Herrschaftszeiten. 44100Samples/s / 4096 Samples = 10.7666015625 Hz.

Das ist ja nun wirklich nicht so schwer, das zu ergänzen.

Eine Wellenlänge?!? Wo kommen denn hier bitteschön Längen ins Spiel?

Falsch. Wenn man einen Algorithmus verwendet, der die komplexen Redundanzen wegläßt, bekommt man 2049 Realwerte und 2047 Imaginärwerte, weil, wie bereits erwähnt, die Werte bei 0 und 2048 keine Phaseninformation haben.

Nö.

Dann solltest Du vielleicht meine Antwort von vorhin nochmal lesen. 0 und 2048 stehen alleine. Ansonsten bildet 1 und 4095 ein Quadraturpaar für 1*(44100/4096)Hz, 2 und 4094 für 2*(44100/4096)Hz,

Du mußt Dir aber auch klar machen, daß die Werte einer FFT nur dann sinnvoll interpretierbar sind, wenn es sich um ein 4096-periodisches Signal handelt, nicht aber um einen Ausschnitt aus Signalen mit anderer Periodenlänge.

Probier einfach den Code aus.

Markus schrieb:

Hallo,

Netzbrumm in Audiodateien muss nicht immer unbedingt 50 Hz haben, Durch die Gleichrichtung im Netzteil entsteht ein Brummanteil mit einer Grundfrequenz von 100 Hz und Vielfachen davon.

Dadurch ist Filterung von Netzbrummen in der Audiotechnik in der Praxis oft nicht trivial. Brummen zu vermeiden (z.B. durch die Schaltung, Pegelwahl usw.) ist zumeist besser.

Bernd Mayer

"David Kastrup"

Ja. Das deckt sich auch mit den Ergebnissen aus meinem DFT-Vektor.

Dort habe ich als Input (als Beispiel mit 8 Samples)

wobei x die Samples sind.

Transformiert ergibt das dann:

r i r i r i r i r i r i r i r i

wobei r Realteil und i Imaginärteil ist.

Was mir allerdings nach wir vor schleierhaft ist, ist die Position der Maxima und des Offsets 0.

Dass es einen Realwert mehr als Imaginärwert gibt klingt erstmal nachvollziehbar, weil der eine Realwert mehr dann der Offset ist. Aber warum ein Imaginärwert weniger?

Was soll das denn sein? Realwert von 0, Imaginärwert von x? Also rein imaginäre Zeitsignale? Das gibt herzlich wenig Sinn.

Das hört sich vollkommen falsch an. Entweder hast Du eine vollständige komplexe DFT, dann ist der Input falsch beschrieben oder gewählt, oder Du hast eine real->komplex DFT, dann hat die Ausgabe eine andere Struktur, weil nur 7 Koeffizienten redundant sind, aber dafür 2 rein reell.

Vielleicht solltest Du mal angeben, was Deine DFT-Routine als Eingabewerte verlangt und als Ausgabewerte liefert.

Sonst kann man hier beliebigen Unsinn rumraten.

Bei dem Zeug, was Du da so zusammengeschrieben hast, wird _jedem_ diese Position herzlich unklar sein.

Wird das Realsatire? Du zitierst doch den Grund:

"David Kastrup"

Nun gut. Da offensichtlich die Definition des DFT-Inputs nicht ganz klar ist, stelle ich an dieser Stelle mal mein Miniprogramm rein, was ich mir gestern gebastelt habe, um alles besser zu verstehen:

Dort kann man 8 Inputwerte einstellen, die dann quasi eine Folge von Samples darstellen. Im unteren Diagramm werden die Fourierdaten on-the-fly angezeit. Real- und Imaginäranteil. Der erste Imaginäranteil ist dort der Offset und der erste Realanteil ist immer gleich Null.

Das deckt sich ebend NICHT mehr mit deiner Beschreibung, dass sich dort irgend etwas konjugiert komplex spiegelt. Mir ist diese konjugiert komplexe Spiegelung wohl vom hörensagen bekannt, ich kenne diese jedoch nur daher, wenn man den Input spiegelt und nur die Imaginärwerte betrachtet, denn dann sind die wirklich symmetrisch gespiegelt.

s.o.

Um den Quelltext anzuschauen, klickt man im Programm auf die entsprechenden Dateikarten (Tabs).

Ja, das denke ich auch.

Sehe ich auch so.

Phaseninformationen ergeben sich doch erst durch das Verhältnis von Realanteil zu maginäranteil. Sprich arctan(imaginär/realteil).

MfG,

Markus

Hopsala. Das nehm ich zurück. Es spiegelt sich alles beim Wert 5 (gezählt von 1 bis 8, also im Diagramm Wert 4). Jetz ergibt das Ganze endlich Sinn *g*

MfG,

Markus

Ich habe kein Windows, und es ist vollkommen schwachsinnig, irgendein übersetztes Executable reinzustellen, daß Du offensichtlich selbst gestrickt hast und das bestenfalls fehlerhaft funktioniert.

Wie soll man denn da, bitteschön, erraten, was Du falsch machst?

Dann zeigt Dein Programm halt Unfug an. Und? Programme, die Unfug anzeigen, kann jeder schreiben.

Wie wäre es denn, wenn Du mal rausrückst, was Du da eigentlich rechnest? So mit Formeln oder Beispiel_quell_code?

Da gibst Du nichts an. Du postest irgendein Executable, das eh nur unter einem unfreien, teuren und virenanfälligen Betriebssystem lauffähig wäre, und daraus kann man im Leben nicht schließen, was das Programm treibt.

Wie wäre es, wenn Du den Quelltext nicht so aufwendig verstecktest?

Es gibt keinen Imaginäranteil bei den Werten 0 und 2048. Die Definition der DFT ist

N - 1 ---- 1 \ -2 pi k n f = - * > exp(---------) a k N / N n ---- n = 0

bei k=0 ist exp(-2 pi k n/N) = 1, bei k=N/2 ist es (-1)^n.

Ich habe keine Ahnung, was Dein Programm treibt, aber es kann nicht besonders sinnvoll sein.

"David Kastrup"

Na dann kannst du wohl ab jetzt nicht mehr mitreden. Dann dürfte die Diskussion jetzt wohl ziemlich einseitig werden *g*

Gut, wenn Du Dich entschieden hast, meine Antworten nicht mehr zu lesen, werde ich Deine Fragen nicht mehr lesen. Mir spart das Zeit.

Danke,

"David Kastrup"

^^^^^^^^^^^^^^^^^^^^^^^^

Ich frage mich, warum du mir antwortest, wenn du das, was ich tue, für "vollkommen schwachsinnig" hälst.

Markus,

Glaubst du wirklich dass dieser "Diskussionsstil" zielführend ist, wenn du etwas (hier Hilfe) von anderen erhalten möchtest?

Bei mir bist du nun jedenfalls die nächsten 90 Tage im Bozo-Bin. Also spar dir eine Antwort auf mein Posting, ich werde sie nicht mehr von meinem Newsreader angezeigt bekommen. Und auch eventuelle künftige Fragen von dir, sehe ich nicht mehr, selbst wenn ich dir dann vielleicht helfen könnte. Und glaube nicht, dass würde nur ich so machen.

cu

Er hält es ja nur für vollkomen schwachsinnig, wenn du eine EXE-Datei erzeugst, in der dann der Source-Code zu lesen ist, statt einfach nur den Source-Code zu veröffentlichen, was dann auch jemand nutzen kann, der Linux installiert hat.

Wahnsinn. Da hast du EINEN gefunden, der ganz offensichtlich weiß, wovon er redet, EINEN der auf deine Postings in immer ca.

"Frank Buss"

Den Code ohne das Programm zu zeigen macht deshalb keinen Sinn, weil man den ohne das Resultat zu sehen kaum nachvollziehen kann. Und ein Linux um es für Unix zu kompilieren habe ich momentan nicht installiert.

unit Unit1;

interface

uses Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms, Dialogs, ActnList, StdCtrls, ExtCtrls, ComCtrls, TeEngine, TeeFunci, Series, TeeProcs, Chart, U_Fourier1;

type TForm1 = class(TForm) PageControl1: TPageControl; TabSheet1: TTabSheet; TabSheet2: TTabSheet; Panel1: TPanel; Panel2: TPanel; ScrollBar1: TScrollBar; ScrollBar2: TScrollBar; ScrollBar3: TScrollBar; ScrollBar4: TScrollBar; ScrollBar5: TScrollBar; ScrollBar6: TScrollBar; ScrollBar7: TScrollBar; ScrollBar8: TScrollBar; ActionList1: TActionList; Action1: TAction; Chart1: TChart; Series1: TLineSeries; TeeFunction1: TAddTeeFunction; Timer_start: TTimer; Label1: TLabel; Label2: TLabel; Chart2: TChart; Series2: TLineSeries; Series3: TLineSeries; TabSheet3: TTabSheet; Memo1: TMemo; Memo2: TMemo; procedure Action1Execute(Sender: TObject); procedure ScrollBar1Scroll(Sender: TObject; ScrollCode: TScrollCode; var ScrollPos: Integer); procedure Timer_startTimer(Sender: TObject); private { Private-Deklarationen } public { Public-Deklarationen } end;

var Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.ScrollBar1Scroll(Sender: TObject; ScrollCode: TScrollCode; var ScrollPos: Integer); begin ActionList1.Actions[0].Execute; end;

procedure TForm1.Action1Execute(Sender: TObject); Var DFT: Tdft8array; begin Series1.Clear; Series1.Add(-127+(255-ScrollBar1.Position)); Series1.Add(-127+(255-ScrollBar2.Position)); Series1.Add(-127+(255-ScrollBar3.Position)); Series1.Add(-127+(255-ScrollBar4.Position)); Series1.Add(-127+(255-ScrollBar5.Position)); Series1.Add(-127+(255-ScrollBar6.Position)); Series1.Add(-127+(255-ScrollBar7.Position)); Series1.Add(-127+(255-ScrollBar8.Position));

DFT[1]:= -127+(255-ScrollBar1.Position); DFT[2] := 0.0; DFT[3]:= -127+(255-ScrollBar2.Position); DFT[4] := 0.0; DFT[5]:= -127+(255-ScrollBar3.Position); DFT[6] := 0.0; DFT[7]:= -127+(255-ScrollBar4.Position); DFT[8] := 0.0; DFT[9]:= -127+(255-ScrollBar5.Position); DFT[10] := 0.0; DFT[11]:= -127+(255-ScrollBar6.Position); DFT[12] := 0.0; DFT[13]:= -127+(255-ScrollBar7.Position); DFT[14] := 0.0; DFT[15]:= -127+(255-ScrollBar8.Position); DFT[16] := 0.0;

DFT8(DFT,8,1);

Series2.Clear; Series2.Add(DFT[2]); Series2.Add(DFT[4]); Series2.Add(DFT[6]); Series2.Add(DFT[8]); Series2.Add(DFT[10]); Series2.Add(DFT[12]); Series2.Add(DFT[14]); Series2.Add(DFT[16]);

Series3.Clear; Series3.Add(DFT[1]); Series3.Add(DFT[3]); Series3.Add(DFT[5]); Series3.Add(DFT[7]); Series3.Add(DFT[9]); Series3.Add(DFT[11]); Series3.Add(DFT[13]); Series3.Add(DFT[15]);

end;

procedure TForm1.Timer_startTimer(Sender: TObject); begin Timer_start.Enabled := FALSE; ActionList1.Actions[0].Execute; end;

end.

unit U_Fourier1;

interface

CONST nn=8; nn2=2*nn; (* 2*nn *) TYPE Tdft8array = ARRAY [1..nn2] OF real;

PROCEDURE DFT8(VAR data: Tdft8array; nn,isign: integer); FUNCTION sngl(x:real):real;

implementation

FUNCTION sngl(x:real):real; BEGIN sngl := x; END;

PROCEDURE DFT8(VAR data: Tdft8array; nn,isign: integer); (* Programs using routine FOUR1 must define type TYPE gldarray = ARRAY [1..nn2] OF real; in the calling routine, where nn2=nn+nn. *) VAR ii,jj,n,mmax,m,j,istep,i: integer; wtemp,wr,wpr,wpi,wi,theta: double; tempr,tempi: real; BEGIN n := 2*nn; j := 1; FOR ii := 1 TO nn DO BEGIN i := 2*ii-1; IF (j > i) THEN BEGIN tempr := data[j]; tempi := data[j+1]; data[j] := data[i]; data[j+1] := data[i+1]; data[i] := tempr; data[i+1] := tempi END; m := n DIV 2; WHILE ((m >= 2) AND (j > m)) DO BEGIN j := j-m; m := m DIV 2 END; j := j+m END; mmax := 2; WHILE (n > mmax) DO BEGIN istep := 2*mmax; theta := 6.28318530717959/(isign*mmax); wpr := -2.0*sqr(sin(0.5*theta)); wpi := sin(theta); wr := 1.0; wi := 0.0; FOR ii := 1 TO (mmax DIV 2) DO BEGIN m := 2*ii-1; FOR jj := 0 TO ((n-m) DIV istep) DO BEGIN i := m + jj*istep; j := i+mmax; tempr := sngl(wr)*data[j]-sngl(wi)*data[j+1]; tempi := sngl(wr)*data[j+1]+sngl(wi)*data[j]; data[j] := data[i]-tempr; data[j+1] := data[i+1]-tempi; data[i] := data[i]+tempr; data[i+1] := data[i+1]+tempi END; wtemp := wr; wr := wr*wpr-wi*wpi+wr; wi := wi*wpr+wtemp*wpi+wi END; mmax := istep END END;

end.

Tja; eventuell hilft Dir ein völlig anderer Ansatz?

Du bereibst offensichtlich so etwas ähnliches wie die Spektroskopiker in der Chemie. Und dort gibt es schon jahrelan eine Berechnungs-Methode um aus der Abfolge von verrauschten äquidistanten Analog-Werten die Peak-Lage herauszu- bekommen. Wenn Du eine (chehmische) Uni-Bibliothek in Reichweite hast, fin- dest Du eventuell den Original-Text:

Savtzky, A. and Golay, M.J.E. 1964, Analytical Chemistry Vol. 36 Pages 1627 - 1639.

Original in FORTRAN programmiert; aber mit etwas Google-Hilfe findet man Quelltexte in C++

Gruss, Holger