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

# In[1]:


from pyFRF import FRF
import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')


# 
# # Showcase of the package ``pyFRF``

# Feb 2018, J. Slavič (janko.slavic@fs.uni-lj.si, [ladisk.si/~slavic](http://ladisk.si/~slavic))
# 
# pyFRF is part of the [www.openmodal.com](http://www.openmodal.com) project.
# 
# The inputs are time signals of excitation and response, the outputs are FRF estimators (H1, H2, Hv, Vector or ODS) and coherence.

# ## Create synthetic FRF

# In[2]:


C = 0.5+0.1j # modal constant
eta = 5e-3   # damping loss factor
f0 = 320     # natural frequency
df = 1      # freq resolution
D = 1e-8*(1-.1j)    # residual

f = 1*np.arange(0, 1400, step=df) # / frequency range
    
w0 = f0 * 2 * np.pi #to rad/s
w = f * 2 * np.pi
H1_syn = C / (w0 ** 2 - w ** 2 + 1.j * eta * w0 ** 2) + \
        +0.5*np.conj(C) / ((w0*2)** 2 - w ** 2 + 1.j * eta * (w0*2)** 2)\
        +0.25*C / ((w0*3)** 2 - w ** 2 + 1.j * eta * (w0*3)** 2)\
        + D


# In[3]:


fig, ax1 = plt.subplots()
ax1.semilogy(f,np.abs(H1_syn), 'b')
ax1.set_xlabel('Frequency [Hz]')
ax1.set_ylabel('H1', color='b')
for tl in ax1.get_yticklabels():
    tl.set_color('b')
    
ax2 = ax1.twinx()
ax2.plot(f,np.angle(H1_syn), 'r', alpha=0.2)
ax2.set_ylabel('Angle', color='r')
for tl in ax2.get_yticklabels():
    tl.set_color('r')


# ## Prepare synthetic impulse response

# In[12]:


h = np.fft.irfft(H1_syn)
l = len(H1_syn)*2-2
t = np.linspace(0, 1, num=l)
exc = np.zeros_like(t)
exc[0] = 1


# In[13]:


fig, ax1 = plt.subplots()
ax1.plot(t, exc, 'b');
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Excitation', color='b')
ax1.set_xlim(left=0, right=1)
for tl in ax1.get_yticklabels():
    tl.set_color('b')
    
ax2 = ax1.twinx()
ax2.plot(t, h, 'r', alpha=0.7)
ax2.set_ylabel('Response', color='r')
for tl in ax2.get_yticklabels():
    tl.set_color('r')


# ## Go back to frequency domain via pyFRF

# In[14]:


frf = FRF(sampling_freq=1/t[1], exc=exc, resp=h, exc_window='None', resp_type='d', resp_window='None')


# In[15]:


len(frf.get_f_axis())


# In[16]:


plt.semilogy(frf.get_f_axis(), np.abs(frf.get_FRF()), '.', label='via FRF')
plt.semilogy(f, np.abs(H1_syn), label='Synthetic')
plt.title('FRF H1')
plt.legend();


# In[17]:


plt.semilogy(frf.get_f_axis(), np.abs(frf.get_FRF(form='accelerance')), '.', label='Accelerance')
plt.semilogy(frf.get_f_axis(), np.abs(frf.get_FRF(form='mobility')), '.', label='Mobility')
plt.semilogy(frf.get_f_axis(), np.abs(frf.get_FRF(form='receptance')), '.', label='Receptance')
plt.title('FRF H1')
plt.legend();


# ## Go back to frequency domain via pyFRF - multiple measurements with noise

# In[18]:


# prepare the instance of FRF 
averages = 10
frf = FRF(sampling_freq=1/t[1], fft_len=len(h), exc_window='None', \
          resp_window='None', resp_type='d', weighting='Linear', n_averages=averages)


# In[19]:


k = 0.1 # rate of noise
for i in range(averages):
    noise = k * (np.random.rand(len(h))-0.5) * np.std(h)
    frf.add_data(exc, h + noise)


# In[20]:


plt.semilogy(frf.get_f_axis(), np.abs(frf.get_H1()), '.', label='via FRF')
plt.semilogy(f, np.abs(H1_syn), label='Synthetic')
plt.title('FRF H1')
plt.legend();


# In[21]:


plt.plot(frf.get_f_axis(), frf.get_coherence(), '.');


# ## Measured data, from time domain to frequency domain - multiple measurements

# Prepared by: Matej Razpotnik, matej.razpotnik@fs.uni-lj.si

# In[22]:


import glob as glob
np.seterr(all='ignore');#ignore numpy warning due to division with 0 at frequency derivation


# In[23]:


measurement_files = glob.glob('data/*')
first_meas = np.load(measurement_files[0])


# In[24]:


fig, ax1 = plt.subplots()
ax1.plot(first_meas[0], first_meas[1], 'b');
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Excitation', color='b')
ax1.set_xlim(left=0, right=0.002)
for tl in ax1.get_yticklabels():
    tl.set_color('b')
    
ax2 = ax1.twinx()
ax2.plot(first_meas[0], first_meas[2], 'r', alpha=0.5)
ax2.set_ylabel('Response', color='r')
for tl in ax2.get_yticklabels():
    tl.set_color('r')
ax1.set_xlim(0, 0.01);


# In[25]:


frf_meas = FRF(sampling_freq=1 / (first_meas[0][1] - first_meas[0][0]), 
               fft_len=len(first_meas[0]), 
               exc_window='Force:0.01', 
               resp_window='Exponential:0.01',
               weighting='Linear', 
               n_averages=len(measurement_files))


# In[26]:


plot_up_to = 4000
for file_name in measurement_files:
    meas_i = np.load(file_name)
    frf_meas.add_data(meas_i[1], meas_i[2])
plt.semilogy(frf_meas.get_f_axis()[:plot_up_to], np.abs(frf_meas.get_FRF(form='accelerance'))[:plot_up_to])
plt.title('FRF H1');


# In[27]:


plt.plot(frf_meas.get_f_axis()[:plot_up_to], frf_meas.get_coherence()[:plot_up_to]);
plt.title('Coherence');