1. Introduction to signal processing

The introductory chapter presents the motivation for the book, defines some terms and outlines the structure of the book.

The chapters of the book are:

• 2. Fourier series,

• 3. The Fourier transform,

• 4. Linear time-invariant systems with one degree of freedom,

• 5. Data quantization,

• 6. The discrete Fourier transform,

• 7. Fundamentals of random processes,

• 8. Response of linear systems to random excitation.

What is signal processing?

We use the term signal processing to denote various methods of processing measured data with the aim of revealing information about the system under consideration. We could potentially also speak of data processing, but such a term is on the one hand too broad and on the other too narrow. Too broad because with signal processing we wish to restrict ourselves to the identification/characterization of various engineering systems and not to address the general processing of arbitrary data. Too narrow because data are always discrete, whereas processing can already begin with signals, which may be either continuous or discrete.

As an example of signal processing we can highlight the identification of bearing faults or the identification of dynamic properties of structures (e.g. damping, stiffness). In both cases we must process the measured signals correctly, even before we sample them in time and then quantize them (see 5. Data quantization). The figure below shows a measurement of the force on an undamaged bearing and on two damaged bearings (Slavič et al. [2011]). In the measurements we see certain differences, but can we identify the type of damage from the measurements? Furthermore, can we identify the size of the damage from the measurements? We will try to answer this and similar questions in this book.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

new = np.load('./data/bearing/new 1A.npy')
radial = np.load('./data/bearing/radial Low 1A.npy')
contamination = np.load('./data/bearing/contam low 2A.npy')
time = np.arange(len(new)) * 1./5000

sel = time<0.5
plt.plot(time[sel], new[sel, 2], alpha=0.5, label='Undamaged bearing')
plt.plot(time[sel], radial[sel, 2], alpha=0.5, label='Damaged bearing 1')
plt.plot(time[sel], contamination[sel, 2], alpha=0.3, label='Damaged bearing 2')
plt.xlabel('Time [s]')
plt.ylabel('Force [N]')
plt.ylim(-0.1, 0.1)
plt.grid()
plt.legend()
plt.show()

As a next example, the figure below shows the frequency response function (FRF), which in the frequency domain relates the excitation of a system to its response (the example is taken from the pyUFF package). Correct identification of the frequency response function is essential for the correct characterization of dynamic systems and is therefore also one of the fundamental learning objectives of this book.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import pyuff

uff = pyuff.UFF('./data/uff/beam.uff').read_sets(setn=4)
fr = uff['x']
H1 = uff['data']

sel = fr<500
plt.semilogy(fr[sel], np.abs(H1[sel]), label='H1')
plt.xlabel('Frequency [Hz]')
plt.ylabel('FRF [m$\\,$s$^{-2}$/N]')
plt.grid()
plt.legend()
plt.show()

What is a system?

We usually try to describe dynamic and other systems on the basis of models, and to identify the model parameters with the help of experimental measurements that we process appropriately. The figure below shows an arbitrary system whose excitation is described in time by \(x(t)\) and whose response by \(y(t)\). As an example of the system excitation \(x(t)\) we can imagine a force acting on some mass, and as the response \(y(t)\) the voltage on an acceleration sensor or a sound-pressure sensor. The response is often measured via an electrical voltage that is usually proportional to the measured quantity (e.g. acceleration, sound pressure); of course, both the excitation and the response may also be related to other quantities that we can measure, e.g. current or the intensity of a picture element in an image (see Gorjup et al. [2021]).

The system, more precisely the linear time-invariant system, will be defined in detail in the chapter 4. Linear time-invariant systems with one degree of freedom.

As an example, let us look at the response of a single-degree-of-freedom dynamic system to a half-sine disturbance with noise. The impulse response function of a single-degree-of-freedom system is defined as:

\[ h(t)=\frac{1}{\omega_{0,d}}\,\mathrm{e}^{-\delta\,\omega_0\,t}\,\sin(\omega_{0,d}\,t). \]

The concrete implementation and the result are shown in the code below.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

T = 1.5
N = 1000
w = 0.1
A = 0.1
t = np.linspace(-T/2, T/2, N, endpoint=False)
dt = t[1] - t[0]
x = A*np.cos(2*np.pi*t/(2*w*T))
x[np.logical_or(-w*T*0.5>t, t>w*T*0.5)] = 0.
x = x + 0.2*A*(np.random.rand(N)-0.5)

delta = 0.1 # damping ratio
omega0 = 2*np.pi*10
omega0d = omega0*np.sqrt(1-delta*delta)
th = dt*np.arange(len(x))
h = np.exp(-delta*omega0*th)*np.sin(omega0d*th)/omega0d

plt.plot(t,x, label='Half-sine disturbance with noise')
plt.plot(th,h, label='Impulse response')
plt.legend()
plt.show()

