#!/usr/bin/env python
# coding: utf-8

# ## 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 = 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


# ### 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]:


get_ipython().run_line_magic('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]:


get_ipython().run_line_magic('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))


# ### 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[ ]: