<jupyter_start><jupyter_text>LTIs, FFTs, FRFs, PSDs, CSDs, TFs, ...Table of contents:[Preamble](Preamble)[LTI system setup](LTI-system-setup)[SS system setup](SS-system-setup)[Check system](Check-system)[Noise generation functions](Noise-generation-functions)[Generate noise](Generate-noise)[Histograms](Histograms)[FFT, PSD](FFT,-PSD)[Coloured noise](Coloured-noise)References:[LabVIEW Frequency analysis](http://zone.ni.com/reference/en-XX/help/371361L-01/TOC102.htm)[LabVIEW FRF](http://zone.ni.com/reference/en-XX/help/371361L-01/lvanlsconcepts/frequency_response_function/)[LabVIEW Coherence function](http://zone.ni.com/reference/en-XX/help/371361L-01/lvanlsconcepts/lvac_coherence_function/)[LabVIEW FRF Mag-Phase](http://zone.ni.com/reference/en-XX/help/371361L-01/lvwave/freq_resp_func_mag_phase/)[LabVIEW FRF Real-Im](http://zone.ni.com/reference/en-XX/help/371361L-01/lvwave/freq_resp_func_real_im/) PreambleComputational environment setup<jupyter_code>from __future__ import division, print_function

import numpy as np
from numpy import linalg as LA
from numpy import random

from scipy import signal

from matplotlib import mlab
import matplotlib.pyplot as plt

%matplotlib inline<jupyter_output><empty_output><jupyter_text>LTI system setup[LTI systems](http://en.wikipedia.org/wiki/LTI_system_theory)<jupyter_code>MM = np.matrix(np.diag([1.,2.])) # mass matrix
KK = np.matrix([[20.,-10.],[-10.,10.]]) # stiffness matrix
C1 = 0.1*MM + 0.02*KK # proportional damping matrix<jupyter_output><empty_output><jupyter_text>SS system setup[SS systems](http://en.wikipedia.org/wiki/State-space_representation)<jupyter_code>A = np.bmat([[np.zeros_like(MM), np.identity(MM.shape[0])], [LA.solve(-MM,KK), LA.solve(-MM,C1)]]) # system matrix

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

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

D = Ca*LA.solve(MM,Bf) # feedthrough matrix

system = signal.lti(A, B, C, D)<jupyter_output><empty_output><jupyter_text>Check system<jupyter_code>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]]))<jupyter_output>Angular frequency: [ 1.48062012  4.77574749]
Damping coefficient: [ 0.04857584  0.05822704]<jupyter_text>Bode plot<jupyter_code>w, mag, phase = signal.bode(system)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8,6), sharex=True)
fig.suptitle('Bode plot')
ax1.semilogx(w, mag)
ax1.grid(True)
ax2.semilogx(w, phase)
ax2.grid(True)
plt.show()<jupyter_output>C:\Users\Paulo\Anaconda3\lib\site-packages\scipy\signal\filter_design.py:1092: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  "results may be meaningless", BadCoefficients)<jupyter_text>Real and imaginary plots<jupyter_code>w, H = signal.freqresp(system)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8,6), sharex=True)
fig.suptitle('Real and imaginary plots')
ax1.semilogx(w, H.real)
ax1.grid(True)
ax2.semilogx(w, H.imag)
ax2.grid(True)
plt.show()<jupyter_output>C:\Users\Paulo\Anaconda3\lib\site-packages\scipy\signal\filter_design.py:1092: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  "results may be meaningless", BadCoefficients)<jupyter_text>Noise generation functions[Uniform white noise](http://zone.ni.com/reference/en-XX/help/371361L-01/lvanls/uniform_white_noise/)[Periodic random noise](http://zone.ni.com/reference/en-XX/help/371361L-01/lvanls/periodic_random_noise/)[NumPy irfft function](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.irfft.htmlnumpy.fft.irfft)<jupyter_code>def UWN(n):
    """Uniform white noise."""
    return random.random_sample(n)

def PRN(n):
    """Periodic random noise."""
    values = np.fft.irfft(np.exp(2j*np.pi*random.random_sample(n//2+1)), n=n)
    maxval = np.max(values)
    minval = np.min(values)
    return (values-minval)/(maxval-minval)<jupyter_output><empty_output><jupyter_text>Generate noise<jupyter_code>N = 8192

uwn = UWN(N)
prn = PRN(N)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,6), sharey=True)
ax1.plot(uwn)
ax2.plot(prn)
plt.show()<jupyter_output><empty_output><jupyter_text>Histograms<jupyter_code>fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16,6), sharex=True)
ax1.hist(uwn, bins=100, normed=True, orientation=u'horizontal')
ax2.hist(prn, bins=100, normed=True, orientation=u'horizontal')
plt.show()<jupyter_output><empty_output><jupyter_text>FFT, PSD<jupyter_code>Fs = 1.
freq = np.fft.fftfreq(N, d=1./Fs)

UU = np.fft.fft(uwn)
PP = np.fft.fft(prn)

(Puu, freqs) = mlab.psd(uwn, NFFT=N, Fs=Fs, window=mlab.window_none, sides='twosided', scale_by_freq=True)
(Ppp, freqs) = mlab.psd(prn, NFFT=N, Fs=Fs, window=mlab.window_none, sides='twosided', scale_by_freq=True)

fig, ax = plt.subplots(2, 2, figsize=(16,12))
ax[0, 0].plot(freq, UU.real)
ax[0, 0].grid(True)
ax[0, 1].plot(freq, PP.real)
ax[0, 1].grid(True)
ax[1, 0].semilogy(freqs, Puu)
ax[1, 0].grid(True)
ax[1, 1].semilogy(freqs, Ppp)
ax[1, 1].grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Coloured noise<jupyter_code>w, H = signal.freqresp(system, w=2.*np.pi*freq)

AAA = UU*H
aaa = np.fft.ifft(AAA)

plt.figure()
plt.plot(w, AAA.real)
plt.show()

plt.figure()
plt.plot(aaa)
plt.show()<jupyter_output>C:\Users\Paulo\Anaconda3\lib\site-packages\scipy\signal\filter_design.py:1092: BadCoefficients: Badly conditioned filter coefficients (numerator): the results may be meaningless
  "results may be meaningless", BadCoefficients)
{'fraction_non_alphanumeric': 0.1684470008216927, 'fraction_numerical': 0.037633525061626955, 'mean_word_length': 6.51358024691358, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '294420', 'license': 'MIT', 'n_tokens_mistral': 2635, 'n_tokens_neox': 2443, 'n_words': 391, 'path': 'Acronyms.ipynb', 'repo_name': 'pxcandeias/py-notebooks', 'size': '6085'}