FFT von Kosinus-Funktionen funktioniert nicht

[pyhill00]

Hi zusammen
ich versuche vier Cosinus-Funktionen zu FFTen aber das Ergebnis ist komisch.

Hier ist das Ergebnis und der Code

https://imgur.com/a/FP3ylak

Code: Alles auswählen

import numpy as np
import scipy.fftpack
from scipy.fftpack import fftfreq
from scipy.fft import fft
from scipy.fft import fft2
from scipy.fft import fftn

from scipy.signal import blackman, flattop, boxcar
import matplotlib.pyplot as plt
import random

plt.rcParams['figure.figsize'] = [9, 5]

gamma = 2.8e6 # short for gamma/2*pi, in Hz, [Hz/G]
tau = 1e-6 # us;
G = np.linspace(0, 1.7e6, 2400) # G_max is 1.7G/um ~ 1.7e6G/m
k = gamma*G*tau # 1/nm


rx = np.linspace(404e-7,416e-7,4)
outer = np.outer(rx,k)

y = 0.5*np.cos(2*np.pi*outer) #+ np.random.normal(0, .1, k.shape)
yy = np.sum(y,axis=0)

plt.figure(0)
#
#plt.xlim(1.0e8, 1.3e8)
plt.plot(k, yy)


f = fftfreq(len(k), np.diff(k)[0])
yf = fft(yy)
plt.figure(1)
plt.xlim(0, 500e-7)
plt.plot(f[:G.size//2], np.abs(yf[:G.size//2]))

r_max = f[:G.size//2][np.argmax(np.abs(yf[:G.size//2]))]
amplitude_max = np.abs(yf[:G.size//2]).max()

Ich habe eher 4 Spitzen erwartet. Ich hab rx zu rx = np.linspace(404e-7,4016e-7,4) geändert und ich bekomme

https://imgur.com/a/ZjTedft

Irgendwie macht der Abstand alles kaputt.

[pyhill00]

Ok ich glaub, dass es Sinn ergibt. Aber kann mir jemand erklären, warum die Amplituden nicht gleich sind? Zb hier:

https://imgur.com/a/92x1nli

Das sind vier Kosinus-Funktionen und die FFT-Amplituden sind nicht gleich. Und wenn ich die Anzahl der Funktionen erhöhe, bekomme ich so ne Art Oszillation der Amplitude https://imgur.com/a/loyG8Gn

[__deets__]

Ist nur geraten, aber ich vermute mal numerische Probleme. Schau mal, ob die Probleme verschwinden oder variieren bei anderen Wertebereichen.

[pyhill00]

Ok danke. Ich habe vier Kosinus-Funktionen mit den Frequenzen 400e-3, 500e-3, 600e-3 und 700e-3, und ich versuche, die FFT von ihnen durchzuführen, aber unter der gebrauchten Zeit können die vier nicht unterschieden werden. Gibt es eine Möglichkeit, die Spitzen zu unterscheiden, ohne die tmax-Zeit von 1,76 und die Frequenzen zu ändern?

Code: Alles auswählen

    import numpy as np
    import scipy.fftpack
    from scipy.fftpack import fftfreq
    from scipy.fft import fft
    import matplotlib.pyplot as plt


    t = np.linspace(0,1.76,2400) 
    f = [400e-3, 500e-3, 600e-3, 700e-3] # these are the frequencies
    yy = 0

    for i in f:
        y = 0.5*np.cos(2*np.pi*i*t)
        yy = yy + y

    plt.figure(0)
    plt.plot(t, yy)


    f = fftfreq(len(t), np.diff(t)[0])
    yf = fft(yy)
    plt.figure(1)
    plt.plot(f[:t.size//2], np.abs(yf[:t.size//2]))
    plt.show()

[__deets__]

Also wenn ich dich richtig verstehe, dann willst du deine 4 Frequenzen im sub-1Hz-Bereich, aber mit nur 1.76 Sekunden Samplezeit? Ausser einer FFT faellt mir da nix ein, und die sieht ja besch..eiden aus bei der kurzen Laufzeit. Insofern: mir ist da nix bekannt.

[markus_]

Mal so ein spontaner Gedanke: Die FFT füllt den Datensatz ja bis zur nächsten Zweierpotenz mit Nullen auf (Zero-Padding). Das könnte das eigentliche Signal hier ziemlich verfälschen, weil ja nur sehr wenige Schwingungen drin sind aufgrund der kurzen Samplezeit. Hast du mal versucht das Ausgangszeitsignal vor der FFT zu Fenstern? Ansonsten würde ich das mal versuchen. Die Fensterung verhindert den abrupten Sprung des Signals auf Null und schafft einen weichen Übergang.

Edit: Ein Versuch wäre noch statt der FFT eine DFT zu verwenden. Diese benötigt nicht einen Datensatz mit Zweierpotenz-Länge, ist dafür aber weit langsamer. Aber Zeit scheint hier ja nicht das Problem zu sein, oder?

[pyhill00]

Mal so ein spontaner Gedanke: Die FFT füllt den Datensatz ja bis zur nächsten Zweierpotenz mit Nullen auf (Zero-Padding). Das könnte das eigentliche Signal hier ziemlich verfälschen, weil ja nur sehr wenige Schwingungen drin sind aufgrund der kurzen Samplezeit. Hast du mal versucht das Ausgangszeitsignal vor der FFT zu Fenstern? Ansonsten würde ich das mal versuchen. Die Fensterung verhindert den abrupten Sprung des Signals auf Null und schafft einen weichen Übergang.

Edit: Ein Versuch wäre noch statt der FFT eine DFT zu verwenden. Diese benötigt nicht einen Datensatz mit Zweierpotenz-Länge, ist dafür aber weit langsamer. Aber Zeit scheint hier ja nicht das Problem zu sein, oder?

Wie geht das mit einer DFT? FFT ist doch ein Algorithmus, um die DFT zu berechnen. Könntest du bitte mehr zu deiner Idee erläutern? Sie interessiert mich sehr.

[markus_]

DFT = diskrete Fouriertransformation
FFT = Fast Fourier Transformation

Der FFT-Algorithmus nutzt einige Symmetrien in der Berechnung der DFT aus und ist somit deutlich schneller als die herkömmliche diskrete FT. Die Bedingung ist aber, dass die Länge der Datenpunkte 2^n entspricht (also ...512,1024, 2048,...). In der Regel füllt der FFT Algorithmus die fehlenden Datenpunkte bis zur nächsten Zweierpotenz mit Nullen auf. Das ist sicher auch bei Python so. Das hat natürlich Einfluß aufs Ergebnis, ob das hier die Ursache sein kann weiß ich allerdings auch nicht, da stecke ich nicht tief genug in der Mathematik drin. Es wäre aber eine Möglichkeit dem nachzugehen. Anzahl der Datenpunkte ändern wolltest du ja nicht. Auch das Thema Fensterung der Daten kann eine Rolle spielen.

Ich werde das morgen aus Interesse mal an deinem Beispiel ausprobieren. Erinnere mich mal dran, falls ich es vergesse. Eine DFT kann man auch einfach von Hand selbst programmieren, macht halt niemand, da die Python-Bibliotheken sicher locker um den Faktor 1000 schneller sind.

[Sirius3]

@markus_: scipy kann zwar zero-padding anwenden, defaultmäßig werden aber keine 0en aufgefüllt. Der verwendete fft-Algorithmus kommt mit beliebigen Samplelängen zurecht

[markus_]

Gibt es eine Möglichkeit, die Spitzen zu unterscheiden, ohne die tmax-Zeit von 1,76 und die Frequenzen zu ändern?

Ich hab mich mal an deinem Beispiel versucht und dabei gemerkt, dass es doch etwas komplizierter wird. Ich glaube unter den gegebenen Randbedingungen wirst du kein zufriedenstellendes Ergebnis bekommen, wie ja auch @__deets__ schon angemerkt hat.

Die minimale Frequenzauflösung der FT ergibt sich zu df = fs/N mit fs = Samplefrequenz und N = Länge der FFT. fs ergibt sich bei dir aus deinem Zeitvektor t (fs = 1/t[1]) und N ist in deinem Fall die Länge des Zeitvektors (in deinem Beispiel 2400). Damit komme ich in deinem Beispiel auf ein df von 0,56 Hz. Das heißt der erste Frequenzpunkt deiner FFT liegt bei 0,56, der zweite dann bei 1,12 Hz usw. Deine Signale liegen genau dazwischen und somit sind die nicht mehr differenzierbar, sondern du bekommst einen breiten Peak.

Um die Auflösung zu verbessern müsstest du z.B. die Länge der FFT vergrößern bei gleichbleibender Abtastfrequenz, was aber dazu führt, dass dein Tmax größer werden wird als die vorgegebenen 1,76 sec. Somit würde ich deine Frage verneinen, du wirst eine längere Beobachtungszeit benötigen. zumindest soweit mein Verständnis von der Sache aber ich bin kein Mathematiker und schließe nicht aus, dass es da ggfs. noch irgendwelche Tricks geben könnte.

[pyhill00]

Danke für die Mühe. Am Ende des Tages muss ich dann doch 1.76 irgendwie erhöhen dann ist das Problem gelöst.