Genaue Frequenzbestimmung nach FFT

Hi

das Verhältnis von Abtastrate und Anzahl der Messpunkte/FFT bestimmt die 
Auflösung der FFT. Mach das Verhältnis passend, dann passt auch deine 
Wunsch-Auflösung.

Gerhard

Detlef, kannst Du das bitte mal genauer erklären?

Hallo Frank,

es gibt verschiedene Möglichkeiten:

- Du bist umso genauer, je länger Dein Meßsignal ist. Im Prinzip kannst 
Du also einfach Nullen an Dein Signal anhängen und dann erneut eine FFT 
durchführen. Dadurch wird der Peak deutlicher aufgelöst. Dieses 
Verfahren heißt "zero padding".

- Eine anderere Möglichkeit besteht darin, daß Du Dein Spektrum im 
Bereich des Peaks interpolierst. Typischerweise verwendest Du dazu ein 
Polynom 2. Ordnung oder eine Gaußfunktion. In beiden Fällen kannst Du 
den Peak aus den Fit-Parametern symbolisch ausrechnen.

- Eine weitere Möglichkeit ist die Messung der Momentanfrequenz mithilfe 
der Hilberttransformation. Dadurch erhältst Du die Phase des Signals, 
die Du nur abzuleiten brauchst, um die Frequenz zu erhalten.
Das funktioniert insbesondere gut bei schmalbandigen Signalen mit wenig 
Rauschen. Unter guten Bedingungen reicht weniger als eine Signalperiode 
zur Frequenzbestimmung aus.

Grundsätzlich gilt:
Jede Änderung (sowohl das zeitliche Begrenzen durch das Ausschalten der 
Messung, als auch das Hinzufügen von Nullen) verändert allerdings Dein 
Signal minimal.


Gruß,
  Michael

Danke fuer die Antworten. Ich haette wohl dazusagen sollen, dass ich 
keine Zeit habe um eine zweite FFT durchzufuehren, gar ein Zero Padding, 
oder mehr Werte Sampeln.
Die Vorschlaege von Michael klingen somit super. Werde ich sofort mal 
ausprobieren.
Soweit ich Detlef _a verstanden habe, soll man die Phasendrehung 
zwischen den beiden Peak-Frequenz-Koeffizienten, die als ergebniss der 
FFTs entstehen, zur genauen Frequenzberechnung heranziehen. Mich wuerde 
es erlichgesagt auch interessieren nach welcher Gleichung.

@ Detlef

Danke fuer das Script. Werd mich da mal einarbeiten.

Ich musste leider schon feststellen, dass eine Interpolaration nicht gut 
funktioniert, sogar im gegenteil, komischerweise liefert mir das 
interpolierte Spektrum eine groessere Abweichung als das nicht 
interpolierte :) Es tendiert grafisch gesehen in die entgegengesetzte 
richtung. Hat das vlt. etwas mit dem Ergebniss einer FFT zu tun, da man 
ja auch ein gespiegeltes Spektrum bekommt? Muss man also die ersten 
Koeffizienten zur Spiegelachse oder die letzteren koeffizienten zur 
interpolaration heranziehen?

@ Detlef

bevor ich mir noch den ganzen Kopf zerbreche frage ich lieber nochmal 
nach, da ich erlichgesagt noch nie mit MatLab gearbeitet habe.
ALso ich mache zwei versetzte FFTs. Die Engangssamples vorher natuerlich 
Fenstern und dann in die FFTs. Nun suche ich jeweils das Maximum bzw. 
den Peak beider FFTs, nachdem ich den Betrag der Koeffizienten bestimmt 
habe. So nun habe ich die zwei entsprechenden FFT koeffizienten. Was 
dann passiert ist mir noch nicht ganz klar.

@ Detlef

Danke für die Idee und das Script. Ich benutze zwar kein Matlab sondern 
Scilab. Dein Script habe ich mit kleineren Änderungen zum laufen 
bekommen.

