<jupyter_start><jupyter_text>Universidade Federal do Rio Grande do Sul (UFRGS)   Programa de Pós-Graduação em Engenharia Civil (PPGEC)    PEC00144: Experimental Methods in Civil Engineering Part II: Instrumentation[8.   Analog signals processing](section_8)  &nbsp;&nbsp;&nbsp;[8.1. Autocovariance function and stationarity](section_81)     &nbsp;&nbsp;&nbsp;[8.2. Fourier series and Fourier transform](section_82)   &nbsp;&nbsp;&nbsp;[8.3. Power spectral density and periodograms](section_83)   &nbsp;&nbsp;&nbsp;[8.4. White noise and pink noise](section_84)   &nbsp;&nbsp;&nbsp;[8.5. Signal derivation and integration](section_85)   &nbsp;&nbsp;&nbsp;[8.6. Basic stationarity test](section_86)   &nbsp;&nbsp;&nbsp;[8.7. Level upcrossing rate and peak factor](section_87)   &nbsp;&nbsp;&nbsp;[8.8. Frequency domain signal filtering](section_88)   &nbsp;&nbsp;&nbsp;[8.9. Cross correlation and coherence functions](section_89)   ---_Prof. <NAME>, Dr.techn._ [(ORCID)](https://orcid.org/0000-0001-5640-1020)  _Porto Alegre, RS, Brazil_<jupyter_code># Importing Python modules required for this notebook
# (this cell must be executed with "shift+enter" before any other Python cell)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from   MRPy import MRPy<jupyter_output><empty_output><jupyter_text>8. Analog signals processing    8.1. Autocovariance function and stationarity   Analog signals are the outcome of some physical measurement, usually some electrical quantity like voltage, current, or electric charge. The magnitude of such quantities are expected to be _analogous_ to some observed physical quantity, like stress, force, displacement, acceleration, etc. An analog signal is generally considered as a random process, $X(t)$. A random processis a time dependent random variable, what requires its statistical properties to be alsoregarded as time dependent.The autocovariance function, $C_X(\tau)$, of a random process, $X(t)$, is defined as thefirst cross moment between the process amplitude at two time instants:$$ C_X(t_1, t_2) = {\rm E}\left\{ X(t_1) X(t_2) \right\} - \mu_1\mu_2$$where $\mu_1$ and $\mu_2$ are the mean value of $X(t)$ at instants $t_1$ and $t_2$,respectively.If a random process is _stationary_, its statistical properties are assumed to be_independent of time_, and the expression above depends only on the time gap, $\tau = t_2 - t_1$:$$ C_X(\tau) = {\rm E}\left\{ X(t) X(t + \tau) \right\} - \mu_X^2 $$The definitions above exclude the mean value of $X(t)$, keeping only the timedependent part of its amplitude. Furthermore, it is evident that the autocovariancefunction is a _pair function_ (symmetric over the axis $\tau = 0$) and that at origin:$$ C_X(0) = {\rm E}\left\{ X^2(t) \right\} - \mu_X^2 = \sigma_X^2$$which is the process _variance_.The autocovariance function can be normalized by $\sigma_X^2$ resulting in theprocess _autocorrelation function_:$$ R^2_X(\tau) = \frac{C_X(\tau)}{\sigma_X^2} $$which is also a pair function such that $-1 \leq R_X(\tau) \leq +1$.  8.2. Fourier series and Fourier transform   An introduction to Fourier analysis can be found in [Class 7](https://nbviewer.jupyter.org/github/mmaiarocha/PEC00025/blob/master/resources/Class_07_FourierTransform.ipynb?flushcache=true) of the course [Introduction to Vibration Theory](https://github.com/mmaiarocha/PEC00025).  8.3. Power spectral density and periodograms   The power spectral density, $S_X(\omega)$, of a stationary random process is defined as the Fourier transform of the autocovariance function, $C_X(\tau)$:$$ S_X(\omega) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} {C_X(\tau) e^{-i\omega\tau} \, d\tau}$$This means that $S_X(\omega)$ and $C_X(\tau)$ are a Fourier transform pair, and consequently:$$ C_X(\tau) = \int_{-\infty}^{+\infty} {S_X(\omega) e^{i\omega\tau} \, d\omega} $$The definitions above implies that:$$ C_X(0) = \int_{-\infty}^{+\infty} {S_X(\omega) \, d\omega} = \sigma_X^2$$what means that the total integral of the spectral density is, by definition,the process variance.The ``MRPy`` module provides a straightforward method for visualizing the spectral density estimator, called _periodogram_, and the corresponding autocorrelationfunction. Below is a short script where a cosine wave is analysed:<jupyter_code>Td   =  2.           # total signal duration (s)
fs   =  512          # sampling rate (Hz)
N    =  int(Td*fs)   # total number of samples

f0   =  8.           # cosine wave frequency (Hz)

X0   =  MRPy.harmonic(NX=1, N=N, fs=fs, X0=1.0, f0=f0, phi=0.0)
 
f00  =  X0.plot_time(fig=0, figsize=(8,3), axis_t=[0, Td, -2.0, 2.0])
f01  =  X0.plot_freq(fig=1, figsize=(8,3), axis_f=[0, 16,  0.0, 2.0])    # periodogram
f02  =  X0.plot_corr(fig=2, figsize=(8,3), axis_T=[0,  1, -2.0, 2.0])    # autocovariance<jupyter_output><empty_output><jupyter_text>It can be seen that the periodogram of a sine signal is a impulse function at the sine frequency. Below it is show how to simulate a process with given autocorrelation function.<jupyter_code>Sx, fs =  X0.periodogram()
X1     =  MRPy.from_periodogram(Sx, fs)

f00    =  X1.plot_time(fig=0, figsize=(8,3), axis_t=[0, Td, -2.0, 2.0])
f01    =  X1.plot_freq(fig=1, figsize=(8,3), axis_f=[0, 16,  0.0, 2.0])    # periodogram
f02    =  X1.plot_corr(fig=2, figsize=(8,3), axis_T=[0,  1, -2.0, 2.0])    # autocovariance

tau    =  X1.T_axis()
kt     = (tau <= 0.5) 
Cx     =  np.zeros(tau.shape)
Cx[kt] =  1 - 2*tau[kt]

X2     =  MRPy.from_autocov(Cx, tau[-1])

f00    =  X2.plot_time(fig=0, figsize=(8,3), axis_t=[0, 8, -4.0,  4.0])
f01    =  X2.plot_freq(fig=1, figsize=(8,3), axis_f=[0, 8,  0.0,  2.0])
f02    =  X2.plot_corr(fig=2, figsize=(8,3), axis_T=[0, 4, -2.0,  2.0])<jupyter_output><empty_output><jupyter_text>8.4. White noise and pink noise  A _Gaussian white noise_ is a signal with random amplitude with normal distribution and a constant power density all over the frequency domain:$$ S_X(\omega) = S_0 $$The associated autocorrelation function is a Dirac's Delta at the origin, with an_infinite_ pulse integral.This signal in practice must have a limited band, otherwise the corresponding variance would be infinite. For practical purposes a signal is considered to be a white noiseif the power is constant (within some statistical error) over some relevantfrequency band $\Delta\omega = \omega_2 - \omega_1$:$$ S_X(\omega) = S_0, \hspace{5mm} {\rm for}                       \hspace{5mm} \omega_1 \leq \omega \leq \omega_2 $$ and zero otherwise. The corresponding autocorrelation function is:$$C_X(\tau) = \frac{4S_0}{\tau} \sin\left( \frac{ \Delta\omega}{2} \tau \right)                               \cos\left( \omega_0 \tau \right) $$where $\omega_0 = (\omega_1 + \omega_2)/2$ is the band center.The corresponding variance is:$$ \sigma_X^2 = C_X(0) = 2\Delta\omega S_0$$As the band width $\Delta\omega$ decreases, the signal above approaches a cosine wave with frequency $\omega_0$, as described in the previous section.Let us take a look on a band-limited Gaussian white noise simulation with ``MRPy``.The simulation uses $S_0 = 1$, hence the standard deviation will be $\sigma_X = \sqrt{2\Delta\omega}$.<jupyter_code>X = MRPy.white_noise(NX=1, N=2**15, fs=512, band=[8, 10])

print('Mean value:         {0:7.4f}'.format(X[0].mean()))
print('Standard deviation: {0:7.4f}'.format(X[0].std()),'\n')

f03 = X.plot_time(fig=3, figsize=(8,3), axis_t=[0, Td, -10.00, 10.00])
f04 = X.plot_freq(fig=4, figsize=(8,3), axis_f=[0, 32,   0.00,  1.00])
f05 = X.plot_corr(fig=5, figsize=(8,3), axis_T=[0,  1,  -1.00,  1.00])<jupyter_output>Mean value:          0.0000
Standard deviation:  1.0000<jupyter_text>The ``MRPy`` module uses a simulation technique that gives an almost perfectlyconstant periodogram, as specified.  8.5. Signal derivation and integration  As a starting point, let us calculate the following derivative:$$ \frac{d}{d\tau}\left[ X(t) X(t+\tau) \right] =                          X(t) \cdot \frac{d X(t + \tau)}{d(t+\tau)} \cdot \frac{d(t+\tau)}{dt} =                         X(t) \dot{X}(t+\tau) $$Now, using the expected value operator and considering autocovariance symmetry:$$ \frac{d C_X(\tau)}{d\tau} = {\rm E}\left\{ X(t) \dot{X}(t+\tau) \right\}                              = {\rm E}\left\{ X(t-\tau) \dot{X}(t) \right\} $$and following the logic to find the second derivative:$$ \frac{d^2 C_X(\tau)}{d\tau^2} = -{\rm E}\left\{\dot{X}(t) \dot{X}(t+\tau) \right\}                                  = - C_{\dot{X}} (\tau)$$With this result at hand, we go back to the relation between power density and autocovariance:$$ C_X(\tau) = \int_{-\infty}^{+\infty} {S_X(\omega) e^{i\omega\tau} \, d\omega} $$and apply double derivative:\begin{align*}\frac{d C_X(\tau)}{d\tau}     &=    \int_{-\infty}^{+\infty} {i\omega  S_X(\omega) e^{i\omega\tau} \, d\omega}   \\\frac{d^2 C_X(\tau)}{d\tau^2} &= - \int_{-\infty}^{+\infty} {\omega^2 S_X(\omega) e^{i\omega\tau} \, d\omega}\end{align*}Now, considering that the following relation is also valid:$$ C_\dot{X}(\tau) = \int_{-\infty}^{+\infty} {S_\dot{X}(\omega) e^{i\omega\tau} \, d\omega} $$It finally results that:\begin{align*}S_\dot{X} (\omega) &= \omega^2 S_{X} (\omega) \\S_\ddot{X}(\omega) &= \omega^4 S_{X}(\omega)\end{align*}These relations allow us to calculate the spectral density of velocity andacceleration processes from the spectral density of displacement process, or vice-versa.They are quite useful for converting signal amplitudes obtained with one type of transducer (for instance, an accelerometer) to amplitudes as they would have beenobtained with other type of transducer (for instance, displacement).The ``MRPy`` class provides the calculation of derivatives and integrals infrequency domain, as demonstrated below. The methods allow the definitionof a passing frequency band, for simultaneously eliminating noise errors.<jupyter_code>acc =  MRPy.white_noise(N=2**15, fs=512) + 0.001*MRPy.harmonic(N=2**15, fs=512, f0=0.5)
Td  =  acc.Td

f03 =  acc.plot_time(fig=3, figsize=(8,3), axis_t=[0,  Td, -10.00, 10.00])
f04 =  acc.plot_freq(fig=4, figsize=(8,3), axis_f=[0,   1,   0.00,  0.01])
f05 =  acc.plot_corr(fig=5, figsize=(8,3), axis_T=[0, 0.1,  -1.00,  1.00])
V1  =  acc.integrate()               # integration without filtering
#V1  =  acc.differentiate()           # integration without filtering

f03 =  V1.plot_time(fig=3, figsize=(8,3))#, axis_t=[0,  Td, -0.50, 0.50])
f04 =  V1.plot_freq(fig=4, figsize=(8,3))#, axis_f=[0,  2,   0.00, 0.10])
f05 =  V1.plot_corr(fig=5, figsize=(8,3))#, axis_T=[0, 32,  -1.00, 1.00])
V2   =  acc.integrate(band=[0.3, 0.7])         # integration with filtering

f06  =  plt.figure(6, figsize=(8,3))
f06a =  plt.plot(V1.t_axis(), V1[0], 'r')
f06b =  plt.plot(V2.t_axis(), V2[0], 'b')

#plt.axis([0, Td, -1, 1])
plt.grid(True)<jupyter_output><empty_output><jupyter_text>Whenever a signal is integrated without a low frequency cut, a _zero drift_ is expected to happen! This can be avoided by setting a lower frequencybound as high as possible without attenuating the useful part of the signal.  8.6. Basic stationarity test   Previously we have demonstrated that:$$ \frac{d C_X(\tau)}{d\tau} = {\rm E}\left\{ X(t) \dot{X}(t+\tau) \right\} =     \int_{-\infty}^{+\infty} {i\omega  S_X(\omega) e^{i\omega\tau} \, d\omega} $$By making $\tau = 0$ in the relations above gives:$$ {\rm E}\left\{ X(t) \dot{X}(t) \right\} =     \int_{-\infty}^{+\infty} {i\omega  S_X(\omega) \, d\omega} $$The power spectral density is a _pair function_. Multiplying it by $i\omega$necessarily yields an _odd function_. Consequently the integral vanishes and:$$ {\rm E}\left\{ X(t) \dot{X}(t) \right\} = 0 $$This relation can be a shortcut for ascertaining the signal stationarity and hencevalidating the constancy of its spectral density.This result can be demonstrated through simulation with the ``MRPy`` module.A Gaussian white noise is simulated and its derivative is calculated. Bothprocedures restrain the signals to the same limited band:<jupyter_code>X    =  MRPy.white_noise(N=2**15, fs=512, band=[8, 16])
Xdot =  X.differentiate(band=[8, 16])

f07  =  plt.figure(7, figsize=(6,6))
f07a =  plt.plot(X[0], Xdot[0], 'b.')

plt.axis([-15, 15, -500, 500])
plt.grid(True)<jupyter_output><empty_output><jupyter_text>The lack of correlation between the process and its derivative observed in the plot above is an evidence of process stationarity.  8.7. Level upcrossing and peak factors   For a stationary process, the expected number of upcrossings, $N^{+}_a(T)$, of a givenamplitude level, $a$, within a given observation time, $T$, is given by:$$ N^{+}_a(T) = \nu^{+}_a T$$where $\nu^{+}_a$ is the _upcrossing rate_ of level $a$. In other words, the number ofupcrossings is proportional to the observation time.The upcrossing rate is calculated by integrating the joint probability distribution of the process amplitude and its derivative, such that the amplitude is fixed atlevel $a$ and only positive values are regarded for its derivative:$$ \nu^{+}_a = \int_0^{\infty} {p_{X\dot{X}}(x=a, \dot{x}) \, x \, d{\dot{x}}} $$This result can be particularized for a Gaussian process, also considering thatstationarity implies that the process and its derivative are uncorrelated:$$ \nu^{+}_a = \frac{1}{2\pi} \, \frac{\sigma_\dot{X}}{\sigma_X} \,               \exp \left( -\frac{a^2}{2\sigma_X^2} \right) $$Recalling from previous results that:$$ \sigma_X^2 = \int_{-\infty}^{+\infty} {S_X(\omega) \, d\omega} $$and that:$$ \sigma_\dot{X}^2 = \int_{-\infty}^{+\infty} {\omega^2 S_X(\omega) \, d\omega} $$allow us to estimate the upcrossing rate for any stationary signal byintegrating the periodogram as indicated above.Para a parcela flutuante da resposta em deslocamento, como a análise foi conduzida no domínio da frequência, resultados no domínio do tempo só pode ser obtidos em termos estatísticos. A abordagem da NBR-6123 consiste em se adotar um fator de pico da resposta modal $g_k = 4$, aplicado sobre o desvio padrão da resposta flutuante.No entanto, já que o espectro da resposta está disponível, é possível obter-se uma estimativa mais precisa deste fator de pico utilizando-se a fórmula de Davenport:$$g_k = \sqrt{2 \ln (\nu^{+}_0 T)} + \frac{0.5772}{\sqrt{2 \ln (\nu^{+}_0 T)}}$$onde $T$ é o tempo de observação, adotado como 600s na NBR-6123, e $\nu_k$ é a taxa de cruzamento do nível zero para o positivo (_zero upcrossing rate_), calculada a partir do espectro como:$$ \nu^{+}_0 = \sqrt{\frac{\int_0^\infty{f^2 S_X(f) \; df}}                    {\int_0^\infty{    S_X(f) \; df}}} $$Observe-se que o denominador dentro da raiz é a variância da resposta modal. Aplicando-se as expressões acima ao exemplo de cálculo tem-se:<jupyter_code>X   =  MRPy.white_noise(N=2**16, fs=256)   # full band white noise
Y   =  X.sdof_Duhamel(4, 0.1)              # mass-spring system response
Y   =  Y/Y[0].std()                        # forces deviation equal to 1

Tm  =  X.Td
gD  =  Y.Davenport(T=600)                 # Gaussian narrow-band process
gS  =  Y.splitmax (T=600)                  # Prof. Rocha's method

Ypk = Y[0].mean() + gD[0]*Y[0].std()      # mean + g.deviation

print('Peak factor from Davenport:  {0:6.3f}'.format(gD[0]))
print('Peak factor from splitmax:   {0:6.3f}'.format(gS[0]))
print('Peak value for displacement: {0:6.3f}'.format(Ypk))

f09  =  plt.figure(9, figsize=(12,4))
f09a =  plt.plot(Y.t_axis(), Y[0], 'r', lw=1)

plt.axis([0, X.Td, -6, 6])
plt.grid(True)<jupyter_output>Peak factor from Davenport:   4.092
Peak factor from splitmax:    3.899
Peak value for displacement:  4.101<jupyter_text>8.8. Frequency domain signal filtering<jupyter_code>Td   =  4.           # total signal duration
fs   =  2048         # sampling rate
N    =  int(Td*fs)   # total number of samples
f60  =  60.          # sinus wave frequency

t    =  np.linspace(0, Td, N) 
xi   =  2*np.random.randn(N) + 1*np.cos(2*np.pi*f60*t)
X    =  MRPy(xi, fs=fs)

Y    =  X #X.filtered(band=[59.8,60.2], mode='pass') # try out also 'stop'

f08  =  plt.figure(8, figsize=(12,4))
f08a =  plt.plot(X.t_axis(), X[0], 'r', lw=0.5)
f08b =  plt.plot(Y.t_axis(), Y[0], 'b', lw=0.8)

plt.axis([0, 1, -10, 10])
plt.grid(True)

Y.plot_freq(fig=9)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1510272071071627, 'fraction_numerical': 0.03960762539330002, 'mean_word_length': 4.003086419753086, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 32, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '208503', 'license': 'MIT', 'n_tokens_mistral': 5989, 'n_tokens_neox': 5605, 'n_words': 1774, 'path': '08_AnalogSignals.ipynb', 'repo_name': 'mmaiarocha/PEC00144', 'size': '16219'}