hx = np.convolve(h, x) # mode='full' -> sum of both arrays
plt.title('Response to a half-sine disturbance with noise')
plt.plot(hx, 'C2')
plt.show()

So that we can understand systems in the time and frequency domains, we will first have to become acquainted with the Fourier series presented later (2. Fourier series) and the Fourier integral transform (3. The Fourier transform).

A system can be uniquely identified if we identify the frequency response function \(H(f)\) or its time-domain counterpart, the impulse response function \(h(t)\). If the excitation \(X(f)\) and the response \(Y(f)\) are known in the frequency domain, the identification is theoretically extremely simple: \(H(f)=Y(f)/X(f)\). In practice, the problem is that the measurement of the excitation and the response is burdened with measurement and other uncertainties. An example of a real system is shown in the figure below; \(n_x(t)\) and \(n_y(t)\) represent the noise on the excitation and response sides, and \(x_m(t)\) and \(y_m(t)\) the measured excitation and the measured response (with noise).

The measured excitation \(x_m(t)\) is obviously different from the one exciting the system, \(x(t)\), and the measured response \(y_m(t)\) is obviously different from the actual response \(y(t)\). Given the measurement and other uncertainties, we do not have access to the actual excitation \(x(t)\) and the actual response \(y(t)\), and therefore we cannot determine the frequency response function \(H(f)\). In the chapter Estimators of the frequency response function we will look at how to estimate \(H(f)\) on the basis of quantities that we can measure (\(x_m(t)\) and \(y_m(t)\)); we will then speak of the estimator of the frequency response function \(\tilde{H}(f)\).

Continuous/discrete signals

Engineering processes are usually continuous and we sense them with the help of various sensors that generate a physically measurable quantity, the so-called signal. A signal is a continuous quantity (sometimes we will also hear the name analog signal). Because we process the data with a computer, we discretize these continuous signals (see the chapter Uniform time sampling). Such discretization is usually done with a constant time step. The figure below shows a continuous signal and a discrete time series.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt


# on a computer we cannot actually generate a continuous time series
# the "continuous" one will just be much denser
t, dt = np.linspace(0,4,1000, retstep=True)
t2, dt2 = np.linspace(0,4,9, retstep=True)

plt.plot(t, np.sin(t**2), label='Continuous signal') 
plt.plot(t2, np.sin(t2**2), 'o',label='Discrete time series $\\Delta t=$0.5 s')
plt.xlabel('Time [s]')
plt.ylabel('$\\sin(t^2)$')
plt.grid()
plt.legend()
plt.show()

Classification of signals

For the treatment that follows, the classification of signals will be important (adapted from Bendat and Piersol [2011], Shin and Hammond [2008]), as shown in the figure below.

We divide signals into:

• deterministic signals \(-\) these are signals whose values in time are uniquely determined, and

• random signals \(-\) these are signals whose value we do not know exactly, but for which we can determine a probability distribution.

An example of a deterministic signal is the function:

\[ x(t) = \sin(t), \]

while an example of a random signal is a measurement of the surface height of a rough sea. When classifying signals we must not forget that we can have processes composed of a deterministic and a random part (e.g. a measurement of a product’s profile also includes a random measurement error).

We further divide deterministic signals into periodic ones, which repeat after a certain time (a period), and those that do not repeat, i.e. aperiodic ones. Periodic signals are divided into harmonic ones, which are defined by a sinusoid of arbitrary phase delay, and complex periodic ones, which are defined by a sum of several harmonic components. For complex periodic signals it is important that the ratio of the frequencies of the harmonic components is a rational number, since in that case the signal repeats with the least common period. If the ratio of the harmonic components is not a rational number, we speak of almost periodic signals, which belong to the group of aperiodic signals. This group also includes transient and chaotic signals.

Understanding the periodicity of signals and the influence of aperiodicity on the Fourier series (2. Fourier series) and the Fourier transform (3. The Fourier transform) is of key importance for correct interpretation in the frequency domain.

We further divide random signals according to whether the statistical distribution that defines the random processes changes with time (non-stationary) or does not (stationary); we will deal with this in more detail in the chapter Moments of the probability density function.

Some examples of deterministic signals

The figure below shows periodic signals: first a triangular signal that repeats with a period of 1 s, then a harmonic signal that likewise repeats with a period of 1 s.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt


t, dt = np.linspace(0,4,1000, retstep=True)

plt.title('Periodic signals')
plt.plot(t, t%1, label='Periodic signal', alpha=0.8) 
plt.plot(t, np.sin(2*np.pi*t+1), label='Harmonic signal', alpha=0.8) 
plt.ylim(-2.1,2.1)
plt.xlabel('Time [s]')
plt.grid()
plt.legend()
plt.show()