Genau wie Frank habe ich allerdings das Verfahren an sich noch nicht 
ganz verstanden. Gibts dafür einen eigenen Fachausdruck wo ich mal 
schauen kann? Oder erklär doch bitte mal noch wie die Weiterverarbeitung 
der beiden FFT-Matrizen funktioniert.

Könnte das Prinzip sehr gut gebrauchen, da ich eine genaue, permanente 
Frequenzbestimmung in Echtzeit auf einem Mikrocontroller umsetzen muss.

@ Detlef

Ich sage mal so: das Script verstehe ich nun. Aber wie bekommste denn 
aus dieser Gleichung eine Frequenz raus? : fg=(dw/d/(2*pi));

Um das in Worte zu fassen: Man teil den Differenzwinkel durch die 
verschobenen Samples. In deinem Beispiel oben hast du ja nur um 2 
Samples verschoben.
Jedenfalls bekomme ich mit dieser Gleichung keine Frequenz raus. Ich 
kann bei dem Script auch nix derartiges entdecken wo die Frequenz sein 
soll. (fg ?!?!)

Benoetigt man das dritte Signal und FFT nur fuer die Phase?


Gruss
Frank

@ Marius

Das dachte ich zuerst auch, aber das klappt irgendwie auch nicht.

Ich sehe es ja auch, dass fg die Frequenz sein muss, aber irgendwie 
bekomme ich da nix richtiges raus.
Also ich habe zwei ueberlappende FFTs. Die zweite FFT wurde um 200 
Samples verschoben. Im script verschiebst man nur um 2 Samples. Die 
verschiebung um x Samples ist doch im prinzip egal, oder?
Jedenfalls liefern mir beide FFTs ein nahezu identisches Spektrum (soll 
ja auch so sein :)). Danach bilde ich den Differenzwinkel. Aber vlt. 
mache ich da generell was falsch: ich nehme den arctan (b/a) beider 
Koeffizienten und bilde die Differenz. Hab auch probiert zuerst den 
Quotienten der Koeffizienten zu bilden und dann den Winkel auszurechnen 
(also so wie im script). Beide Methoden ergeben den gleichen Winkel. Ich 
berechne alles in Bogenmass.
Die Formel fg leuchtet mir ja auch ein. Um das zu konkretisieren:
Ich habe eine Sinusfrequenz von 394631 Hz simuliert. Die beiden 
versetzten FFTs ergeben ein richtiges Spektrum. Durch interpolation habe 
ich sogar nur ne Abweichung von gerade mal 0,07%, das aber nur am Range. 
Der Differenzwinkel der beiden FFT Peak Koeffizienten eingesetzt in die 
fg formel im Script und multipliziert mit der Samplefrequenz liefert mir 
aber ein Ergebnis von 5368 Hz.

hmm hab gerade versucht nur um 2 samples zu verschieben und bekomme 
tatsaechlich die richtige frequenz raus :). Lag es wirklich dadran? Muss 
es unbedingt nur um 2 Samples verschoben sein??

Bei den Simuationen sehe ich nun auch, das es mit hohen Frequenzen nicht 
funktioniert. Also ich habe versucht ein Signal von ca. 4 MHz und einer 
Samplerate von 10 MHz zu berechnen, jedoch funktioniert das 
komischerweise nicht. Wenn ich die 4 MHz gegen 300 KHz tausche, dann 
funktioniert diese Methode wieder fantastisch. Gibt es da eine Regel 
welche Frequenzen mit dieser Methode betrachtet werden koennen?

@ Detlef

Habs mal durchsimuliert: Bis Nyquist/2 braucht man den Differenzwinkel 
nicht umzurechnen. Wenn der Winkel PI ueberschreitet, dann geht er 
wieder ins negative. Man muss also nur die Differenz zu PI wieder mit PI 
addieren und schon kann man auch die Frequenzen ueberhalb Nyquist/2 
berechnen.

Vielen Dank fuer deine Hilfe und deine ausfuerlichen Erklaerungen. Das 
hat mir sehr sehr geholfen.

