## Example: FFT Analysis of Hot Wire data via Python¶

Murat's Data used.

100 datasets, each has of 131072 reasilation. Mean velocity is 6.842 m/s.

In [5]:
import numpy as np
#import pandas as pd
from scipy.fftpack import fft, ifft,fftshift,fftfreq

#data = pd.read_csv('../../lfs/murats_velocity_data.out') #will not work that with bare version
print("Data loaded on Ipython environment, and it's shape is", data.shape)

mean_data = np.mean(data)
print("Mean velocity is", mean_data)

fluc_data= data - mean_data

Data loaded on Ipython environment, and it's shape is (131072, 100)
Mean velocity is 6.84244558364


### Function to calculate spectrum¶

Often magnitude of FT or FFT is not important, for us it is!! So normalization is important, it is not just the shape.

There are 3 steps to follow;

1. Take the FFT of the signal
2. Calculate spectrum
3. Normalization

Without understanding all these 3, it is difficult get correct transform.

It is always possible to check that if your transform is correct or not. To do so,

1. Check the integral of the are under the estimation of S(f), that should be equalt to var(signal)
2. IFFT of spectrum should give the correlation, and correlation can be also calculated with brute force, so two correlation can be compared (tidious FFT has windowing effect always, it needs taken out/handled properly)
3. In case having both time and space resolved data, spectrums of these to can be compared, of course within the limits of Taylor Hypothesis. Moreover, energy levels should reasonably match, it is a good way to double check correctness of wavenumber/wavelenght conversations.
In [9]:
# Calculates the spectrum of a given signal with fft
def calc_spectrum(mysignal, size, fs, typeWindow):
# Window
# Add window or periodization here
if typeWindow == 'Hanning':
mysignal    = mysignal[0:size]*np.hanning(size)
integWindow = np.trapz(np.hanning(size),np.linspace(0,1,size))
elif typeWindow == 'None':
mysignal    = mysignal[0:size]*1.0
integWindow = 1.0

# Take the FFT of the signal, -1 important
# Remember you have the -1 here if you periodize your signal outside of this routine
# Multiplying by 2 or not, and having -fs/2:fs/2 frequencies should be compatible
# following FFT double sided no need to multiply by 2
# Skip the division by N-1 as well, for the moment, normalization is done at the end
myfft  = fft(mysignal[0:size-1]) #*2.0 #/float(size-1)

# Calculate spectrum
# (real**2.0+imag**2.0) is kind of variance
# if you take sqrt() that makes it kind of standart deviation, same unit as velocity,
# nice for Fluid Mechanics applications
# But let's keep (standart deviation)**2.0, so no need take sqrt()
myspec = (myfft.real[0:int((size-1)/2)]**2.0+myfft.imag[0:int((size-1)/2)]**2.0) #**0.5
# OR :: Common way, attention it is missing the sqrt() as well, good
#myspec = myfft[0:(size)/2]*np.conj(myfft)[0:(size)/2]

# Normalization
# EXAMPLE :: if the data collected with 30kHz, and you want to have first 256 points, at every 2 points
# First 256 points that gives the new size
# New fs 30kHz/2 = 15kHz
myspec = myspec/(float(size-1)*fs)/(integWindow)

return myspec

In [11]:
N=256              #Choose record length, maximum is 131072
#N=2000

fs=30000                                                 # Define fs
freq_3 = fftfreq(N, d=1./fs)[0:int((N/2)-1)]                  # Create frequencies, needed to plot, pick the half of freq.
spec_3 = calc_spectrum(fluc_data[:,0], N, fs, 'None')         # Call function

fs=15000                                                 # Define fs
freq_4 = fftfreq(int(N/2), d=1./fs)[0:int((N/2)/2)-1]              # Create frequencies, needed to plot, pick the half of freq.
spec_4 = calc_spectrum(fluc_data[0::2,0], int(N/2), fs,'None')       # Call function

# Loop over the packages (note: this is unefficient)
for i in range(1,100):
fs=30000
spec_3 += calc_spectrum(fluc_data[:   ,i], N, fs,'None')
fs=15000
spec_4 += calc_spectrum(fluc_data[0::2,i], int(N/2), fs,'None')

# Do not forget the averaging packages
spec_3 = spec_3/100.
spec_4 = spec_4/100.

In [13]:
#N=256              #Choose record length, maximum is 131072
N=13072

fs=30000                                                 # Define fs
freq_1 = fftfreq(N, d=1./fs)[0:int((N/2)-1)]                  # Create frequencies, needed to plot, pick the half of freq.
spec_1 = calc_spectrum(fluc_data[:   ,0], N, fs, 'None')         # Call function

# Loop over the packages (note: this is unefficient)
for i in range(1,100):
fs=30000
spec_1 += calc_spectrum(fluc_data[:   ,i], N, fs,'None')

# Do not forget the averaging packages
spec_1 = spec_1/100.

In [14]:
%matplotlib notebook
import matplotlib.pyplot as plt

fig = plt.figure()
ax = plt.subplot(111)

#Put data into the plot
ax.loglog(freq_1,spec_1, "k*-")
ax.loglog(freq_3,spec_3, "r.-")
ax.loglog(freq_4,spec_4, "g-")
# Shrink current axis's height by 10% on the bottom
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1,
box.width, box.height * 0.9])

plt.grid()
plt.title("N=256")
plt.xlabel("Frequency")
plt.ylabel("Spectrum")
plt.show()

In [15]:
%matplotlib notebook
import matplotlib.pyplot as plt

fig = plt.figure()
ax = plt.subplot(111)

#Put data into the plot
ax.semilogy(freq_1,spec_1, "k*-")
ax.semilogy(freq_3,spec_3, "r.-")
ax.semilogy(freq_4,spec_4, "g-")
# Shrink current axis's height by 10% on the bottom
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1,
box.width, box.height * 0.9])

plt.grid()
plt.title("N=256")
plt.xlabel("Frequency")
plt.ylabel("Spectrum")
plt.show()

fig = plt.figure()
ax = plt.subplot(111)

#Put data into the plot
ax.semilogx(freq_1,spec_1, "k*-")
ax.semilogx(freq_3,spec_3, "r.-")
ax.semilogx(freq_4,spec_4, "g-")
# Shrink current axis's height by 10% on the bottom
box = ax.get_position()
ax.set_position([box.x0, box.y0 + box.height * 0.1,
box.width, box.height * 0.9])

plt.grid()
plt.title("N=256")
plt.xlabel("Frequency")
plt.ylabel("Spectrum")
plt.show()

In [16]:
'''
Let's check the integral of the are under the estimation of S(f)
'''

# The way calculated it to spectrum is similar to the what is given in Bill's books
# So we need multiply integral by 2. (we took the half of S(f) only)
I1 = np.trapz(spec_1,freq_1)*2.0
I3 = np.trapz(spec_3,freq_3)*2.0
I4 = np.trapz(spec_4,freq_4)*2.0
print('integral of S(f) 1', I1)
print('integral of S(f) 3', I3)
print('integral of S(f) 4', I4)
print('variance of u', np.var(fluc_data))

integral of S(f) 1 0.604697071356
integral of S(f) 3 0.617515167469
integral of S(f) 4 0.615620322208
variance of u 0.593587975241


### Try yourself¶

• Please, test with different windows, any satisfactory result?
• Normalization issues are solved, try to add into this document correlation calculations as well, and then compare with IFFT(S(f)) with brute force correlation results.
In [ ]: