Forum: Digitale Signalverarbeitung / DSP / Machine Learning Genaue Frequenzbestimmung nach FFT


von Frank (Gast)


Lesenswert?

Hi,

welche Moeglichkeiten gibt es nach einer FFT die excacte Frequenz 
auszurechnen. Ich habe das Problem, dass ich bei der FFT berechnung 
Zeitlich eingeschrenkt bin und somit nur ein ungenauses Spektrum 
bekomme. Ich arbeite zurzeit mit einer 4096 Punkte FFT und einer 
Samplerate von 10 MSPS.
Aus der FFT bekomme ich nur die grobe Frequenz, muss nun danach aber die 
Frequenz hochgenau bestimmen.
Hat jemand Vorschlaege?

Gruss
Frank

von Detlef _. (detlef_a)


Lesenswert?

Du machst zwei FFTs, die überlappend versetzt sind, meinetwegen 200 
samples oder so. Dann kuckst Du an der Frequenz, deren Betrag am 
höchsten ist: Dort die Phasendrehung ansehen, damit kannst Du die genaue 
Frequenz berechnen. Die FFTs brauchen nen Fenster, Kaiser 16 z.B..

Cheers
Detlef

von Gerhard (Gast)


Lesenswert?

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

von Kleinster Gemeiner Nenner (Gast)


Lesenswert?

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

von Michael L. (Gast)


Lesenswert?

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

von Frank (Gast)


Lesenswert?

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.

von Detlef _. (detlef_a)


Lesenswert?

Das untenstehende Matlab-script bestimmt mittels versetzter FFTs die 
genaue Frequenz und deren Phase. Die Frequenz ist dabei nicht an das 
Frequenzraster der FFT angepaßt. Man benötigt also keine hochauflösenden 
FFTs.

Nix zero padding, nix interpolieren, das geht besser.

Cheers
Detlef

clear
n=1152;
ff= 941/8000;
%sigsoll=5000*cos(n*ff*2*pi*(0:n-1)/n)+ ...
%        5000*cos(n*ff1*2*pi*(0:n-1)/n);
sigsoll=5000*cos(n*ff*2*pi*(0:n-1)/n));
n=length(sigsoll);
plot(abs(fft(sigsoll)),'.-')

%plot(sigsoll,'b.-');

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);

%plot(20*log10(abs(sp1)),'.-')
[m,ind] = max(abs(sp1(1:length(sp1)/2)));

%angle(sp1(ind))
dw=angle(sp2(ind)/sp1(ind));

%wb=ff/2*2*pi/512
fg=(dw/d/(2*pi));

sp=fft(cos(2*nn*fg*pi*(0:nn-1)/nn).*wind);
%angle(sp(ind))
ph=angle(sp1(ind)/sp(ind));

return

von Frank (Gast)


Lesenswert?

@ 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?

von Frank (Gast)


Lesenswert?

@ 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.

von Achim (Gast)


Lesenswert?

@ 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.

von Detlef _. (detlef_a)


Lesenswert?

Ich habe die entscheidenden Punkte nochmal kopiert und kommentiert

[m,ind] = max(abs(sp1(1:length(sp1)/2)));
% ind ist der Index, an dem das Betragsspektrum maximal wird

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

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

sp=fft(cos(2*nn*fg*pi*(0:nn-1)/nn).*wind);
% Das ist ein Signal mit der genauen Frequenz, gefenstert und 
transformiert

ph=angle(sp1(ind)/sp(ind));
% Die Phase des Signals ist der Differenzwinkel, wieder am 
Betragsmaximum
% Das geht auch mit Mischen der Signale, das spart die dritte FFT

Fachausdruck für das Verfahren ist mir nicht bekannt, in einem Buch habe 
ich das so noch nicht gesehen. Es ist aber nicht so genial, dass außer 
mir noch niemand drauf gekommen wäre ;-))


Cheers
Detlef

von Frank (Gast)


Lesenswert?

@ 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

von Marius W. (mw1987)


Lesenswert?

Soweit ich das Skript verstehe, arbeitet es mit normierten Frequenzen. 
Das bedeutet du bekommst als Frequenz etwas auf die Abtastfrequenz 
normiertes heraus. Multipliziere dein Ergebnis mit der Abtastfrequenz 
und du hast die Frequenz die du suchst.

MfG
Marius

von Frank (Gast)


Lesenswert?

@ Marius

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

von Detlef _. (detlef_a)


Lesenswert?

>>Jedenfalls bekomme ich mit dieser Gleichung keine Frequenz raus. Ich
>>kann bei dem Script auch nix derartiges entdecken wo die Frequenz sein
>>soll. (fg ?!?!)

klar ist fg die gesuchte Frequenz, das siehst Du schon an dieser Zeile:

sp=fft(cos(2*nn*fg*pi*(0:nn-1)/nn).*wind);

Und diese ist gleich der reingesteckten Frequenz, das ist ff=941/8000 . 
Wie Marius richtig sagt, ist das auf die Abtastfrequenz normiert, das 
ist bei digitalen Signalen immer der Fall.

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

Sach doch mal, was irgendwie nicht klappt.

Cheers
Detlef

von Frank (Gast)


Lesenswert?

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.

von Frank (Gast)


Lesenswert?

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??

von Detlef _. (detlef_a)


Lesenswert?

>>Die zweite FFT wurde um 200 Samples verschoben.