Das sollte nach meiner Vorstellung eigentlich nicht allein von der 
Abtastfrequenz abhängen. Die Zielfrequenz muss natürlich kleiner Nyquist 
sein, zusätzlich sollte es von der Verschiebung dt abhängen. Die 
Periodendauer der gesuchten Frequenz muss länger sein als dt (im 
Spezialfall genau gleich, siehe unten). Auf deutsch: Die FFT's müssen 
sich überlappen (wenigestens aneinander angrenzen). Anderenfalls bekommt 
man - wie Detlef schon schrieb - die Unsicherheit von (N*)2*Pi rein.

Beispiel: Im Extremfall sucht man die Frequenz genau bei Nyquist, hat 
also nur zwei Abtastpunkte pro Periode; die Verschiebung dürfte dann 
maximal 2 Samples betragen, das Ermitteln der Frequenz sollte auch dann 
klappen. Der besagte Spezialfall: 2-Punkt-FFT. Etwa so:

Frank schrieb:
> Hat jemand Vorschlaege?

In der Artikelserie von Lothar Wenzel in der Elektronik gab es in Teil 5 
der Serie eine Formel um die Auflösung für Sinusförmige Signale exakt 
auszurechnen.

Gruß Anja

Sorry habe den Link vergessen:
http://digital.ni.com/worldwide/germany.nsf/3528e4c55a5f22d4852568e100532f33/5f0bb6c216c0df648025680100348c48?OpenDocument

Gruß Anja

Noch ein Zusatz: Die Methode von Detlef funktioniert erst ab dem dritten 
oder vierten FFT Koeffizienten sehr genau. Vorher kann man es leider 
nicht anwenden. Das bringt mich zu einem weiteren Problem. Da ich bei 
einer 512 Punkte FFT und einer Samplerate von 10 MSPS erst bei ca. 40 
kHz die Methode von Detlef anwenden kann, ich brauche aber auch einer 
sehr sehr genau Frequenzerkennung unterhalb 40 kHz. Kennt jemand eine 
Möglichkeit auch diese Frequenzen (also < 40 kHz, bei einer Samplerate 
von 10 MSPS) zu erfassen ohne die Samplerate oder die Auflösung der FFT 
zu verändern?

@ Detlef

Wenn ich aber die 10 Samples zusammenfuege und somit die Samplerate 
zehntel, dann brauche ich doch 10 mal mehr samples um die FFT 
durchzufuehren? Das bringt mich ja auch nicht weiter. Oder habe ich da 
was falsch verstanden?

Gruss
Frank

@ Detlef

Die Abhaengigkeit moegen ja stimmen, aber mit was soll ich denn die FFT 
fuettern? Wenn ich 512 Samples habe und dann ein zehntel wegschmeisse 
oder aufaddiere, dann erhalte ich ja letztendlich 51 Samples. Damit 
luege ich ja vor, dass ich mit einem zehntel der echten Samplerate 
gesamplet habe. Aber da ich ja nur noch 51 Samples habe, bleibt die 
Aufloesung der FFT gleich. Ich kann den Peak nur ins "hochfrequente" der 
FFT ruecken wenn ich mehr Samples habe, aber woher soll ich diese 
herbekommen?
Und es ist ebend nicht unabhaengig von der laenge der FFT. Wenn ich mit 
10 MSPS eine 512 Punkte FFT durchfuehre, habe ich die gleiche aufloesung 
wie mit 1 MSPS und einer 51 Punkte FFT. Das heisst, dass auch der 
gesuchte Peak (niedriger Frequenz) in den gleichen Koeffizienten faellt. 
Und falls dieser halt in den ersten 3 Koeffizienten faellt, dann ist der 
Fehler halt zu gross. In meinerm Fall verliere ich die ersten 60 KHz, 
die ich unbedingt auch brauche.
Auch wenn ich die FFT aufloesung mit Zero Padding vergroessere habe ich 
das grundlegende Problem, dass die Nyquist bedingung nicht erfuellt 
wird. Soweit ich das verstanden habe, muss ich ja nicht nur mit 
midestens doppelter Frequenz sampeln (was in diesem fall erfuellt ist), 
sondern auch mindestens eine Periode des Signal erfassen. Genau das wir 
ja nicht erfuellt unterhalb der ersten 40 kHz (deshalb auch der hohe 
fehler).
Deshalb habe ich die letzten Tage vergeblich versucht deine Vorschlaege 
(seoweit ich sie verstanden habe :)) durchzusaetzen.