The last signal represents a complex periodic signal \(-\) these signals are composed as a sum of several harmonic components; the period of the one in the figure is 2 s.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt


t, dt = np.linspace(0,10,1000, retstep=True)

plt.title('Deterministic signals')
plt.plot(t, np.sin(2*np.pi*t+1)+np.sin(3*np.pi*t+1), label='Complex periodic signal', alpha=0.8)
plt.plot(t, np.sin(2*np.pi*t+1)+np.sin(3.075*np.pi*t+1), label='Almost periodic signal', alpha=0.8)
plt.ylim(-2.1,2.1)
plt.xlabel('Time [s]')
plt.grid()
plt.legend(loc=8)
plt.show()

For aperiodic signals, we will show as an example a transient signal that arises when a dynamic system is excited by an impulse disturbance.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt


t, dt = np.linspace(0,4,1000, retstep=True)
omega = 10*np.pi
damping = 0.05

plt.title('Aperiodic signal')

plt.plot(t, np.exp(-omega*damping*t)*np.sin(omega*t+1), label='Transient phenomenon', alpha=0.8) 
plt.ylim(-1.1,1.1)
plt.xlabel('Time [s]')
plt.grid()
plt.legend()
plt.show()

We will not treat chaotic aperiodic signals in more detail here; for them it holds that, because of their high degree of dynamics, they can be reliably predicted only for a relatively short time ahead, after which their behavior seems increasingly random.

When exciting dynamic systems we often use a frequency sweep (also called a chirp); below, a 4 s long linear (the frequency increases monotonically) chirp from 1 to 5 Hz is shown.

Show code cell source

Hide code cell source

import numpy as np
import matplotlib.pyplot as plt
import pyExSi as es

fs = 44000
T = 4
t = np.linspace(0,T,fs*T+1)
x, _, fr = es.sine_sweep(time=t, freq_start=1, freq_stop=5, mode='linear', freq_return=True)

fig, ax1 = plt.subplots()
plt.title('Aperiodic signal: frequency sweep/chirp')

ax1.set_xlabel('Time [s]')
ax1.set_ylabel('$x(t)$ [-]', color='C0')
ax1.plot(t, x, color='C0', alpha=0.8)
ax1.tick_params(axis='y', labelcolor='C0')

ax2 = ax1.twinx() 
ax2.set_ylabel('Frequency [Hz]', color='C1')
ax2.plot(t, fr, color='C1', alpha=0.8) 
ax2.tick_params(axis='y', labelcolor='C1')
plt.show()

When analyzing signals it sometimes helps to hear them; below we can hear a 4 s long logarithmic chirp from 100 to 16000 Hz.

Show code cell source

Hide code cell source

import numpy as np
import pyExSi as es
from IPython.display import Audio, display

fs = 44000
T = 4
t = np.linspace(0,T,fs*T+1)
x, _, fr = es.sine_sweep(time=t, freq_start=100, freq_stop=16000, mode='logarithmic', freq_return=True)
display(Audio(data=x, rate=fs))

plt.plot(t, fr, label='Logarithmic frequency change')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.grid()
plt.legend()
plt.show()

Some examples of random signals

The figure below shows a random signal with a uniform distribution. The left-hand figure shows the first 50 ms of the data in time, while the right-hand figure shows the power spectral density (see Power spectral density (PSD)). We can observe that, in the frequency domain as well, the uniform distribution has a relatively constant power across the entire analyzed frequency range (the data are generated with a sampling frequency of 44 kHz).

Show code cell source

Hide code cell source

import numpy as np
import pyExSi as es
from IPython.display import Audio, display

fs = 44000
T = 4
t = np.linspace(0,T,fs*T+1)
N = len(t)
x = es.uniform_random(N=N)
sel = t<0.05
display(Audio(data=x, rate=fs))


fig, axs = plt.subplots(1,2)
fig.suptitle('Stationary signal: uniform distribution')
plt.tight_layout(pad=1, w_pad=3)
axs[0].plot(t[sel], x[sel])
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Response in time (first 50 ms only)')
axs[1].psd(x, Fs=fs)
axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power spectral density [dB/Hz]')
axs[1].set_ylim(-80, -20)
axs[1].set_yticks(np.arange(-80,-20,5))
plt.show()

Because the energy is distributed across a (too) wide frequency range, exciting dynamic systems with a uniform distribution is rarely appropriate. We usually use a normal (or Gaussian) distribution with a defined frequency spectrum. The figure below shows a normal distribution with uniform power density in the frequency range from 1000 to 4000 Hz.

Show code cell source

Hide code cell source

import numpy as np
import pyExSi as es
from IPython.display import Audio, display

