Frequency Response Functions (FRFs) plots

This notebook is about frequency response functions (FRFs) and the various ways they can be plotted.

Preamble

Dynamic system setup

Frequency response function

Preamble

We will start by setting up the computational environment for this notebook. Since it was created with Python 2.7, we will import a few things from the "future". Furthermore, we will need numpy and scipy for the numerical simulations and matplotlib for the plots:



In [1]:

    
from __future__ import division, print_function

import sys
import numpy as np
import scipy as sp
import matplotlib as mpl

print('System: {}'.format(sys.version))
print('numpy version: {}'.format(np.__version__))
print('scipy version: {}'.format(sp.__version__))
print('matplotlib version: {}'.format(mpl.__version__))









    



System: 3.5.2 |Anaconda custom (64-bit)| (default, Jul  5 2016, 11:41:13) [MSC v.1900 64 bit (AMD64)]
numpy version: 1.11.2
scipy version: 0.18.1
matplotlib version: 1.5.3

We will also need some specific modules and a litle "IPython magic" to show the plots:



In [2]:

    
from numpy import linalg as LA
from scipy import signal
import matplotlib.pyplot as plt

%matplotlib inline

Dynamic system setup

In this example we will simulate a two degree of freedom system (2DOF) as a LTI system. For that purpose, we will define a mass and a stiffness matrix and use proportional damping:



In [3]:

    
MM = np.asmatrix(np.diag([1., 2.]))
print(MM)









    



[[ 1.  0.]
 [ 0.  2.]]



In [4]:

    
KK = np.asmatrix([[20., -10.],[-10., 10.]])
print(KK)









    



[[ 20. -10.]
 [-10.  10.]]



In [5]:

    
C1 = 0.1*MM+0.02*KK
print(C1)









    



[[ 0.5 -0.2]
 [-0.2  0.4]]

For the LTI system we will use a state space formulation. For that we will need the four matrices describing the system (A), the input (B), the output (C) and the feedthrough (D):



In [6]:

    
A = np.bmat([[np.zeros_like(MM), np.identity(MM.shape[0])], [LA.solve(-MM,KK), LA.solve(-MM,C1)]])
print(A)









    



[[  0.    0.    1.    0. ]
 [  0.    0.    0.    1. ]
 [-20.   10.   -0.5   0.2]
 [  5.   -5.    0.1  -0.2]]



In [7]:

    
Bf = KK*np.asmatrix(np.ones((2, 1)))
B = np.bmat([[np.zeros_like(Bf)],[LA.solve(MM,Bf)]])
print(B)









    



[[  0.]
 [  0.]
 [ 10.]
 [  0.]]



In [8]:

    
Cd = np.matrix((1,0))
Cv = np.asmatrix(np.zeros((1,MM.shape[1])))
Ca = np.asmatrix(np.zeros((1,MM.shape[1])))
C = np.bmat([Cd-Ca*LA.solve(MM,KK),Cv-Ca*LA.solve(MM,C1)])
print(C)









    



[[ 1.  0.  0.  0.]]



In [9]:

    
D = Ca*LA.solve(MM,Bf)
print(D)









    



[[ 0.]]

The LTI system is simply defined as:



In [10]:

    
system = signal.lti(A, B, C, D)

To check the results presented ahead we will need the angular frequencies and damping coefficients of this system. The eigenanalysis of the system matrix yields them after some computations:



In [11]:

    
w1, v1 = LA.eig(A)
ix = np.argsort(np.absolute(w1)) # order of ascending eigenvalues
w1 = w1[ix] # sorted eigenvalues
v1 = v1[:,ix] # sorted eigenvectors
zw = -w1.real # damping coefficient time angular frequency
wD = w1.imag # damped angular frequency
zn = 1./np.sqrt(1.+(wD/-zw)**2) # the minus sign is formally correct!
wn = zw/zn # undamped angular frequency
print('Angular frequency: {}'.format(wn[[0,2]]))
print('Damping coefficient: {}'.format(zn[[0,2]]))









    



Angular frequency: [ 1.48062012  4.77574749]
Damping coefficient: [ 0.04857584  0.05822704]