Ich konkretisier mal meine Versuche: Ich habe 512 Samples mit 10 MSPS 
aufgenommen. Meine FFT Aufloesung (512 Punkte) entspricht somit knapp 20 
kHz. Der erste Koeffizient liefert mir keine Phase --> ich verliere 
deshalb die ersten 20 kHz. Der zweite und dritte Koeffizient ist leider 
sehr ungenau. Deshalb habe ich nur ab ca. 60 kHz eine sehr gutes 
Ergebnis. Wenn ich nun deinem Vorschlag folge nehme ich von den Samples 
jedes zehnte und zehntel somit meine Samplerate. Wenn ich aber nun eine 
FFT durchfuehre (egal ob mit 51 Samples oder 512 Samples (Zero Padding 
ab Sample 51)), faellt mir das Ergebnis bei einem Signal unterhalb 20 
kHz, immer in den ersten Koeffizienten und ist somit nicht verwertbar.

Habe ich deine Vorschlaege richtig verstanden?

Gruss
Frank

Warum hasst du die Länge der FFT von 4096 auf 512 verringert? Damit das 
Verfahren funktioniert muss die FFT mindestens drei Perioden deines 
Eingangssignal enthalten.

@ ralf

was heisst verringert. Mit 512 Samples kann ich von 60 kHz bis 4 MHz die 
Signalfrequenz unheimlich genau bestimmen. Natuerlich sehe ich ein, dass 
wenn ich mehr Samples nehme, also laenger Sample, ich auch die unteren 
Frequenzen sehr genau bestimmen kann. Jedoch auch bei 4096 Samples kann 
ich Frequenzen unterhalb 6 kHz nicht bestimmen. Aber diese brauche ich 
auch.
Da muss es doch bestimmt einen weg geben, wie man ohne ehoehung der 
Samples und ohne das Aufnehmen einer oder mehreren Periode die Frequenz 
errechnen kann? Weil wenn ich z.B. 50 Hz Signal habe, dann muesste ich 
ja theoretisch millionen von Samples nehmen und eine FFT durchfuehren, 
und das bringt mir nix.
Deswegen versuche ich ja vergeblich eine Loesung zu finden. Was mir 
bislang in den Sinn kommt, ist nach dem Vorschlag von Detlef, Samples 
wegzuschmeissen und damit die Samplerate verringern. Aber es ist 
irgendwie immer noch nicht das wahre.
Deshalb, falls jemand noch eine andere Methode kennt, waere ich sehr 
dankbar.

Wenn du keine Annahmen für dein Eingangssignal treffen kannst kommst du 
nicht am Zeit-Bandbreite-Gesetz vorbei.

http://storage.sk.uni-bonn.de/Milca/ssv/content/ssv_s431.xhtml

http://gdr-isis.org/tftb/tutorial/node25.html

Das bedeutet das bei vorgegebener minimaler Frequenz das Eingangssignal 
zwangsläufig auch über einen minimalen Zeitbereich betrachtet werden 
muss.

Hallo,

es gibt eine prinzipielle Grenze zur messabren Frequenzgenauigkeit in 
Abhängigkeit der Messzeit. Diese nennt man Heisenbergsche 
Unschärferelation. Und deswegen bringt auch Zero-Padding oder 
Interpolation des Frequenzspektrums nur bis zu einem gewissen Punk (IMHO 
das anhängen von der selben Anzahl von 0en an das Messsignal) etwas.
Es gibt darüber ein sehr schönes Buch von Rüdiger Hoffmann mit dem Titel 
"Signalanalyse und-erkennung: Eine Einführung für 
Informationstechniker". Dort ist auch beschrieben wie man ggf. einen 
Kompromiss mit Hilfe der Wavelettransformation finden kann.

Ansonsten für den mathematisch interessierten:

