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 $ y = V (t) = \ sin (t) $. The trigonometric function with period $ can be written as $ \ sin (2 \ pi ft) $, so the period $ f = 1/2 \ pi $$ is set. The amplitude is 1.0.

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.

scipy_rawdata_beforefft.png

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

fft_norm_check.png

(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 $ amplitude = \ sqrt {2} \ sqrt {\ Sigma_i psd (f_i) ^ 2 \ Delta f} $. This value is 1 in this example

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

How to confirm the Persival theorem using the Fourier transform (FFT) of matplotlib and scipy
[Python] How to specify the window display position and size of matplotlib
How to avoid the cut-off label of the graph created by the plot module using matplotlib
FFT (Fast Fourier Transform): Formulas and Implementation Examples for Implementation
Signal processing in Python (1): Fourier transform
[Scientific / technical calculation by Python] 1D-3D discrete fast Fourier transform, scipy
Drawing with Bezier curve and Fourier transform
Fourier transform the wav file read by Python, reverse transform it, and write it again.
Image processing from scratch with python (5) Fourier transform
How to confirm the Persival theorem using the Fourier transform (FFT) of matplotlib and scipy
How to add new data (lines and plots) using matplotlib
Installation of SciPy and matplotlib (Python)
How to get followers and followers from python using the Mastodon API
[EC2] How to install chrome and the contents of each command
[Python] How to get the first and last days of the month
I summarized how to change the boot parameters of GRUB and GRUB2
Properties of the discrete Fourier transform
How to copy and paste the contents of a sheet in Google Spreadsheet in JSON format (using Google Colab)
I tried to transform the face image using sparse_image_warp of TensorFlow Addons
How to divide and process a data frame using the groupby function
I tried to extract and illustrate the stage of the story using COTOHA
How to check the version of Django
I tried to notify the update of "Hamelin" using "Beautiful Soup" and "IFTTT"
[Circuit x Python] How to find the transfer function of a circuit using Lcapy
How to count the number of elements in Django and output to a template
[Numpy, scipy] How to calculate the square root of a semi-fixed definite matrix
How to assign multiple values to the Matplotlib colorbar
How to calculate the volatility of a brand
How to find the area of the Voronoi diagram
I tried to notify the update of "Become a novelist" using "IFTTT" and "Become a novelist API"
[Blender] How to get the selection order of vertices, edges and faces of an object
[For beginners] How to display maps and search boxes using the GoogleMap Javascript API
How to find out which process is using the localhost port and stop it
Save the text of all Evernote notes to SQLite using Beautiful Soup and SQLAlchemy
How to know the port number of the xinetd service
How to write a GUI using the maya command
How to get the number of digits in Python
I tried to summarize how to use matplotlib of python
Add information to the bottom of the figure with Matplotlib
FFT (Fast Fourier Transform): Formulas and Implementation Examples for Implementation
How to visualize the decision tree model of scikit-learn
How to use the grep command and frequent samples
[Blender] How to dynamically set the selection of EnumProperty
How to use argparse and the difference between optparse
[Python] Summary of how to specify the color of the figure
How to hit the document of Magic Function (Line Magic)
How to access the global variable of the imported module
[Introduction to Python] How to stop the loop using break?
[Introduction to Python] Basic usage of the library matplotlib
The story of using circleci to build manylinux wheels
[Selenium] How to specify the relative path of chromedriver?
[Linux] [C / C ++] How to get the return address value of a function and the function name of the caller
[Super easy! ] How to display the contents of dictionaries and lists including Japanese in Python
How to know the number of GPUs from python ~ Notes on using multiprocessing with pytorch ~
How to easily draw the structure of a neural network on Google Colaboratory using "convnet-drawer"
How to plot a lot of legends by changing the color of the graph continuously with matplotlib
How to insert a specific process at the start and end of spider with scrapy
A story about porting the code of "Try and understand how Linux works" to Rust