Wenn Du zuviel verschiebst, dreht der Frequenzzeiger in dieser Zeit 
mehrmals einen vollen Kreis, Du bekommst also eine Unsicherheit von 2pi 
rein. Es müssen nicht zwei samples sein, aber der Zeiger darf auch nicht 
mehr als eine volle Runde drehen (oder Du kennst Deine Frequnz ungefähr 
und kannst das einbeziehen)

Cheers
Detlef

von Frank (Gast)


Lesenswert?

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?

von Detlef _. (detlef_a)


Lesenswert?

Bei 10MHz geht das bis 2.5MHz, also Nyquist/2, wenn ich mich recht 
entsinne. Ab dann kriegst Du eine 'negative' Frequenz, die man noch 
entsprechend umlügen muss.

Cheers
Detlef

von Frank (Gast)


Lesenswert?

@ 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.

von pumpkin (Gast)


Lesenswert?

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:


von Anja (Gast)


Lesenswert?

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

von Anja (Gast)


Lesenswert?


von Frank (Gast)


Lesenswert?

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?

von Detlef _. (detlef_a)


Lesenswert?

>>Die Methode von Detlef funktioniert erst ab dem dritten oder vierten FFT 
Koeffizienten sehr genau.

Ja, ist möglicherweise nen Leckeffekt des Fensters, muß ich nochmal 
draufkucken.

Du kannst aber Deine samples für die niedrigen Frequenzen einfach 
dezimieren: Zehn Abstastwerte zu einem addieren und damit die samplerate 
zehnteln, z.B. . Damit sollten dann die niedrigen Frequenzen auch 
bestimmbar sein.

Cheers
Detlef

von Frank (Gast)


Lesenswert?

@ 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

von Detlef _. (detlef_a)


Lesenswert?

Wenn Du die Überabtastung reduzierst, also samples wegwirfst oder sie 
zusammenaddierst, rückt der peak ins Hochfrequente und vermeidet damt 
das Problem bei den tiefen Frequenzen. Das hat mit der Länge der FFTs 
erstmal nix zum Tun, wenn Du das Verhältnis Signalfrequenz zu Abtastrate 
größer machst wirds in der FFT hochfrequenter, unabhängig von deren 
Länge. Oder habe ich da noch was nicht verstanden?

Cheers
Detlef

von Frank (Gast)


Lesenswert?

@ 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

von ralf (Gast)


Lesenswert?

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.

von Frank (Gast)


Lesenswert?

@ 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.

von ralf (Gast)


Lesenswert?

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.

von Detlef _. (detlef_a)


Lesenswert?

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

Ja, so is das. Die Abtastrate runtersetzen und die Meßzeit hochsetzen 
sollte die Bestimmung für niedrige Frequnzen ermöglichen.

Ich habe das nochmal mit ner Hilbert Transformation probiert, das geht 
auch, wird aber auch zu tiefen Frequenzen schlechter und ich habe keine 
Erfahrung, wie das auf z.B. Meßrauschen reagiert.

Cheers
Detlef

clear
a=sin(1.2*2*pi*(0:511)/512);
c=hilbert(a);
cd=c(2:end)./c(1:end-1);
d=0;
length(a)*mean(angle(cd(256-d:256+d)))/2/pi
return

von Martin L. (Gast)


Lesenswert?

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.

von Hannes (Gast)


Lesenswert?

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 ? 
:-)

von dspgast (Gast)


Angehängte Dateien:

Lesenswert?

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

von Detlef _. (detlef_a)


Lesenswert?

Danke für die Blumen.

>>Eine Idee hätte ich noch:

Ja, so sollte das auch gehen, der Vektor des runtergemischten Signals 
dreht sich mit der Differenzfrequenz


Das Verfahren oben war mein Erkenntnisstand damals, Frequenzbestimmung 
eines unbekannten Signals geht noch intuitiver und einfacher:

Beitrag "Re: Phasenverschiebung messen (us Bereich)"

Cheers
Detlef

von Gue5t87 (Gast)


Lesenswert?

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!

von Christoph db1uq K. (christoph_kessler)


Lesenswert?

Auf der UKW-Tagung gab es einen Vortrag zum Thema. Eine Parabel wird 
durch die drei nächstliegenden Punkte der FFT um die gesuchte Frequenz 
gelegt und deren Maximum genomman.
Der Verfasser Paul Boven zititert dazu eine Doktorarbeit am CERN:
http://mgasior.web.cern.ch/mgasior/pap/index.html
M. Gasior, Improving Frequency Resolution of Discrete Spectra

von Gue5t87 (Gast)


Lesenswert?

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

von Detlef _. (detlef_a)


Lesenswert?

>>Der Code von Detlef ist wirklich spitze.
Yo, THX :))))

Wenn Du ne höhere Frequnz als Nyquist/Halbe hast, hat der Zeiger zu weit 
gedreht: Die beiden FFTs nur um einen Abstawert verschieben, gemoddeter 
Code anbei.

Cheers
Detlef


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

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

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

nn=512; % mit zwei 512er FFTs
d=1;
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

von Gue5t87 (Gast)


Lesenswert?

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

von T. (Gast)


Lesenswert?

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.,

von Christoph db1uq K. (christoph_kessler)


Lesenswert?

ja guten Tag erstmal
Das war zwar schon vor drei Jahren, aber zum Nachlesen: ISBN 
9783540634430
Es gibt auch Matlab-Bücher von einem Josef Hoffmann, aber hier stimmt 
der Vorname tatsächlich.
Wobei Amazon die Ü-Tüpfelchen wegläßt

: Bearbeitet durch User
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.