http://en.wikipedia.org/wiki/Uncertainty_principle#Uncertainty_theorems_in_harmonic_analysis

Die Frage taucht ja aber in abgewandelter Form immer wieder mal auf. 
Vielleicht sollte man mal eine Seite dazu schreiben....

Viele Grüße,
 Martin L.

Hätte ich garnicht gedacht dass die HUR sogar in die Mathe reinspielt.

Martin Laabs schrieb:
> Rüdiger Hoffmann
und dass der solche Bücher schreibt (?) Macht der sonst nicht Comedy ? 
:-)

Danke an Deltef_a für diese geniale Methode!
Im Anhang übrigens meine Scilab-Übersetzung davon, falls es jemand 
braucht.
Ich bin sehr überrascht von der Genauigkeit dieser Methode: Man kann die 
Frequenz auf  941.001 Hz setzen und bekommt das immer noch auf die 
letzte Stelle genau heraus.
Eine Idee hätte ich noch: Kann man, wenn die Frequenz hinreichend genau 
bekannt ist und es nur auf den genauen Wert ankommt, nicht auch einfach 
eine "fft mit nur einem Bin" machen, also die Eingangsdaten mit I und Q 
runtermischen (mit der bekannten, groben Frequenz), jeweils aufsummieren 
und dann den Winkel berechnen?

Gruß, dspgast

Der Code von Detlef ist wirklich spitze. Jedoch bei einer Abtastfrequenz 
von fs = 48 kHz, funktioniert dieser nur bis eine Frequenz von f = 12 
kHz.

Woran liegt das? Ich benötige eine genaue Frequenzbestimmung für f = 20 
kHz.


Hier nochmal der Code:

clear; close all;
fs = 8000;
t = 0:1/fs:1;
f = 3654 % moeglich bis 12kHz

sigsoll=sin(2*pi*f*t);
n=length(sigsoll);

figure
plot(abs(fft(sigsoll)),'.-')

nn=512; % mit zwei 512er FFTs
d=2;
wind=kaiser(nn,16).'; % 3e-7

sp1=fft(sigsoll(1   : nn).*wind);
sp2=fft(sigsoll(1+d:nn+d).*wind);


figure
plot(20*log10(abs(sp1)),'.-')
[m,ind] = max(abs(sp1(1:length(sp1)/2)));
               % ind ist der Index, an dem das Betragsspektrum maximal 
wird


%angle(sp1(ind))
dw=angle(sp2(ind)/sp1(ind));
                % dw ist der Differenzwinkel, um den sich der komplexe
                % Frequenzzeiger bei einer
                % Verschiebung um d Abtastwerte gedreht hat


%wb=ff/2*2*pi/512
fg=(dw/d/(2*pi));
            % dw/d ist der Drehwinkel bezogen auf einen Abtastwert.
            % Haette man einen Differenzwinkel von pi, wäre man genau 
auf der
            % Nyquistfrequenz (2 Abtastwerte für 2*pi)
            % Bei 2*pi Differenzwinkel wäre man genau auf der 
Abtastfrequenz.
            % fg normiere ich entgegen der Konvention auf 2*pi, also auf
            % die Abtastfrequenz

freque = fg * fs



Hoffe Ihr könnt mir weiterhelfen. Vielen Dank im voraus!

Ok dankeschön. Aber irgendwie hilft mir das nicht wirklich weiter, sorry 
:-/

Tatsächlich! - Jetzt funktionierts! Dickes Dankeschön! :)

Martin Laabs schrieb:
> Es gibt darüber ein sehr schönes Buch von Rüdiger Hoffmann mit dem Titel
> "Signalanalyse und-erkennung:
Da würde ich gerne mal den Einleitungstext sehen:
"Hallo erstmal, ich weiss garnicht, ob Sie's wussten, aber 
wavelet-Transformationen lassen sich genial nutzen, um FFTs 
auszurechnen. Meine Bekannte, die Birte, hatte mal wieder das Problem, 
dass sie nach einer Fast Food Transaktion einen fetten Arsch bekommen 
hat, u.s.w.,