[PYTHON] Comment confirmer le théorème de Persival en utilisant matplotlib et la transformée de Fourier de Scipy (FFT)
Comment calculer le spectre de puissance avec python
Comment calculer la FFT avec matplotlib et scipy
Matplotlib et scipy ont tous deux des méthodes pour calculer le spectre de puissance, mais ils ont des définitions différentes des fonctions par défaut. La plupart de l'utilisation est la suivante.
--Si vous voulez simplement calculer le spectre de puissance, mlab.psd est le plus simple. -Utilisez scipy.fftpack pour obtenir les parties réelles et imaginaires. Il est possible de faire des choses détaillées telles que le calcul de phase, ou vous devez vous normaliser.
Quoi qu'il en soit, la puissance totale sera la même que celle du [théorème de Perseval](https://ja.wikipedia.org/wiki/ théorème de Persebal), mais si vous ne considérez pas la différence de standardisation entre les deux, le spectre de puissance L'axe vertical n'est pas exactement le même.
Voici un exemple où une simple fonction triangulaire est transformée de Fourier par les deux méthodes et le résultat est exactement le même dans matplotlib et scipy FFT. En prime, j'expliquerai comment trouver l'amplitude d'une onde triangulaire à partir du spectre de puissance et la rendre identique à l'amplitude de la fonction triangulaire d'origine, et quel type d'effet elle produit lorsqu'un filtre est ajouté.
Problème de réglage
Fonction avant transformée de Fourier
L'axe vertical n'a pas de signification particulière, mais considérons une seule fonction triangulaire avec la tension de signal V (t). L'expression la plus simple est
Il s'agit de mlab.psd et scipy.fftpack /generated/scipy.fftpack.fft.html) et transformée de Fourier pour calculer la puissance. Calculez chaque cas avec filtre suspendu et voyez la différence.
La résolution temporelle dt est de 0,1 seconde. [Fréquence de Nyquist](https://ja.wikipedia.org/wiki/Fréquence de Nyquist) f_n est la moitié de l'inverse de celle, qui est de 5 Hz.
Exemple de code
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()
La méthode d'exécution consiste simplement à ne rien spécifier ou à modifier la largeur de temps avec -t.
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
Données de séries temporelles avant transformée de Fourier
Lorsqu'elle est exécutée sans spécifier d'options, la transformation de Fourier suivante des données de séries temporelles est exécutée.
Lorsque le filtre suspendu est appliqué, le filtre suspendu (bleu) est appliqué aux données de série chronologique d'origine (rouge) et les données vertes sont transformées de Fourier. Si aucun filtre n'est utilisé, le rouge sera transformé de Fourier sous la condition aux limites périodique. Lors de l'application d'un filtre, s'il y a un décalage (fréquence = 0 composante), un bruit artificiel basse fréquence peut être ajouté par le filtre, donc lors de l'utilisation d'un filtre, il est préférable de soustraire la composante CC avant de l'utiliser. Bien.
Résultat du calcul du spectre de puissance
(Supérieur) Données chronologiques à transformer de Fourier. Les amplitudes des données chronologiques calculées à partir du spectre de puissance ont été superposées et affichées. (Milieu) Le rouge est lorsqu'il n'y a pas de filtre et le bleu lorsqu'il y a un filtre suspendu. La ligne pleine est calculée par matplotlib, et les points sont calculés par scipy. La fine ligne rouge est la valeur de puissance estimée à partir des données de la série chronologique. (Inférieur) Version linéaire à axe vertical au milieu.
Confirmation du [Théorème de Perséval](https://ja.wikipedia.org/wiki/ Théorème de Persebal)
D'après le [Théorème de Perséval](https://ja.wikipedia.org/wiki/Perseval's Theorem), les puissances des valeurs intégrées dans l'espace de fréquence et l'espace et le temps correspondent. Dans le cas d'une fonction triangulaire, la puissance de fluctuation (RMS) et l'amplitude sont reliées par $ \ sqrt {2} $, donc l'amplitude et la RMS peuvent être converties l'une à l'autre. Ici, l'amplitude des données de série temporelle (1 dans cet exemple) a été calculée à partir de la puissance obtenue en intégrant complètement l'espace de fréquence, et la correspondance est indiquée par la ligne en pointillé bleue dans l'image ci-dessus.
La partie qui s'intègre dans l'espace des fréquences est le code
python
amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp = " + str(amp))
Correspond à. Il est calculé par
python
plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)
Il est tracé avec une ligne pointillée bleue dans la rangée supérieure.
Effet du filtre suspendu
Sans filtres, mlab.psd et [scipy.fftpack](https://docs.scipy.org/doc/ scipy / reference / generated / scipy.fftpack.fft.html) est une correspondance parfaite. Si vous utilisez hanning filter, la puissance diminuera à un certain rythme, vous devez donc corriger ce montant, mais la définition est-elle différente entre matplotlib et scipy? (Non fait), si vous utilisez la valeur par défaut, les résultats ne correspondront pas.
Supplément de code
Comment utiliser 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)
Ici, c est les données acquises en même temps que la largeur de détection que vous souhaitez FFT. La FFT suppose que les échantillons ont été prélevés exactement au même intervalle de temps. (Si l'intervalle de temps change, la transformée de Fourier discrète est effectuée directement, mais le calcul prend beaucoup de temps. Dans certains cas, comme dans une expérience précise.) FNum est le nombre d'échantillons et 1 / timebin est l'inverse de la largeur de temps. Si vous spécifiez None pour window, rien n'est fait, et il y a différents filtres dans les options de mlab, mais voici un exemple de hanning. Côtés signifie qu'un côté est normalisé de sorte que la puissance totale corresponde à la puissance de la série temporelle lorsqu'elle est intégrée d'un côté (valeur positive de la fréquence). Utilisez scale_by_freq pour le corriger par fréquence. En conséquence, la valeur de retour devient une unité de puissance ^ 2 / Hz. Ici, à la fin, psd = np.sqrt (psd2) est défini sur l'unité de puissance / sqrt (Hz).
Comment utiliser 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)]
Il calcule la fréquence en utilisant une fonction appelée fftfreq. C'est en fait une caractéristique très importante. La raison en est qu'ici dans scipy fft, une transformation de Fourier complexe est effectuée, mais la valeur de retour est le coefficient de la partie réelle et de la partie imaginaire, et il n'est pas évident dans quel ordre de fréquence il est renvoyé à l'utilisateur. Nous avons besoin d'une fonction qui renvoie les fréquences correspondantes.
python
real = fft.real
imag = fft.imag
psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))
Ensuite, sortez la partie imaginaire de la partie réelle et standardisez-la d'un côté (la fréquence est positive) pour qu'il y ait une queue de livre.
Recommended Posts