=== Step 6: Make FIR filters and save into a npz file ===
 * We assume a simple DARM model with
   * Actuation (A) with f^{-2} (fc=0.7, Q=100) and delay of 61 us x2
   * Sensing (C) with frequency independent and delay of 61 us
 * By importing [[https://docs.scipy.org/doc/scipy/reference/signal.html|scipy.signal]] we can construct and manipulate LTI models easily. 
{{{#!python
import scipy.signal as sg
import numpy as np
import matplotlib.pyplot as plt
}}}
 * Here is a small piece of python code to represent a f^{-2} pole 
{{{#!python
# Convert f0 and Q to a pair of complex numbers                                 
def fq2c(f0, Q, rad=False):
    w0 = -f0 if rad else -2*np.pi*f0
    a  = w0/(2*Q)
    b  = w0/(2*Q)*np.sqrt(4*Q**2-1)
    return np.array([a+1j*b, a-1j*b])

# Complex-conjugate pair of poles at f with Q and the DC gain = gn              
def pole2(f, Q, gn=1, dt=0, rad=False):
    (a,b) = fq2c(f, Q, rad);
    w1 = 2j*np.pi
    g0 = abs(1/(w1-a)/(w1-b))
    return sg.lti([], [a,b], gn/g0, dt=dt)
}}}
 * Here is the part to define DARM model parameters
{{{#!python
fNyq = 8192           # Nyquist frequency                                       

# Actuation model                                                               
f0   = 0.7            # pole frequency                                          
Q    = 100            # pole Q                                                  
gn   = 0.917e-10      # actuation gain (m/ct)                                   
tda  = 2*(1./2)**14   # actuation delay (2 digits)                              
Tla  = 2.             # Length of time-domain filter in second                  
fca  = 4.             # High-pass filter frequency                              
dla  = Tla*fNyq
act  = pole2(f0, Q, gn)

# Inverse sensing model                                                         
C    = 9.356e-12      # inverse optical gain (m/ct)                             
tdc  = (1./2)**14     # sensing delay (1 digit)                               
Tlc  = 0.5            # Length of time-domain filter in second                  
fcc  = 8.             # High-pass filter frequency                              
dlc  = Tlc*fNyq
sen  = sg.lti([], [], C)
}}}
 * Here is the part to generate the FIR filter and save to npz (numpy zipped) file
{{{#!python
# Generate time-domain filters                                                  
fla = gentd.gentd(act, tda, fca, Tla, fNyq)
flc = gentd.gentd(sen, tdc, fcc, Tlc, fNyq)

np.savez('filt.npz', a=np.array(fla), c=np.array(flc), dla=dla, dlc=dlc)
}}}
 * Another python script gentd.py contains the actual part to generate time-domain FIR filter
 * Here is a core part of gentd.py. It is just an inverse-FFT from frequency response data into time domain
 * Please note that high-pass filter is applied at the frequency specified by fcut. This is needed to minimize the FIR filter length to get reasonable computation time
{{{#!python
# Generate time domain FIR filter of LTI system                                 
def gentd(sys, td=0, fcut=4, Tl=2, fNyq=8192):
...
    # Apply Hann window to roll off low frequencies                             
    h[0:2*hp] = h[0:2*hp]*cond[0:2*hp]

    h[0]    = 0   # Zero out the DC component                                   
    h[Nf-1] = 0   # Zero out the Nyquist component                              

    # Apply delay                                                               
    delay = np.exp(-1j*w*(np.floor(Tl*fNyq)*dt+td))
    h     = h*delay

    # Fill out negative frequencies                                             
    for k in np.arange(0, Nf-2):
        h = np.append(h, np.conj(h[Nf-k-2]))

    # Take inverse Fourier transform                                            
    td = np.fft.ifft(h)
    td = td.real-td[0].real
}}}
 * The python script is available at git repository
{{{
> kagra-cal/Observation/Phase1/tutorials/genfl
}}}
 * The bode plot shows the the frequency response of the actuation filter
[[attachment:fresp.png|{{attachment:fresp.png|Bode plot to show the frequency response of the actuation filter|width="600"}}]]

 * Check the contents of the output file (filt.npz)
{{{
> python
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> fig=plt.figure()
>>> fl=np.load('filt.npz')
>>> plt.plot(fl['a'])
>>> fig.show()
>>> plt.savefig('afilt.png')
}}}
 * The plot shows the time domain FIR coefficients (= impulse response) for the actuation filter
[[attachment:afilt.png|{{attachment:afilt.png|Time domain FIR coefficients for the actuation filter|width="600"}}]]
 * [[KAGRA/Subgroups/CAL/GstLAL/tutorials/step7|(Next) Step 7: Apply FIR filters and make a strain frame file]]
 * [[KAGRA/Subgroups/CAL/GstLAL/tutorials/step5|(Prev) Step 5: Demodulate a calibration line]]
 * [[KAGRA/Subgroups/CAL/GstLAL/tutorials|Tutorials]]