[PYTHON] How to confirm the Persival theorem using the Fourier transform (FFT) of matplotlib and scipy
How to calculate the power spectrum in python
How to calculate FFT with matplotlib and scipy
Both matplotlib and scipy have methods to calculate the power spectrum, but they have different definitions from the default functions. Most of the usage is as follows.
--If you simply want to calculate the power spectrum, mlab.psd is the easiest. -Use scipy.fftpack to get the real and imaginary parts. It is possible to do detailed calculations such as phase calculation, or you need to standardize yourself.
In either case, the total power will be the same from [Parseval's theorem](https://ja.wikipedia.org/wiki/Parseval's theorem), but the power spectrum will be The vertical axis is not exactly the same.
Here is an example where a simple trigonometric function is Fourier transformed by both methods and the result is exactly the same in matplotlib and scipy FFT. As a bonus, I will explain how to find the amplitude of a triangular wave from the power spectrum, make it the same as the amplitude of the original trigonometric function, and what effect it has when a filter is added.
Problem setting
Function before Fourier transform
The vertical axis has no particular meaning, but consider a single trigonometric function with a signal of voltage V (t). The simplest formula is
These are mlab.psd and scipy.fftpack /generated/scipy.fftpack.fft.html) and Fourier transform to calculate the power. Calculate each case with hanning filter and see the difference.
The time resolution dt is 0.1 seconds. [Nyquist frequency](https://ja.wikipedia.org/wiki/Nyquist frequency) f_n is half the reciprocal of that, which is 5Hz.
Sample 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()
The execution method is simply to specify nothing or change the time width with -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
Time series data before Fourier transform
When executed without specifying any options, the following Fourier transform of time series data is executed.
When the hanning filter is applied, the hanning filter (blue) is applied to the original time series data (red), and the green data is Fourier transformed. If no filter is used, the red will be Fourier transformed under the periodic boundary condition. When applying a filter, if there is an offset (frequency = 0 component), artificial low frequency noise may be added by the filter, so when using a filter, it is better to subtract the DC component before using it. Good.
Power spectrum calculation result
(Upper) Time series data to be Fourier transformed. The amplitudes of the time series data calculated from the power spectrum are superimposed and displayed. (Middle) Red is when there is no filter, and blue is when there is a hanning filter. The solid line is calculated by matplotlib, and the dots are calculated by scipy. The thin red line is the power value estimated from the time series data. (Lower) Vertical axis linear version in the middle.
Confirmation of [Parseval's Theorem](https://ja.wikipedia.org/wiki/Parseval's Theorem)
From [Parseval's theorem](https://ja.wikipedia.org/wiki/Parseval's theorem), the powers of the integral values of frequency space and space-time match. In the case of trigonometric functions, the power of fluctuation (RMS) and the amplitude are connected by $ \ sqrt {2} $, so the amplitude and RMS can be converted to each other. Here, the amplitude of the time series data (1 in this example) was calculated from the power obtained by fully integrating the frequency space, and the agreement is shown by the blue dotted line in the above picture.
The part of the code that integrates in the frequency space is
python
amp = np.sqrt(2) * np.sqrt(np.sum(psd * psd) * Fbin)
print("Amp = " + str(amp))
Corresponds to. It is calculated by
python
plt.errorbar(t, amp * np.ones(len(c)), fmt='b--', label="Amplitude from FFT", alpha=0.5, lw=1)
It is plotted with a blue dotted line in the upper row.
Effect of hanning filter
Without filters, mlab.psd and [scipy.fftpack](https://docs.scipy.org/doc/ scipy / reference / generated / scipy.fftpack.fft.html) is a perfect match. If you use hanning filter, the power will decrease at a certain rate, so you need to correct that amount, but is the definition different between matplotlib and scipy? (Not done), if you use the default, the results will not match.
Code supplement
How to use 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)
Here, c is the data acquired with the same time sense width that you want to FFT. The FFT assumes that the samples were taken at exactly the same time interval. (If the time interval changes, the discrete Fourier transform is performed directly, but it takes a long time to calculate. In some cases, such as in a precise experiment.) FNum is the number of samples, and 1 / timebin is the inverse of the time width. If you specify None for window, nothing is done, and there are various filters in the mlab options, but here is an example of hanning. Sides is onesided and means that the total power is normalized so that it matches the power in the time series when integrated on one side (the frequency is a positive value). Use scale_by_freq to fix it per frequency. As a result, the return value becomes the unit of power ^ 2 / Hz. Here, at the end, psd = np.sqrt (psd2) is set to the unit of power / sqrt (Hz).
How to use 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)]
It calculates the frequency by using a function called fftfreq. This is actually a very important feature. The reason is that here in scipy fft, the complex Fourier transform is executed, but the return value is the coefficient of the real part and the imaginary part, and it is not obvious in what frequency order it is returned to the user. We need a function that returns the corresponding frequencies.
python
real = fft.real
imag = fft.imag
psd = np.abs(fft)/np.sqrt(df)/(bin/np.sqrt(2.))
Then, take out the imaginary part of the real part and standardize it on one side (frequency is positive) so that there is a book tail.
Recommended Posts