Murat's Data used.
100 datasets, each has of 131072 reasilation. Mean velocity is 6.842 m/s.
import numpy as np
#import pandas as pd
from scipy.fftpack import fft, ifft,fftshift,fftfreq
data = np.loadtxt('velocity_data.out') #too slow, http://akuederle.com/stop-using-numpy-loadtxt/
#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
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;
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,
# 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
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.
#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.
%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()
%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()
'''
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