[PYTHON] Wie man den Satz von Persival mit Matplotlib und der Fourier-Transformation (FFT) von scipy bestätigt
So berechnen Sie das Leistungsspektrum mit Python
Wie berechnet man die FFT mit Matplotlib und Scipy?
Sowohl matplotlib als auch scipy haben Methoden zur Berechnung des Leistungsspektrums, aber sie haben andere Definitionen als die Standardfunktionen. Der Großteil der Verwendung ist wie folgt.
- Wenn Sie einfach das Leistungsspektrum berechnen möchten, ist mlab.psd am einfachsten.
- Verwenden Sie scipy.fftpack, um die Real- und Imaginärteile abzurufen. Es ist möglich, detaillierte Dinge wie die Phasenberechnung durchzuführen, oder Sie müssen sich selbst standardisieren.
Was auch immer Sie tun, die Gesamtleistung ist dieselbe wie in [Persevals Theorem](https://ja.wikipedia.org/wiki/Persebals Theorem), aber wenn Sie den Unterschied in der Standardisierung zwischen beiden nicht berücksichtigen, das Leistungsspektrum Die vertikale Achse ist nicht genau gleich.
Hier ist ein Beispiel, bei dem eine einfache Dreiecksfunktion mit beiden Methoden Fourier-transformiert wird und das Ergebnis in matplotlib und scipy FFT genau das gleiche ist. Als Bonus zeige ich Ihnen, wie Sie die Amplitude einer Dreieckswelle aus dem Leistungsspektrum ermitteln, sie der Amplitude der ursprünglichen Dreiecksfunktion anpassen und welchen Effekt sie hat, wenn Sie einen Filter hinzufügen.
Problemstellung
Funktion vor Fourier-Transformation
Die vertikale Achse hat keine besondere Bedeutung, betrachtet jedoch eine einzelne Dreiecksfunktion mit der Signalspannung V (t). Die einfachste Formel lautet
Dies sind mlab.psd und scipy.fftpack /generated/scipy.fftpack.fft.html) und Fourier-Transformation zur Berechnung der Leistung. Berechnen Sie jeden Fall mit hanning filter und sehen Sie den Unterschied.
Die Zeitauflösung dt beträgt 0,1 Sekunden. [Nyquist-Frequenz](https://ja.wikipedia.org/wiki/Nyquist Frequenz) f_n ist die Hälfte der Umkehrung von 5 Hz.
Beispielcode
qiita_fft_matplotlib_scipy_comp.py
#!/usr/bin/env python
"""
#2013-03-12 ; ver 1.0; First version, using Sawada-kun's code
#2013-08-28 ; ver 1.1; refactoring
#2020-03-31 ; ver 1.2; updated for python3
"""
__version__= '1.2'
import numpy as np
import scipy as sp
import scipy.fftpack as sf
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import optparse
def scipy_fft(inputarray,dt,filtname=None, pltfrag=False):
bin = len(inputarray)
if filtname == None:
filt = np.ones(len(inputarray))
elif filtname == 'hanning':
# filt = sp.hanning(len(inputarray))
filt = np.hanning(len(inputarray))
else:
print('No support for filtname=%s' % filtname)
exit()
freq = sf.fftfreq(bin,dt)[0:int(bin/2)]
df = freq[1]-freq[0]
fft = sf.fft(inputarray*filt,bin)[0:int(bin/2)]
real = fft.real
imag = fft.imag
psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))
if pltfrag:
binarray = range(0,bin,1)
F = plt.figure(figsize=(10,8.))
plt.subplots_adjust(wspace = 0.1, hspace = 0.3, top=0.9, bottom = 0.08, right=0.92, left = 0.1)
ax = plt.subplot(1,1,1)
plt.title("check fft of scipy")
plt.xscale('linear')
plt.grid(True)
plt.xlabel(r'number of bins')
plt.ylabel('FFT input')
plt.errorbar(binarray, inputarray, fmt='r', label="raw data")
plt.errorbar(binarray, filt, fmt='b--', label="filter")
plt.errorbar(binarray, inputarray * filt, fmt='g', label="filtered raw data")
plt.legend(numpoints=1, frameon=False, loc=1)
plt.savefig("scipy_rawdata_beforefft.png ")
plt.show()
return(freq,real,imag,psd)
usage = u'%prog [-t 100] [-d TRUE]'
version = __version__
parser = optparse.OptionParser(usage=usage, version=version)
parser.add_option('-f', '--outputfilename', action='store',type='string',help='output name',metavar='OUTPUTFILENAME', default='normcheck')
parser.add_option('-d', '--debug', action='store_true', help='The flag to show detailed information', metavar='DEBUG', default=False)
parser.add_option('-t', '--timelength', action='store', type='float', help='Time length of input data', metavar='TIMELENGTH', default=20.)
options, args = parser.parse_args()
argc = len(args)
outputfilename = options.outputfilename
debug = options.debug
timelength = options.timelength
timebin = 0.1 # timing resolution
# Define Sine curve
inputf = 1/ (2 * np.pi) # f = 1/2pi -> sin (t)
inputt = 1 / inputf # T = 1/f
t = np.arange(0.0, timelength*np.pi, timebin)
c = np.sin(t)
fNum = len(c)
Nyquist = 0.5 / timebin
Fbin = Nyquist / (0.5 * fNum)
fftnorm = 1. / (np.sqrt(2) * np.sqrt(Fbin))
print("................. Time domain ................................")
print("timebin = " + str(timebin) + " (sec) ")
print("frequency = " + str(inputf) + " (Hz) ")
print("period = " + str(inputt) + " (s) ")
print("Nyquist (Hz) = " + str(Nyquist) + " (Hz) ")
print("Frequency bin size (Hz) = " + str(Fbin) + " (Hz) ")
print("FFT norm = " + str(fftnorm) + "\n")
# Do FFT.
# Hanning Filtered
psd2, freqlist = mlab.psd(c,
fNum,
1./timebin,
window=mlab.window_hanning,
sides='onesided',
scale_by_freq=True
)
psd = np.sqrt(psd2)
spfreq,spreal,spimag,sppower = scipy_fft(c,timebin,'hanning',True)
# No Filter
psd2_nofil, freqlist_nofil = mlab.psd(c,
fNum,
1./timebin,
window=mlab.window_none,
sides='onesided',
scale_by_freq=True
)
psd_nofil = np.sqrt(psd2_nofil)
spfreq_nf,spreal_nf,spimag_nf,sppower_nf = scipy_fft(c,timebin,None)
# Get input norm from FFT intergral
print("................. FFT results ................................")
amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp = " + str(amp))
peakval = amp / (np.sqrt(2) * np.sqrt(Fbin) )
print("Peakval = " + str(peakval))
# (1) Plot Time vs. Raw value
F = plt.figure(figsize=(10,8.))
plt.subplots_adjust(wspace = 0.1, hspace = 0.3, top=0.9, bottom = 0.08, right=0.92, left = 0.1)
ax = plt.subplot(3,1,1)
plt.title("FFT Norm Check")
plt.xscale('linear')
plt.xlabel(r'Time (s)')
plt.ylabel('FFT input (V)')
plt.errorbar(t, c, fmt='r', label="raw data")
plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)
plt.ylim(-2,2)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)
# (2a) Plot Freq vs. Power (log-log)
ax = plt.subplot(3,1,2)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'Frequency (Hz)')
plt.ylabel(r'FFT Output (V/$\sqrt{Hz}$)')
plt.errorbar(freqlist_nofil , psd_nofil, fmt='r-', label="Window None, matplotlib", alpha=0.8)
plt.errorbar(freqlist, psd, fmt='b-', label="Window Hanning, matplotlib", alpha=0.8)
plt.errorbar(spfreq_nf, sppower_nf, fmt='mo', label="Window None, scipy", ms = 3, alpha=0.8)
plt.errorbar(spfreq, sppower, fmt='co', label="Window Hanning, scipy", ms = 3, alpha=0.8)
plt.errorbar(freqlist, fftnorm * np.ones(len(freqlist)), fmt='r--', label=r"1./($\sqrt{2} \sqrt{Fbin}$)", alpha=0.5, lw=1)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)
# (2b) Plot Freq vs. Power (lin-lin)
ax = plt.subplot(3,1,3)
plt.xscale('linear')
plt.xlim(0.1,0.2)
plt.yscale('linear')
plt.xlabel(r'Frequency (Hz)')
plt.ylabel(r'FFT Output (V/$\sqrt{Hz}$)')
plt.errorbar(freqlist_nofil , psd_nofil, fmt='r-', label="Window None, matplotlib", alpha=0.8)
plt.errorbar(freqlist, psd, fmt='b-', label="Window Hanning, matplotlib", alpha=0.8)
plt.errorbar(spfreq_nf, sppower_nf, fmt='mo', label="Window None, scipy", ms = 3, alpha=0.8)
plt.errorbar(spfreq, sppower, fmt='co', label="Window Hanning, scipy", ms = 3, alpha=0.8)
plt.errorbar(freqlist, fftnorm * np.ones(len(freqlist)), fmt='r--', label=r"1./($\sqrt{2} \sqrt{Fbin}$)", alpha=0.5, lw=1)
plt.legend(numpoints=1, frameon=True, loc="best", fontsize=8)
plt.savefig("fft_norm_check.png ")
plt.show()
Die Ausführungsmethode besteht einfach darin, nichts anzugeben oder die Zeitbreite mit -t zu ändern.
python
$ python qiita_fft_matplotlib_scipy_comp.py -h
Usage: qiita_fft_matplotlib_scipy_comp.py [-t 100] [-d TRUE]
Options:
--version show program's version number and exit
-h, --help show this help message and exit
-f OUTPUTFILENAME, --outputfilename=OUTPUTFILENAME
output name
-d, --debug The flag to show detailed information
-t TIMELENGTH, --timelength=TIMELENGTH
Time length of input data
Zeitreihendaten vor der Fourier-Transformation
Bei Ausführung ohne Angabe von Optionen wird die folgende Fourier-Transformation von Zeitreihendaten ausgeführt.
Wenn der Hanning-Filter angewendet wird, wird der Hanning-Filter (blau) auf die ursprünglichen Zeitreihendaten (rot) angewendet, und die grünen Daten werden Fourier-transformiert. Wenn kein Filter verwendet wird, wird Rot unter der periodischen Randbedingung Fourier-transformiert. Wenn beim Anwenden eines Filters ein Offset (Frequenz = 0-Komponente) vorhanden ist, kann das Filter künstliches Niederfrequenzrauschen hinzufügen. Bei Verwendung eines Filters ist es daher besser, die DC-Komponente vor der Verwendung zu subtrahieren. Gut.
Ergebnis der Leistungsspektrumberechnung
(Obere) Zeitreihendaten, die Fourier-transformiert werden sollen. Die aus dem Leistungsspektrum berechneten Amplituden der Zeitreihendaten wurden überlagert und angezeigt. (Mitte) Rot ist, wenn es keinen Filter gibt, und Blau ist, wenn es einen Hanning-Filter gibt. Die durchgezogene Linie wird mit matplotlib berechnet, und die Punkte werden mit scipy berechnet. Die dünne rote Linie ist der aus den Zeitreihendaten geschätzte Leistungswert. (Untere) Lineare Version der vertikalen Achse in der Mitte.
Bestätigung des [Satzes von Perseval](https://ja.wikipedia.org/wiki/Persebals Satz)
Aus [Persevals Theorem](https://ja.wikipedia.org/wiki/Persevals Theorem) stimmen die Potenzen der integrierten Werte in Frequenzraum und Raum und Zeit überein. Bei einer Dreiecksfunktion sind die Schwankungsleistung (RMS) und die Amplitude durch $ \ sqrt {2} $ verbunden, sodass Amplitude und RMS ineinander umgewandelt werden können. Hier wurde die Amplitude der Zeitreihendaten (1 in diesem Beispiel) aus der Leistung berechnet, die durch vollständige Integration des Frequenzraums erhalten wurde, und die Übereinstimmung wird durch die blau gepunktete Linie im obigen Bild gezeigt.
Der Teil, der in den Frequenzraum integriert wird, ist der Code
python
amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp = " + str(amp))
Entspricht. Sie wird berechnet durch
python
plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)
Es ist mit einer blau gepunkteten Linie in der oberen Reihe dargestellt.
Wirkung des Hanning-Filters
Ohne Filter mlab.psd und [scipy.fftpack](https://docs.scipy.org/doc/ scipy / reference / generate / scipy.fftpack.fft.html) passt perfekt zusammen. Wenn Sie hanning filter verwenden, nimmt die Leistung mit einer bestimmten Geschwindigkeit ab, sodass Sie diesen Betrag korrigieren müssen. Unterscheidet sich die Definition jedoch zwischen matplotlib und scipy? (Nicht erledigt) Wenn Sie die Standardeinstellung verwenden, stimmen die Ergebnisse nicht überein.
Code-Ergänzung
Verwendung von mlab.psd
python
# Do FFT.
# Hanning Filtered
psd2, freqlist = mlab.psd(c,
fNum,
1./timebin,
window=mlab.window_hanning,
sides='onesided',
scale_by_freq=True
)
psd = np.sqrt(psd2)
Hier ist c die Daten, die zur gleichen Zeit erfasst werden, wie Sie die Breite erfassen möchten. Die FFT geht davon aus, dass die Proben genau im gleichen Zeitintervall entnommen wurden. (Wenn sich das Zeitintervall ändert, wird die diskrete Fourier-Transformation direkt durchgeführt, die Berechnung dauert jedoch lange. In einigen Fällen, z. B. in einem präzisen Experiment.) FNum ist die Anzahl der Abtastwerte und 1 / timebin ist die Umkehrung der Zeitbreite. Wenn Sie für das Fenster Keine angeben, wird nichts unternommen, und die Optionen von mlab enthalten verschiedene Filter. Hier ist jedoch ein Beispiel für das Hanning. Seiten bedeutet, dass einseitig so standardisiert ist, dass die Gesamtleistung bei einseitiger Integration mit der Leistung in der Zeitreihe übereinstimmt (positiver Frequenzwert). Verwenden Sie scale_by_freq, um es pro Frequenz festzulegen. Infolgedessen wird der Rückgabewert zu einer Leistungseinheit ^ 2 / Hz. Hier wird am Ende psd = np.sqrt (psd2) auf die Einheit Leistung / sqrt (Hz) gesetzt.
Verwendung von scipy.fftpack
python
freq = sf.fftfreq(bin,dt)[0:int(bin/2)]
df = freq[1]-freq[0]
fft = sf.fft(inputarray*filt,bin)[0:int(bin/2)]
Es berechnet die Frequenz mithilfe einer Funktion namens fftfreq. Dies ist eigentlich ein sehr wichtiges Merkmal. Der Grund ist, dass hier in scipys fft eine komplexe Fourier-Transformation durchgeführt wird, der Rückgabewert jedoch der Koeffizient des Realteils und des Imaginärteils ist und es nicht offensichtlich ist, in welcher Frequenzreihenfolge er an den Benutzer zurückgegeben wird. Wir brauchen eine Funktion, die die entsprechenden Frequenzen zurückgibt.
python
real = fft.real
imag = fft.imag
psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))
Nehmen Sie dann den Imaginärteil des Realteils heraus und standardisieren Sie ihn auf einer Seite (Frequenz ist positiv), sodass ein Buchschwanz entsteht.
Recommended Posts