fs = 44000
T = 4
t = np.linspace(0,T,fs*T+1)
N = len(t)

M = N // 2 + 1 # number of data points of frequency vector
freq = np.arange(0, M, 1) * fs / N # frequency vector
freq_lower = 1000 # PSD lower frequency limit  [Hz]
freq_upper = 4000 # PSD upper frequency limit [Hz]

PSD = es.get_psd(freq, freq_lower, freq_upper) # one-sided flat-shaped PSD
x = es.random_gaussian(N, PSD, fs)

sel = t<0.05
display(Audio(data=x, rate=fs))


fig, axs = plt.subplots(1,2)
fig.suptitle(f'Stationary signal: normal distribution in the range from {freq_lower} to {freq_upper} Hz')
plt.tight_layout(pad=1, w_pad=3)
axs[0].plot(t[sel], x[sel])
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Response in time (first 50 ms only)')
axs[1].psd(x, Fs=fs)
axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power spectral density [dB/Hz]')
axs[1].set_xlim(0, 5000)
axs[1].set_ylim(-80, -20)
axs[1].set_yticks(np.arange(-80,-20,5))
plt.show()

For comparison, let us also look at a normal distribution with uniform power density in the frequency range from 100 to 400 Hz.

Show code cell source

Hide code cell source

import numpy as np
import pyExSi as es
from IPython.display import Audio, display

fs = 44000
T = 4
t = np.linspace(0,T,fs*T+1)
N = len(t)

M = N // 2 + 1 # number of data points of frequency vector
freq = np.arange(0, M, 1) * fs / N # frequency vector
freq_lower = 100 # PSD lower frequency limit  [Hz]
freq_upper = 400 # PSD upper frequency limit [Hz]

PSD = es.get_psd(freq, freq_lower, freq_upper) # one-sided flat-shaped PSD
x = es.random_gaussian(N, PSD, fs)

sel = t<0.05
display(Audio(data=x, rate=fs))


fig, axs = plt.subplots(1,2)
fig.suptitle(f'Stationary signal: normal distribution in the range from {freq_lower} to {freq_upper} Hz')
plt.tight_layout(pad=1, w_pad=3)
axs[0].plot(t[sel], x[sel], label='Uniform distribution')
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Response in time (first 50 ms only)')
axs[1].psd(x, Fs=fs)
axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power spectral density [dB/Hz]')
axs[1].set_xlim(0, 5000)
axs[1].set_ylim(-80, -20)
axs[1].set_yticks(np.arange(-80,-20,5))
plt.show()

Let us also look at a non-stationary signal. We start by preparing a stationary signal with uniform power density in the frequency range from 1000 to 4000 Hz and a non-normal distribution (with a kurtosis parameter = 5). We then modulate the stationary signal with a carrier signal in the frequency range from 10 to 20 Hz in order to obtain a non-stationary signal. Even though the power density is distributed similarly to above, the non-normal and non-stationary signal usually sounds considerably more unpleasant; in the case of vibration fatigue, non-stationary loading is usually also considerably more damaging.

Show code cell source

Hide code cell source

import matplotlib.pyplot as plt
import numpy as np
import pyExSi as es
from IPython.display import Audio, display

fs = 44000
T = 4
t = np.linspace(0,T,fs*T+1)
N = len(t)

M = N // 2 + 1 # number of data points of frequency vector
freq = np.arange(0, M, 1) * fs / N # frequency vector
freq_lower = 1000 # PSD lower frequency limit  [Hz]
freq_upper = 4000 # PSD upper frequency limit [Hz]
freq_lower_mod = 10 # modulating signals's PSD lower frequency limit  [Hz]
freq_upper_mod = 20 # modulating signals's PSD upper frequency limit [Hz]

PSD = es.get_psd(freq, freq_lower, freq_upper) # one-sided flat-shaped PSD
PSD_modulating = es.get_psd(freq, freq_lower_mod, freq_upper_mod) # one-sided flat-shaped PSD

k_u = 5
x = es.nonstationary_signal(N,PSD,k_u=k_u, fs=fs, modulating_signal=('PSD',PSD_modulating))

sel = t<0.25
display(Audio(data=x, rate=fs))


fig, axs = plt.subplots(1,2)
fig.suptitle(f'Non-stationary signal: non-normal distribution in the range from {freq_lower} to {freq_upper} Hz')
plt.tight_layout(pad=1, w_pad=3)
axs[0].plot(t[sel], x[sel])
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel(f'Response in time (first 250 ms only)')
axs[1].psd(x, Fs=fs)
axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Power spectral density [dB/Hz]')
axs[1].set_xlim(0, 5000)
axs[1].set_ylim(-80, -20)
axs[1].set_yticks(np.arange(-80,-20,5))
plt.show()