This lecture builds directly on Lecture 4 (sampling and aliasing) and Lecture 5 (LTI systems, Bode plots, sampling as multiplication by an impulse train), so revisit those if the terms Nyquist frequency and frequency response are not yet second nature. The code examples use NumPy and SciPy from the course environment, and the case study reads the accelerometer records in data/SeismicResponse_noNoise.txt. Those records were shared with me for use in this course, not for public posting, so the file is not on this site: if you want to run the case study yourself, ask me or the TA for a copy.
Every code block in this chapter assumes that the setup below has run first, in the same Python session: it imports the libraries and defines the names the blocks share (see About the code in the Welcome chapter).
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
dB = lambda x: 20*np.log10(x)
fs_case = 200 # case-study sampling rate (Hz)
seis = np.loadtxt("data/SeismicResponse_noNoise.txt") # 5000 samples x 8 channels, a1 ... a8
t_case = np.arange(len(seis))/fs_caseLectures 4 and 5 kept asking “what frequencies does this signal contain?” — to pick a sampling rate, to place an anti-aliasing filter, to predict where an alias would land — without ever saying how one would find out. Lecture 4 even drew a spectrum (Figure 5.4) and asked you to take it on faith. This section builds the tool that produces such a picture: a way to take any measured signal apart into sinusoids, starting from what a signal is and what we actually want to know about it.
The same accelerometer trace can be summarized by a single number, by a spread, or by a full account of how it moves in time, and which of these you need decides which analysis you do. We classify signals first, then sort their characteristics into three questions — where the signal sits, how far it strays, and how it strays — of which only the last is sensitive to the time ordering of the samples.
A signal, for our purposes, is the physical information about a measured variable in transit: from the process to the sensor, between the stages of the measurement chain, or out of the chain as a record (Figliola and Beasley 2019). Its shape in time is its waveform, and the waveform carries two kinds of information — magnitude (how big) and frequency (how fast it changes) — which is why understanding waveforms comes before choosing an instrument or interpreting its output. Three classifications, two of which you already own from Lecture 4, sort the signals we meet:
Now the three questions, asked of a record \(y(t)\) observed from \(t_1\) to \(t_2\) — or of \(N\) samples of it.
(a) Where does it sit? The mean, or dc component,
\[ \bar{y} = \frac{1}{t_2 - t_1}\int_{t_1}^{t_2} y(t)\,\mathrm{d}t \qquad\text{or}\qquad \bar{y} = \frac{1}{N}\sum_{n=0}^{N-1} y[n], \tag{7.1}\]
is the static portion of the signal over the window. It says nothing about how much the signal moves.
(b) How much does it deviate from that static portion? The root-mean-square value answers with a single number,
\[ y_{\mathrm{rms}} = \sqrt{\frac{1}{t_2 - t_1}\int_{t_1}^{t_2} y(t)^2\,\mathrm{d}t} \qquad\text{or}\qquad y_{\mathrm{rms}} = \sqrt{\frac{1}{N}\sum_{n=0}^{N-1} y[n]^2}, \tag{7.2}\]
and it earns its definition from power: a current \(I(t)\) through a resistor dissipates \(\int I^2 R\,\mathrm{d}t\) (Lecture 3), so the constant current that would dissipate the same energy over the window is exactly \(I_{\mathrm{rms}}\). Once the mean has been subtracted, the RMS of what remains — the ac component — is the standard deviation of the record, the size of its fluctuations. The variance, the peak, the minimum, and the peak-to-peak range are answers of the same kind. Take Figliola and Beasley’s example, \(I(t) = 10\sin t\) A: over a half period (\(0\) to \(\pi\)) the mean is \(20/\pi \approx 6.37\) A and over a full period it is zero, while the RMS is \(\sqrt{50} \approx 7.07\) A over either — the power delivered does not care that the current reverses (Figliola and Beasley 2019).
(c) How does it deviate? This is the question the first two cannot touch. The steepest slope in the record, \(\max|\mathrm{d}y/\mathrm{d}t|\), is one answer; how strongly the value now predicts the value a second later is another; and “which frequencies is it made of?” is the answer the rest of this lecture develops. What these summaries share, and what separates them from (a) and (b), is that they depend on the order of the samples. Shuffle a record and its mean, RMS, peak, and every other answer to (a) and (b) are untouched, because they are computed from the collection of values without regard to when each occurred (Figure 7.1). Every answer to (c) is destroyed by the same shuffle. Static summaries are blind to time; dynamic ones are made of it.
seg = seis[1600:2600, 6] # channel a7, t = 8 to 13 s
rng = np.random.default_rng(12778)
shuf = rng.permutation(seg)
tt = np.arange(len(seg))/fs_case
fig, axs = plt.subplots(1, 2, figsize=(7.2, 2.8), sharey=True)
for ax, y, ttl in ((axs[0], seg, "as recorded"), (axs[1], shuf, "same samples, shuffled")):
m, r = y.mean(), y.std()
ax.plot(tt, y, lw=0.6, color="C0")
ax.axhline(m, color="k", lw=1); ax.axhspan(m - r, m + r, color="C1", alpha=0.25)
ax.set_title(f"{ttl}\nmean {m:.2f}, RMS about the mean {r:.2f}", fontsize=9)
ax.set_xlabel("time (s)")
axs[0].set_ylabel("acceleration $a_7$")
plt.tight_layout(); plt.show()
You are given two records of 1000 samples each: a vibration record from a footbridge, and the same record with its samples shuffled into a random order. Using only mean, RMS, minimum, maximum, and any other summary you can compute without looking at the plot, decide which is which. When you conclude that no such summary can tell them apart, name a quantity that could — and notice that every candidate you propose involves the order in which the samples arrived.
One key interesting fact that will support most of the rest of this section is that most signals (even non-deterministic ones) can be approximated via an infinite sum of periodic signals. This view allows us to estimate all of the above properties by analyzing the values of the periodic signals that comprise this approximation: the mean turns out to be the zero-frequency term, the RMS the total size of all the terms together, and the time-evolving behavior — question (c) — the distribution of that size across frequency. Figliola and Beasley’s analogy is a prism: white light in, a spectrum of pure colors out, each color a simple periodic function and its intensity an amplitude (Figliola and Beasley 2019). Fourier analysis is the prism for signals.
A rectified 60 Hz supply, a footbridge under a marching crowd, and mains hum in a strain-gauge line all repeat, and for a signal that repeats the question “which frequencies does it contain?” has an exact answer. A periodic signal satisfies \(y(t + T) = y(t)\) for some period \(T\), and Fourier’s claim is that any such signal of practical interest can be written as a constant plus a sum of sines and cosines whose frequencies are integer multiples of \(1/T\).
Imagine we have a periodic signal \(y(t)\) with period \(2\pi\), and suppose it can be written as the trigonometric series
\[ y(t) = A_0 + \sum_{n=1}^{\infty}\left(A_n \cos nt + B_n \sin nt\right). \tag{7.3}\]
The whole question is how to find the coefficients, and the answer rests on three integrals over one period — the orthogonality of sines and cosines: for positive integers \(m\) and \(n\),
\[ \int_{-\pi}^{\pi} \cos nt\,\cos mt\,\mathrm{d}t = \pi\,\delta_{nm},\qquad \int_{-\pi}^{\pi} \sin nt\,\sin mt\,\mathrm{d}t = \pi\,\delta_{nm},\qquad \int_{-\pi}^{\pi} \sin nt\,\cos mt\,\mathrm{d}t = 0, \tag{7.4}\]
where \(\delta_{nm}\) is \(1\) when \(n = m\) and \(0\) otherwise, and every one of these integrands also integrates to zero against a constant. Integrate Equation 7.3 over a period: every sine and cosine term vanishes and \(\int_{-\pi}^{\pi} y\,\mathrm{d}t = 2\pi A_0\). Multiply Equation 7.3 by \(\cos mt\) and integrate: by Equation 7.4 every term dies except the one with \(n = m\), which contributes \(\pi A_m\). The same with \(\sin mt\) picks out \(B_m\). Together these are the Euler formulas,
\[ A_0 = \frac{1}{2\pi}\int_{-\pi}^{\pi} y(t)\,\mathrm{d}t,\qquad A_n = \frac{1}{\pi}\int_{-\pi}^{\pi} y(t)\cos nt\,\mathrm{d}t,\qquad B_n = \frac{1}{\pi}\int_{-\pi}^{\pi} y(t)\sin nt\,\mathrm{d}t, \tag{7.5}\]
and the series with these coefficients is the Fourier series of \(y\) (Figliola and Beasley 2019). Each coefficient is a projection: how much of the signal “points along” a given sinusoid, in exactly the sense that a dot product measures how much one vector points along another.
Now imagine we make the period be an arbitrary \(T\). Nothing changes but the scale of time: with the fundamental circular frequency \(\omega = 2\pi/T\) the series and its coefficients become
\[ y(t) = A_0 + \sum_{n=1}^{\infty}\left(A_n \cos n\omega t + B_n \sin n\omega t\right), \tag{7.6}\]
\[ A_0 = \frac{1}{T}\int_{-T/2}^{T/2} y(t)\,\mathrm{d}t,\qquad A_n = \frac{2}{T}\int_{-T/2}^{T/2} y(t)\cos n\omega t\,\mathrm{d}t,\qquad B_n = \frac{2}{T}\int_{-T/2}^{T/2} y(t)\sin n\omega t\,\mathrm{d}t. \tag{7.7}\]
The \(n = 1\) term is the fundamental, at the lowest frequency in the series, \(f_1 = 1/T\); the terms for \(n = 2, 3, \ldots\) are the harmonics, at exact multiples of it. A sine and a cosine at the same frequency are one sinusoid with a phase, so the series is equally well written with one amplitude and one phase per harmonic,
\[ y(t) = A_0 + \sum_{n=1}^{\infty} C_n \cos\left(n\omega t - \phi_n\right),\qquad C_n = \sqrt{A_n^2 + B_n^2},\qquad \tan\phi_n = \frac{B_n}{A_n}, \tag{7.8}\]
and a plot of \(C_n\) against \(n f_1\) is the signal’s amplitude spectrum: the prism’s output.
A symmetry shortcut saves half the work. A function is even if \(g(-t) = g(t)\), like \(\cos\), and odd if \(h(-t) = -h(t)\), like \(\sin\). The product of an even function with an odd one is odd, and an odd function integrates to zero over a symmetric interval, so an even signal has \(B_n = 0\) for every \(n\) — a cosine series — and an odd signal has \(A_0 = A_n = 0\) — a sine series. A signal that is neither needs both. The rectified sine of a power supply is even, so its series is cosines only; the square wave below is odd, so it is sines only.
Take Figliola and Beasley’s square wave, \(y(t) = -5\) for \(-\pi < t < 0\) and \(+5\) for \(0 < t < \pi\), repeating with period \(2\pi\). It is odd, so only the \(B_n\) survive, and by Equation 7.5
\[ B_n = \frac{1}{\pi}\left[\int_{-\pi}^{0}(-5)\sin nt\,\mathrm{d}t + \int_{0}^{\pi} 5\sin nt\,\mathrm{d}t\right] = \frac{10}{n\pi}\left(1 - \cos n\pi\right) = \begin{cases} 20/(n\pi) & n \text{ odd} \\ 0 & n \text{ even,}\end{cases} \tag{7.9}\]
so that \(y(t) = \frac{20}{\pi}\left(\sin t + \tfrac{1}{3}\sin 3t + \tfrac{1}{5}\sin 5t + \cdots\right)\). Only odd harmonics, with amplitudes falling as \(1/n\): the spectrum of a square wave is a picket fence with every other picket missing. Figure 7.2 shows the partial sums closing in on the waveform — and one detail worth remembering when we get to leakage: the jump is the hard part. A discontinuity needs infinitely many harmonics, and every truncated series overshoots at the edge, however many terms are kept.
tq = np.linspace(-4, 4, 1600)
square = 5*np.sign(np.sin(tq))
partial = lambda t, K: sum(20/(n*np.pi)*np.sin(n*t) for n in range(1, 2*K, 2))
fig, (a1, a2) = plt.subplots(1, 2, figsize=(7.4, 3.0), gridspec_kw={"width_ratios": [1.6, 1]})
a1.plot(tq, square, color="0.6", lw=1, label="square wave")
for K, c in ((1, "C0"), (2, "C1"), (4, "C2"), (20, "C3")):
a1.plot(tq, partial(tq, K), color=c, lw=1.2, label=f"{K} term{'s' if K > 1 else ''}")
a1.set_xlabel("$t$"); a1.set_ylabel("$y(t)$"); a1.set_ylim(-7, 10); a1.legend(fontsize=7, ncol=3, loc="upper center")
nn = np.arange(1, 12)
Cn = np.where(nn % 2 == 1, 20/(nn*np.pi), 0)
a2.stem(nn, Cn, basefmt=" ")
a2.set_xlabel("harmonic $n$"); a2.set_ylabel("$C_n$"); a2.set_xticks(nn)
plt.tight_layout(); plt.show()
Build the square wave yourself: the slider sets how many odd harmonics are kept.
Now let’s imagine (as is the case in practice, mostly) that we don’t know the actual period of the signal we are measuring. The most that this period could be would be infinity — a signal that never repeats is one whose period is infinite — so the natural question is what happens to the series when \(T\) is taken to that limit. Letting the period grow without bound turns the sum of harmonics into an integral over every frequency, and the coefficients into a single complex function \(Y(\omega)\) whose magnitude says how much of each sinusoid is present and whose phase says where it sits in time.
Start from Equation 7.7 and pull the \(2/T\) outside: \(A_n = (2/T)\,A_T(n\omega)\) with \(A_T(\omega') = \int_{-T/2}^{T/2} y(t)\cos\omega' t\,\mathrm{d}t\), and likewise for \(B_n\). Since \(2/T = \omega/\pi\), the series Equation 7.6 reads
\[ y(t) = A_0 + \frac{1}{\pi}\sum_{n=1}^{\infty}\left[A_T(n\omega)\cos n\omega t + B_T(n\omega)\sin n\omega t\right]\omega . \]
As \(T \to \infty\) three things happen at once: the harmonic spacing \(\omega = 2\pi/T\) shrinks toward an infinitesimal \(\mathrm{d}\omega\); the harmonic frequencies \(n\omega\) fill in every value of a continuous \(\omega\); and the sum, which has the form \(\sum g(n\omega)\,\omega\), becomes the integral \(\int g(\omega)\,\mathrm{d}\omega\). (The constant \(A_0\) goes to zero for any pulse of finite energy, since it is a finite integral divided by \(T\).) The Fourier coefficients \(A_n\) and \(B_n\) have become continuous functions of frequency,
\[ A(\omega) = \int_{-\infty}^{\infty} y(t)\cos\omega t\,\mathrm{d}t,\qquad B(\omega) = \int_{-\infty}^{\infty} y(t)\sin\omega t\,\mathrm{d}t, \tag{7.10}\]
and the signal is recovered from them as \(y(t) = \frac{1}{\pi}\int_0^{\infty}\left[A(\omega)\cos\omega t + B(\omega)\sin\omega t\right]\mathrm{d}\omega\). Two real functions are one complex function, and Euler’s identity \(e^{-j\theta} = \cos\theta - j\sin\theta\) (with \(j = \sqrt{-1}\), the engineer’s letter, since \(i\) is a current) packs the pair into the Fourier transform
\[ Y(\omega) \equiv A(\omega) - jB(\omega) = \int_{-\infty}^{\infty} y(t)\,e^{-j\omega t}\,\mathrm{d}t . \tag{7.11}\]
For a real signal \(Y(-\omega)\) is the complex conjugate of \(Y(\omega)\), which lets the recovery integral run over negative frequencies too and absorb the sine and cosine into one exponential. Written in hertz, with \(\omega = 2\pi f\), the pair is
\[ Y(f) = \int_{-\infty}^{\infty} y(t)\,e^{-j2\pi f t}\,\mathrm{d}t, \qquad y(t) = \int_{-\infty}^{\infty} Y(f)\,e^{+j2\pi f t}\,\mathrm{d}f, \tag{7.12}\]
the forward transform and the inverse Fourier transform, and with \(f\) as the variable no stray factors of \(2\pi\) appear anywhere. The transform is two-sided — it runs over negative frequencies as well — and for a real signal the negative half is the mirror image of the positive half, so it carries no new information; plots show \(f \geq 0\).
What is \(Y(f)\)? At each frequency it is a complex number, and its polar form is the useful one:
\[ Y(f) = |Y(f)|\,e^{j\phi(f)},\qquad |Y(f)| = \sqrt{A(f)^2 + B(f)^2},\qquad \phi(f) = -\tan^{-1}\frac{B(f)}{A(f)} . \tag{7.13}\]
The magnitude \(|Y(f)|\) is the amplitude spectrum of Equation 7.8 with the pickets pushed together into a continuous curve: how much of the sinusoid at \(f\) the signal contains. The phase \(\phi(f)\) says where that sinusoid sits in time — shift the whole signal later by \(t_0\) and every \(|Y(f)|\) is unchanged while every phase turns by \(-2\pi f t_0\). A signal is, in other words, described exactly as well by \(Y(f)\) as by \(y(t)\); the transform is a change of point of view, not a loss of information, and it is the point of view from which Lecture 5’s filters and this lecture’s structures are easy to read.
Figure 7.3 shows the limit at work on a rectangular pulse. Repeated with period \(T\), the pulse has a line spectrum whose pickets sit \(1/T\) apart and whose heights, scaled by \(T\), all lie on one curve. Double the period and the pickets crowd twice as densely under the same curve; the curve is \(|Y(f)|\), the transform of the single pulse, and the pickets are its samples.
tau = 1.0
fig, axs = plt.subplots(3, 2, figsize=(7.2, 5.0), gridspec_kw={"width_ratios": [1, 1.3]})
ff = np.linspace(0, 4, 800)
Ycont = tau*np.abs(np.sinc(ff*tau)) # |Y(f)| of the single pulse (numpy's sinc is sin(pi x)/(pi x))
for row, T in zip(axs, (2, 4, 8)):
tt = np.linspace(-8, 8, 3200)
pulse = (np.abs((tt + T/2) % T - T/2) <= tau/2).astype(float)
row[0].plot(tt, pulse, lw=1); row[0].set_yticks([0, 1]); row[0].set_ylabel(f"$T = {T}$")
nn = np.arange(1, int(4*T) + 1)
An_ = (2/T)*tau*np.sinc(nn*tau/T) # cosine coefficients of the even pulse train
row[1].plot(ff, Ycont, color="0.6", lw=1)
row[1].stem(nn/T, T*np.abs(An_)/2, basefmt=" ", markerfmt="C0.", linefmt="C0-")
row[1].stem([0], [T*tau/T], basefmt=" ", markerfmt="C0.", linefmt="C0-")
row[1].set_ylim(0, 1.1); row[1].set_yticks([0, 1])
axs[-1][0].set_xlabel("$t$"); axs[-1][1].set_xlabel("$f$")
axs[0][1].set_title("$T\\,C_n / 2$ against $n/T$, with $|Y(f)|$", fontsize=9)
plt.tight_layout(); plt.show()
A periodic signal has period \(T = 0.5\) s. Where are its harmonics on the frequency axis, and what happens to that picture if the period is doubled, then doubled again, then made infinite?
The harmonics sit at multiples of \(1/T = 2\) Hz: \(2, 4, 6, \ldots\) Hz. Doubling \(T\) to 1 s puts them 1 Hz apart, and to 2 s puts them 0.5 Hz apart — the pickets crowd together, and (scaled by \(T\)) they keep lying on the same curve. In the limit of infinite period the pickets are dense on the axis and the curve itself, \(|Y(f)|\), is the description: the Fourier transform is not a different tool from the series but its limit, which is why a signal that does not repeat still has a spectrum.
In practice, we will be computing this transform on a discrete signal \(y[n]\) given that we will apply it after the ADC has digitized the signals from the sensors. Thus, we need a discrete version of the transform — and the version we get carries a fine print about periodicity that is the source of every practical difficulty in the rest of the lecture.
Start from the forward transform in Equation 7.12 and make it computable. The record is \(N\) samples, \(y[n] = y(n\,\delta t)\) for \(n = 0, \ldots, N-1\), so the integral becomes a sum with \(\mathrm{d}t \to \delta t\) and \(t \to n\,\delta t\). The record is finite, so the transform can only be asked about frequencies that fit a whole number of cycles into it — \(f_k = k/(N\,\delta t)\), or equivalently \(k\,f_s/N\), for integer \(k\). Substituting both, \(e^{-j2\pi f_k n\,\delta t} = e^{-j2\pi kn/N}\), and dropping the overall factor \(\delta t\) (a convention; it only scales the result), we have the discrete Fourier transform
\[ Y[k] = \sum_{n=0}^{N-1} y[n]\,e^{-j2\pi kn/N},\qquad k = 0, 1, \ldots, N-1, \tag{7.14}\]
with its inverse
\[ y[n] = \frac{1}{N}\sum_{k=0}^{N-1} Y[k]\,e^{+j2\pi kn/N},\qquad n = 0, 1, \ldots, N-1 . \tag{7.15}\]
\(N\) numbers in, \(N\) complex numbers out, and Equation 7.15 says nothing was lost: the \(N\) values \(Y[k]\) are just another set of coordinates for the same record, exactly as \(Y(f)\) was for \(y(t)\). This is the convention NumPy’s fft uses, and the one we adopt; Figliola and Beasley’s version divides by \(N/2\) so that the result reads directly in amplitude, which is the scaling in Equation 7.17 below (Figliola and Beasley 2019).
What did this do to our assumption of an infinite period? It threw it away. The exponential \(e^{-j2\pi kn/N}\) repeats every \(N\) samples, so the only signals Equation 7.15 can build are ones that repeat every \(N\) samples — the DFT treats your \(N\) samples as one period of a periodic signal of period \(N\,\delta t\), whether or not the physical signal ever repeated. This is the fine print. It is also the reason the frequencies come out where they do: the lowest nonzero frequency the DFT can distinguish is one full cycle per record,
\[ \delta f = \frac{1}{N\,\delta t} = \frac{f_s}{N}, \tag{7.16}\]
Lecture 4’s frequency resolution (Equation 5.5), now with its origin visible. Resolution is bought with record length, \(N\,\delta t\), and with nothing else: sampling faster with the same \(N\) shortens the record and makes \(\delta f\) coarser, not finer (Figure 7.4). Three of the \(N\) bins deserve names. Bin \(k = 0\) is \(Y[0] = \sum y[n] = N\bar{y}\), the mean of Equation 7.1 scaled by \(N\) — the static component, question (a), is the first coefficient of the spectrum. Bin \(k = N/2\) (for even \(N\)) sits at \(f_s/2\), the Nyquist frequency. And the bins above \(N/2\) are the negative frequencies in disguise: for a real record \(Y[N-k]\) is the complex conjugate of \(Y[k]\), the mirror image Lecture 5 saw in the replica picture, so the upper half repeats the lower half and NumPy’s rfft simply omits it. Reading an amplitude off the lower half requires the scaling
\[ C_k = \frac{2}{N}\,|Y[k]|,\qquad 0 < k < \frac{N}{2}, \tag{7.17}\]
the factor \(2\) collecting the mirror bin’s share and the \(1/N\) undoing the sum’s growth with record length — this is the amplitude \(C_n\) of Equation 7.8, for the harmonic that fits \(k\) cycles into the record.
N_ = 16
fig, (a1, a2) = plt.subplots(2, 1, figsize=(6.6, 3.2))
n_ = np.arange(N_)
a1.stem(n_, 0.6 + 0.4*np.sin(2*np.pi*3*n_/N_), basefmt=" ", markerfmt="C0.", linefmt="C0-")
a1.annotate("", xy=(1, 1.25), xytext=(0, 1.25), arrowprops=dict(arrowstyle="<->", lw=0.8)); a1.text(0.5, 1.32, r"$\delta t = 1/f_s$", ha="center", fontsize=8)
a1.annotate("", xy=(N_ - 1, -0.25), xytext=(0, -0.25), arrowprops=dict(arrowstyle="<->", lw=0.8)); a1.text(N_/2, -0.5, r"record length $N\,\delta t$", ha="center", fontsize=8)
a1.set_ylim(-0.7, 1.5); a1.set_yticks([]); a1.set_xlabel("sample $n$", fontsize=8); a1.set_title("time domain: $N$ samples", fontsize=9)
k_ = np.arange(N_//2 + 1)
a2.stem(k_, np.exp(-((k_ - 3)/0.5)**2), basefmt=" ", markerfmt="C1.", linefmt="C1-")
a2.annotate("", xy=(1, 1.15), xytext=(0, 1.15), arrowprops=dict(arrowstyle="<->", lw=0.8)); a2.text(0.5, 1.22, r"$\delta f = f_s/N$", ha="center", fontsize=8)
a2.axvline(N_//2, color="k", ls="--", lw=0.8); a2.text(N_//2 - 0.3, 0.6, r"$f_s/2$ at $k = N/2$", ha="right", fontsize=8)
a2.set_ylim(0, 1.4); a2.set_yticks([]); a2.set_xlabel("bin $k$", fontsize=8); a2.set_title("frequency domain: $N/2 + 1$ distinct bins", fontsize=9)
plt.tight_layout(); plt.show()
You keep \(N\) fixed at 1024 samples and double the sampling rate from 100 to 200 Hz. Does the frequency resolution of your DFT improve?
No — it gets worse. By Equation 7.16, \(\delta f = f_s/N\) goes from about 0.098 Hz to 0.195 Hz, because the record now lasts 5.12 s instead of 10.24 s and the lowest frequency that fits one cycle into it has doubled. What the higher rate bought is range: the Nyquist frequency moved from 50 to 100 Hz. Resolution comes from record length, range from sampling rate, and only more samples can buy both.
A call to numpy.fft.fft returns a list of \(N\) complex numbers indexed from 0 to \(N-1\), and nothing in that list says hertz. The conversion is the one we just derived, \(f_k = k\,f_s/N\), and a worked example makes it concrete: a sinusoid of 120 Hz and unit amplitude, sampled at 360 Hz. (120 Hz is the fundamental of a full-wave-rectified 60 Hz supply and the flicker frequency of mains-powered lamps — it turns up uninvited in a great many lab records, including the one Lecture 5’s exercises aliased to DC.)
Three samples per cycle, since \(360/120 = 3\), so the sampled tone takes only three distinct values, \(\sin(0)\), \(\sin(2\pi/3)\), and \(\sin(4\pi/3)\) — that is, \(0\), \(0.866\), and \(-0.866\) — over and over. For the tone to land on a bin the record must hold a whole number of cycles, so \(N\) must be a multiple of 3; take \(N = 36\), a tenth of a second. Then \(\delta f = 360/36 = 10\) Hz and the tone falls at \(k = 120/10 = 12\), exactly. Figure 7.5 shows the result: \(|Y[12]| = 18 = N/2\), its mirror at \(k = 24\) (that is, \(N - 12\)), and nothing anywhere else. The amplitude scaling of Equation 7.17 gives \(C_{12} = 2 \times 18/36 = 1\), the amplitude we started with.
fs120, f120, N120 = 360, 120, 36
n120 = np.arange(N120)
y120 = np.sin(2*np.pi*f120*n120/fs120)
Y120 = np.fft.fft(y120)
fig, (a1, a2) = plt.subplots(1, 2, figsize=(7.2, 2.9), gridspec_kw={"width_ratios": [1, 1.4]})
a1.stem(n120/fs120*1e3, y120, basefmt=" ", markerfmt="C0.", linefmt="C0-")
a1.plot(np.linspace(0, N120/fs120, 600)*1e3, np.sin(2*np.pi*f120*np.linspace(0, N120/fs120, 600)), color="0.8", lw=0.8, zorder=0)
a1.set_xlabel("time (ms)"); a1.set_ylabel("$y[n]$")
a2.stem(n120, np.abs(Y120), basefmt=" ", markerfmt="C1.", linefmt="C1-")
a2.axvspan(N120/2, N120 - 0.5, color="0.9", zorder=0)
a2.set_xlabel("bin $k$"); a2.set_ylabel("$|Y[k]|$")
a2.annotate("$k = 12$: $|Y| = N/2 = 18$", (12, 18), xytext=(2, 15.5), fontsize=8)
a2.annotate("mirror, $k = 24$", (24, 18), xytext=(25, 15.5), fontsize=8)
top = a2.secondary_xaxis("top", functions=(lambda k: k*fs120/N120, lambda f: f*N120/fs120))
top.set_xlabel("$f_k$ (Hz)", fontsize=8)
plt.tight_layout(); plt.show()
Computing Equation 7.14 literally costs \(N\) multiplications for each of \(N\) bins, \(N^2\) in all, which for a one-second record at 200 kS/s is \(4 \times 10^{10}\) operations — too slow to be routine in 1965 and too slow for real-time instruments today. The fast Fourier transform (FFT) is an algorithm, published by Cooley and Tukey, that computes exactly the same \(N\) numbers in about \(N \log_2 N\) operations by splitting the sum recursively into even- and odd-indexed halves (Cooley and Tukey 1965). For \(N = 2^{20}\) (about a million samples) the saving is a factor of \(N/\log_2 N \approx 52{,}000\), which is the difference between minutes and milliseconds — and the reason spectral analysis went from a specialist’s calculation to a function call. Its one habit worth knowing is a preference for \(N\) that is a power of 2, which is why records are often padded or truncated to \(1024\), \(4096\), and so on. The FFT is not a different transform: numpy.fft.fft returns Equation 7.14 to within rounding error, for any \(N\).
Since our assumption after the DFT is that the signals are periodic, we need to be careful about how we interpret the results if we are only carrying out the analysis on a subset (a window) of the full signal we’re measuring. Unless the record holds a whole number of cycles of every component in it, its periodic extension has a jump at the seam, and that jump smears energy into neighboring bins — leakage — so that a pure tone reads as a smaller peak with a skirt. Tapering the record with a window softens the seam, trading a little resolution for much less smearing.
Figure 7.6 draws the assumption the DFT makes without asking. The signal continues before and after the record (top); we captured a finite piece of it (middle); and Equation 7.15 describes not that piece but the piece repeated end to end forever (bottom). If the record held a whole number of cycles the copies would join seamlessly and the extension would be the original signal. Here it does not, and at every seam the extended signal jumps from the last sample’s value to the first’s — a discontinuity the physical signal never had.
fpe, fs_pe, Tpe = 2.5, 100, 1.0
tt = np.linspace(-1.5, 2.5, 2000)
tr = np.arange(0, Tpe, 1/fs_pe)
fig, axs = plt.subplots(3, 1, figsize=(6.8, 4.6), sharex=True)
axs[0].plot(tt, np.cos(2*np.pi*fpe*tt), lw=1); axs[0].set_title("the actual signal", fontsize=9)
axs[1].plot(tr, np.cos(2*np.pi*fpe*tr), lw=1); axs[1].set_xlim(-1.5, 2.5); axs[1].set_title("the captured record ($N\\,\\delta t = 1$ s)", fontsize=9)
ext = np.cos(2*np.pi*fpe*(tt % Tpe))
axs[2].plot(tt, ext, lw=1, color="C3"); axs[2].set_title("what the DFT analyzes: the record repeated", fontsize=9)
from matplotlib.patches import Ellipse
for s in (-1, 0, 1, 2):
axs[2].add_patch(Ellipse((s, 0), width=0.16, height=2.5, fill=False, color="k", lw=0.8))
for ax in axs: ax.set_yticks([-1, 0, 1]); ax.set_ylim(-1.35, 1.35); ax.axvspan(0, Tpe, color="0.93", zorder=0)
axs[2].set_xlabel("time (s)")
plt.tight_layout(); plt.show()
Why should a jump at the seam trouble the spectrum? Because the square wave already told us: a discontinuity is the most expensive feature a Fourier series can represent, needing every harmonic with amplitudes that decay only as \(1/n\) (Equation 7.9). The DFT of the extended signal must therefore spend energy on harmonics the true signal never contained, and it takes that energy from the tone’s own bin. The tone’s peak comes out lower than the truth and a skirt of spurious content spreads to either side — content that leaked out of the bin it belonged in. The same mechanism is responsible for Lecture 4’s amplitude ambiguity (Figure 5.4), which we can now state precisely: the DFT reports a tone’s amplitude correctly only when the tone fits a whole number of cycles into the record.
Return to the 120 Hz tone, and change one thing: sample \(N = 40\) points instead of 36. Now \(\delta f = 9\) Hz, and 120 Hz sits between bins 13 (117 Hz) and 14 (126 Hz) — the record holds \(13\tfrac{1}{3}\) cycles. Figure 7.7 shows the damage: the largest bin reads an amplitude of \(0.80\) instead of \(1\), its neighbor reads \(0.44\), and every bin above dc is nonzero. About 11% of the tone’s energy has left the immediate neighborhood of its true frequency.
The remedy is to soften the seam before transforming. A window \(w[n]\) is a sequence that is \(1\) in the middle of the record and tapers smoothly to \(0\) at both ends; multiply the record by it,
\[ y_w[n] = w[n]\,y[n], \tag{7.18}\]
and the periodic extension of \(y_w\) now joins gently at every seam, because both ends have been brought to zero. The unwindowed DFT is itself the windowed DFT with the rectangular window — all ones — so the choice is never whether to window, only which window. The workhorse is the Hann window, \(w[n] = \tfrac{1}{2}\left[1 - \cos\!\big(2\pi n/(N-1)\big)\right]\), a raised cosine. Figure 7.7 shows what it does to the \(N = 40\) record: the skirt collapses by an order of magnitude and the leaked energy outside the tone’s neighborhood falls from 11% to under 1%, at the price of a peak that is now two bins wide rather than one. (The window also scales the record — it averages \(\tfrac{1}{2}\) — so amplitudes read from a Hann-windowed spectrum are multiplied by 2, or in general divided by the window’s mean.)
def amp_spectrum(y, w):
A = np.abs(np.fft.rfft(y*w))/w.sum() # corrected for the window's mean
A[1:(len(y) + 1)//2] *= 2 # fold in the mirror bins, 0 < k < N/2 (@eq-dft-amp); dc and Nyquist have none
return A
cases = [("$N = 36$, rectangular", 36, "rect"), ("$N = 40$, rectangular", 40, "rect"), ("$N = 40$, Hann", 40, "hann")]
fig, axs = plt.subplots(1, 3, figsize=(7.6, 2.9), sharey=True)
for ax, (ttl, N_, kind) in zip(axs, cases):
n_ = np.arange(N_); y_ = np.sin(2*np.pi*f120*n_/fs120)
w_ = np.ones(N_) if kind == "rect" else np.hanning(N_)
A_ = amp_spectrum(y_, w_); fk = np.fft.rfftfreq(N_, 1/fs120)
ax.stem(fk, np.maximum(dB(A_ + 1e-12), -60), basefmt=" ", markerfmt="C1.", linefmt="C1-")
ax.axvline(120, color="0.6", ls=":", lw=0.8); ax.set_title(ttl, fontsize=9); ax.set_xlabel("$f_k$ (Hz)")
axs[0].set_ylabel("amplitude (dB re 1)"); axs[0].set_ylim(-62, 3)
plt.tight_layout(); plt.show()
What kind of operation is windowing? It is linear — window a sum and you get the sum of the windowed records — but it is not time-invariant: slide the record by ten samples and the taper falls on different samples, so windowing is a time-varying gain and not one of Lecture 5’s LTI blocks (Section 6.2.1). Its natural description is the mirror image of a filter’s. A filter convolves in time and multiplies in frequency; a window multiplies in time, Equation 7.18, and therefore convolves in frequency,
\[ Y_w(f) = \left(Y * W\right)(f) = \int_{-\infty}^{\infty} Y(\nu)\,W(f - \nu)\,\mathrm{d}\nu, \tag{7.19}\]
where \(W\) is the window’s own transform. Every line in the true spectrum is replaced by a copy of \(W\) centered on it, so the spectrum you measure is the true spectrum seen through the window’s transform, and choosing a window means choosing the shape of \(W\) (Figure 7.8). The rectangular window’s \(W\) has a narrow main lobe — good resolution — but side lobes only 13 dB down that decay slowly, which is the skirt in Figure 7.7. The Hann window’s main lobe is twice as wide, and its first side lobe is 31 dB down and the rest fall away quickly. That is the whole trade: a narrow main lobe resolves close tones, low side lobes stop a strong tone from burying a weak one, and no window has both, which is why the catalog of windows runs to dozens (Harris 1978). For finding modal peaks in vibration data, the Hann window is the default and is rarely wrong.
Nw, pad = 64, 64*32
n_ = np.arange(Nw)
fig, (a1, a2) = plt.subplots(1, 2, figsize=(7.2, 2.9), gridspec_kw={"width_ratios": [1, 1.5]})
for w_, name, c in ((np.ones(Nw), "rectangular", "C0"), (np.hanning(Nw), "Hann", "C1")):
a1.plot(n_, w_, color=c, lw=1.4, label=name)
W_ = np.abs(np.fft.fft(w_, pad))/w_.sum()
fbins = np.fft.fftfreq(pad, 1/Nw) # frequency in bins
order = np.argsort(fbins)
a2.plot(fbins[order], dB(W_[order] + 1e-12), color=c, lw=1.2, label=name)
a1.set_xlabel("sample $n$"); a1.set_ylabel("$w[n]$"); a1.legend(fontsize=8)
a2.set_xlim(-10, 10); a2.set_ylim(-80, 3); a2.set_xlabel("frequency (bins)"); a2.set_ylabel("$|W|$ (dB)")
a2.axhline(-13.3, color="C0", ls=":", lw=0.8); a2.axhline(-31.5, color="C1", ls=":", lw=0.8); a2.legend(fontsize=8)
plt.tight_layout(); plt.show()
Try it: the 120 Hz tone at \(f_s = 360\) Hz again, with the record length and the window under your control. Watch what happens to the peak and to the skirt as \(N\) passes through multiples of 3.
Set the widget above (or the code in Figure 7.7) aside for a moment. The 120 Hz tone is sampled at 360 Hz for \(N = 42\) samples, then for \(N = 44\). For each, say before computing whether the tone lands on a bin, which bin or pair of bins it falls nearest, and what amplitude the largest bin will read — then check. Explain, in one sentence each, what \(N = 42\) has in common with \(N = 36\), and why \(N = 44\) leaks about as badly as \(N = 40\) despite being the longer record.
A colleague says that applying a Hann window “removes the leakage” from a spectrum. Is that right?
Not quite: it redistributes it. By Equation 7.19 every spectral line is convolved with the window’s transform, and the Hann window’s transform has a main lobe twice as wide as the rectangular window’s. The far skirt drops by tens of decibels, but the energy near the tone spreads over two bins instead of one — resolution is traded for side-lobe suppression. The only way to have no leakage at all is a record holding a whole number of cycles of every component, which for real data with unknown frequencies never happens.
This section is optional: it is not among the lecture’s learning objectives, and nothing after it depends on it. Read it when you meet a signal — wind, traffic, machinery — that never settles into a repeating pattern and you need a spectrum of the process rather than of one record.
For a nondeterministic signal — traffic-induced vibration of a bridge deck, wind pressure on a façade — the DFT magnitude of a single record never settles down no matter how long the record, so we need a spectrum that describes the process rather than one realization of it. The power spectral density tells how the signal’s variance is distributed over frequency; it is estimated by averaging the periodograms of windowed segments (Welch’s method), and its integral is the variance — the RMS from Section 7.2, now resolved by frequency.
Everything so far took the signal to be deterministic: a tone has an amplitude, and a longer record pins it down more finely. A random signal has no amplitude at any frequency; it has a distribution of amplitudes, and what is fixed — for a stationary process, one whose statistics do not drift over the record — is the average power at each frequency. The natural estimator is the squared DFT, scaled so that it integrates to the variance: the periodogram
\[ \hat{S}[k] = \frac{2}{N f_s}\,\left|Y[k]\right|^2,\qquad 0 < k < \frac{N}{2}, \tag{7.20}\]
with units of (signal units)\(^2\) per hertz — acceleration in \((\mathrm{m/s^2})^2/\mathrm{Hz}\), pressure in \(\mathrm{Pa^2/Hz}\) — so that the area under it between two frequencies is the variance contributed by that band. (The factor 2 folds the mirror half in, as in Equation 7.17; the dc bin and the Nyquist bin get half of it.) That the area is the variance is Parseval’s theorem,
\[ \sum_{n=0}^{N-1} y[n]^2 = \frac{1}{N}\sum_{k=0}^{N-1}\left|Y[k]\right|^2, \tag{7.21}\]
energy in the time domain equals energy in the frequency domain: for a zero-mean record, dividing both sides by \(N\) makes the left the variance, \(y_{\mathrm{rms}}^2\) of Equation 7.2, and the right the sum of the periodogram times \(\delta f\). So the RMS is the square root of the area under the PSD — question (b) again — and the PSD is what question (c) adds: where on the frequency axis that area sits.
The periodogram has one crippling defect. Each of its bins is, for a random signal, an estimate from what amounts to one sample, and its scatter is as large as its mean regardless of \(N\): a longer record gives more bins, each as noisy as before. The left panel of Figure 7.9 shows it — a spectrum that looks like grass. Welch’s remedy is to spend the record on averaging rather than on resolution (Welch 1967). Split the \(N\) samples into \(L\) segments of length \(M\) (overlapping by half, which recovers some of what the window taper throws away), window each segment, compute its periodogram, and average:
\[ \hat{S}_W[k] = \frac{1}{L}\sum_{\ell=1}^{L}\hat{S}_\ell[k],\qquad \delta f = \frac{f_s}{M}. \tag{7.22}\]
The scatter of the average falls as \(1/\sqrt{L}\), and the price is resolution: \(\delta f\) is now \(f_s/M\), not \(f_s/N\). Choosing \(M\) is choosing where to stand between a spectrum too noisy to read and one too coarse to resolve the peaks you care about. SciPy’s signal.welch implements exactly this, with a Hann window and 50% overlap by default; you supply \(M\) as nperseg.
The example in Figure 7.9 is a building mode: a lightly damped oscillator at 3 Hz driven by broadband (white) excitation, sampled at 100 Hz for 200 s. Its true PSD is known — it is the excitation’s flat spectrum shaped by the oscillator’s frequency response, a peak at 3 Hz whose half-power width is \(2\zeta f_n = 0.3\) Hz for 5% damping. The periodogram scatters around it wildly; the Welch estimate with \(M = 1024\) (\(\delta f \approx 0.1\) Hz, 38 segments) follows it closely and still resolves the peak. And the areas agree: the Welch estimate integrates to the record’s variance to within 1%.
fs_ex, fn_ex, zeta_ex = 100, 3.0, 0.05
wn = 2*np.pi*fn_ex
b_ex, a_ex = signal.bilinear([wn**2], [1, 2*zeta_ex*wn, wn**2], fs_ex) # the oscillator as a digital filter
rng_ex = np.random.default_rng(1)
x_ex = signal.lfilter(b_ex, a_ex, rng_ex.normal(0, 1, 20000)) # 200 s of response to unit white noise
f_p, S_p = signal.periodogram(x_ex, fs=fs_ex)
f_w, S_w = signal.welch(x_ex, fs=fs_ex, nperseg=1024)
f_t = np.linspace(0.01, 50, 5000); _, H_t = signal.freqz(b_ex, a_ex, worN=f_t, fs=fs_ex)
S_t = 2*np.abs(H_t)**2/fs_ex # true one-sided PSD: white input shaped by |H|^2
fig, (a1, a2) = plt.subplots(1, 2, figsize=(7.4, 3.0), sharey=True)
a1.semilogy(f_p, S_p, color="0.7", lw=0.5, label="periodogram"); a1.semilogy(f_t, S_t, "k", lw=1, label="true PSD")
a2.semilogy(f_w, S_w, "C1", lw=1.4, label="Welch, $M = 1024$"); a2.semilogy(f_t, S_t, "k", lw=1, label="true PSD")
for ax in (a1, a2): ax.set_xlim(0, 20); ax.set_xlabel("frequency (Hz)"); ax.legend(fontsize=8)
a1.set_ylabel("PSD (units$^2$/Hz)"); a1.set_ylim(1e-6, 1e2)
plt.tight_layout(); plt.show()
Two habits keep PSDs honest. First, subtract the mean before estimating: a dc offset is a spike at \(f = 0\) whose leakage can bury the low-frequency content you want. Second, read the units: a PSD in \((\mathrm{m/s^2})^2/\mathrm{Hz}\) is not an amplitude, and comparing “peak heights” across records with different \(\delta f\) is meaningless unless both are densities. Integrate over a band and take the square root, and you have an RMS in \(\mathrm{m/s^2}\) that can be compared with anything.
The PSD of a floor’s vertical acceleration is flat at \(2 \times 10^{-4}\ (\mathrm{m/s^2})^2/\mathrm{Hz}\) from 1 to 4 Hz and at \(5 \times 10^{-5}\ (\mathrm{m/s^2})^2/\mathrm{Hz}\) from 4 to 10 Hz, and negligible elsewhere. The specification limits the RMS acceleration in the 1–10 Hz band to \(0.05\ \mathrm{m/s^2}\). Does the floor comply?
The variance is the area: \(2 \times 10^{-4} \times 3 + 5 \times 10^{-5} \times 6 = 6 \times 10^{-4} + 3 \times 10^{-4} = 9 \times 10^{-4}\ (\mathrm{m/s^2})^2\), so the RMS is \(\sqrt{9 \times 10^{-4}} = 0.030\ \mathrm{m/s^2}\), below the limit. Note that two-thirds of the variance comes from the lower band even though it is half as wide — density times bandwidth, not density alone, is what counts.
Lecture 5 drew the Bode plot of an RC filter but never said how to compute what that filter does to an entire measured signal, only to one sinusoid at a time. Now that we can compute \(Y(\omega)\) for any signal \(y(t)\), every block in the measurement chain becomes a multiplication in the frequency domain, which makes it easy to reason about how the chain affects the magnitude and phase of each frequency component throughout.
Recall from Section 6.2.1 that an LTI block is described by its impulse response \(h(t)\), and that its output for any input is the convolution \(y_{\mathrm{out}}(t) = \int h(\tau)\,y(t - \tau)\,\mathrm{d}\tau\) — every input sample launches a scaled, delayed copy of \(h\), and the output is the sum of the copies. Take the transform of both sides. Substituting the convolution into Equation 7.12, swapping the order of integration, and changing variable \(t \to t + \tau\) in the inner integral gives \(\int h(\tau)e^{-j2\pi f\tau}\,\mathrm{d}\tau\) times \(\int y(t)e^{-j2\pi f t}\,\mathrm{d}t\), that is,
\[ Y_{\mathrm{out}}(f) = H(f)\,Y(f), \tag{7.23}\]
the convolution theorem: convolution in time is multiplication in frequency, and \(H(f)\) — the transform of the impulse response — is exactly the frequency response of Lecture 5, whose Bode plot is \(|H|\) and \(\angle H\). This one line is why the frequency domain exists as an engineering tool. It explains the sinusoid’s privilege (a sinusoid is one line in \(Y\), and multiplying one line by \(H(f)\) can only scale and shift it, never change its shape); it reads a Bode plot as an instruction (“multiply the spectrum by this”); and it is the dual of Equation 7.19, with the roles of time and frequency exchanged. In polar form the multiplication separates: magnitudes multiply and phases add, so for two blocks in cascade
\[ \left|H_{\mathrm{chain}}\right| = \left|H_1\right|\left|H_2\right|,\qquad \angle H_{\mathrm{chain}} = \angle H_1 + \angle H_2, \tag{7.24}\]
and in decibels the magnitudes add — which is the reason the decibel was invented. Figure 7.10 draws the measurement chain both ways.
The sensor itself is the first block, and the ADXL335 makes the point with a single capacitor. Its outputs leave through a \(32\ \mathrm{k\Omega}\) resistor, and the bandwidth capacitor \(C\) the designer adds at each output pin forms with it a first-order RC low-pass — Lecture 5’s anti-aliasing filter (Section 6.4) — with
\[ H(f) = \frac{1}{1 + j f/f_c},\qquad f_c = \frac{1}{2\pi \cdot 32\ \mathrm{k\Omega} \cdot C}, \tag{7.25}\]
so that \(|H| = 1/\sqrt{1 + (f/f_c)^2}\) and \(\angle H = -\tan^{-1}(f/f_c)\). With \(C = 0.1\ \mu\mathrm{F}\), \(f_c \approx 50\) Hz, as in Lecture 5. With \(C = 1\ \mu\mathrm{F}\), \(f_c \approx 5\) Hz, and Figure 7.11 shows what Equation 7.23 then does to a structure with modes at 2.5 and 25 Hz: the first mode passes at 89% of its amplitude and lagging by 27°, the second arrives at 20% (−14 dB) — attenuated before it ever reaches the ADC, and, if nobody looked at the capacitor, silently absent from the spectrum. The lesson is not that filtering is bad; it is that a modal peak that is missing from a spectrum may be missing because of the chain, not the structure, and Equation 7.23 is how you check.
fc_1u = 1/(2*np.pi*32e3*1e-6)
fs_v, Tv = 200, 20
tv = np.arange(0, Tv, 1/fs_v)
xv = np.sin(2*np.pi*2.5*tv) + np.sin(2*np.pi*25*tv)
fv = np.fft.rfftfreq(len(tv), 1/fs_v)
Xv = amp_spectrum(xv, np.hanning(len(tv)))
Hv = 1/(1 + 1j*fv/fc_1u)
fig, ax = plt.subplots(figsize=(6.6, 3.0))
ax.plot(fv, Xv, color="0.7", lw=1.2, label="input spectrum $|Y|$")
ax.plot(fv, np.abs(Hv)*Xv, color="C1", lw=1.4, label="output spectrum $|H\\,Y|$")
ax.plot(fv, np.abs(Hv), "k--", lw=0.9, label="$|H(f)|$, $f_c = 5$ Hz")
ax.set_xlim(0, 40); ax.set_ylim(0, 1.15); ax.set_xlabel("frequency (Hz)"); ax.set_ylabel("amplitude"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
The same equation runs backwards, and that is how sensors and structures are characterized. If you know the input as well as the output — a shaker table with an instrumented base, a calibrator driving a sensor — then \(H(f) = Y_{\mathrm{out}}(f)/Y_{\mathrm{in}}(f)\), frequency by frequency, is a measured transfer function: the sensor’s frequency response, or the structure’s, with its resonances standing up as peaks. Without an input channel one is left with the output spectrum alone, which mixes the structure’s \(|H|\) with whatever spectrum the excitation had; that is the situation in the case study that follows, and it is worth knowing what you are giving up.
A signal component passes through two first-order low-pass filters, each with its corner exactly at that component’s frequency. By how much is it attenuated, and what happens to its phase?
Each filter contributes \(-3\) dB and \(-45°\) at its own corner. By Equation 7.24 the magnitudes multiply — \((1/\sqrt{2})^2 = 1/2\), that is \(-6\) dB — and the phases add, to \(-90°\). In decibels the attenuations simply add, \(-3 + (-3) = -6\); what one must not do is add the linear magnitudes or multiply the phases.
A two-story frame sits on a shaker table with eight accelerometers on its floors, and the first question a structural engineer asks of such a record is “what are the natural frequencies?” — the numbers that decide whether a design resonates with traffic, wind, or an earthquake. We apply the windowed FFT to the eight channels, read off the dominant peaks, and use the sensor layout to attribute each peak to translation or rotation of a story.
The data were collected by Matteo Pozzi and are shared here with permission. The record holds 25 s of response at \(f_s = 200\) Hz — 5000 samples on each of eight channels, \(a_1\) through \(a_8\) — from a two-story frame excited through its base. Figure 7.12 shows where the sensors sit: on each story, two accelerometers measure in the \(Y\) direction at the two ends of one edge and two measure in the \(X\) direction at the two ends of a perpendicular edge, so that each floor’s motion in its own plane — two translations and a rotation — is captured by four sensors. Only the response was recorded; there is no channel for the shaker’s own motion, so this is output-only analysis, with the limitation just described in Section 7.5.
Look at the record before transforming it — always. Figure 7.13 shows channel \(a_7\): quiet for the first two seconds, a strong burst of motion through the middle of the record, and a decay. This is a transient, not a stationary signal, so there is no steady process whose average power we could estimate — the power spectral density of the optional callout in Section 7.4 would describe a process, not this event. What we want is the amplitude spectrum of Equation 7.17 of the whole record, Hann-windowed since the record certainly does not hold a whole number of cycles of anything, with the mean removed first. The full 25 s buys \(\delta f = 0.04\) Hz, resolution enough to separate peaks a few tenths of a hertz apart.
fig, ax = plt.subplots(figsize=(6.8, 2.4))
ax.plot(t_case, seis[:, 6], lw=0.5)
ax.set_xlabel("time (s)"); ax.set_ylabel("$a_7$")
plt.tight_layout(); plt.show()
The analysis is four lines, and they are worth reading rather than folding away: remove the mean, window, transform, scale.
N_case = len(seis)
w_case = np.hanning(N_case)
f_case = np.fft.rfftfreq(N_case, 1/fs_case)
def amp_spec(x):
x = x - x.mean() # remove the dc component
return amp_spectrum(x, w_case) # window, transform, scale to amplitude
spectra = np.column_stack([amp_spec(seis[:, c]) for c in range(8)])
dominant = f_case[np.argmax(spectra, axis=0)]
fig, axs = plt.subplots(4, 2, figsize=(7.4, 6.4), sharex=True, sharey=False)
for c, ax in enumerate(axs.T.flat):
ax.plot(f_case, spectra[:, c], lw=0.8)
ax.axvline(dominant[c], color="C3", lw=0.8, ls=":")
ax.set_title(f"$a_{c + 1}$: dominant peak {dominant[c]:.2f} Hz", fontsize=8)
ax.set_xlim(0, 8); ax.set_yticks([])
for ax in axs[-1]: ax.set_xlabel("frequency (Hz)")
plt.tight_layout(); plt.show()
Three things stand out in Figure 7.14. First, the dominant frequency depends on where you put the sensor: the \(Y\)-direction channels all peak at 2.48 Hz, the two \(X\) channels of story B at 3.68 Hz (\(a_7\)) and 4.24 Hz (\(a_8\)), and \(a_4\) at 4.24 Hz as well — a single accelerometer reports the mode it happens to sit on, which is why one never instruments a structure with one. Second, the same frequencies recur from channel to channel with different strengths — 2.48, 2.96, 3.20, 3.68, 4.24, and 5.04 Hz — which is the signature of modes: a natural frequency belongs to the whole structure, and each sensor sees it in proportion to how much that mode moves its location and direction. Third, the peaks are sharp and distinct because the record is long; the exercise at the end of the lecture shows what happens to them at a fifth and a twenty-fifth of the length.
The layout in Figure 7.12 lets us go one step further and separate the motions before transforming. Rotation of a rigid floor about its vertical axis moves the two ends of an edge in opposite senses, translation moves them together; so on story A the average \((a_1 + a_2)/2\) is the \(Y\)-acceleration of the floor with the rotation cancelled, and the difference \(a_2 - a_1\) is the rotation with the translation cancelled — and likewise \((a_3 + a_4)/2\) and \(a_4 - a_3\) for the \(X\) direction and the same rotation. Figure 7.15 shows the spectra of these combinations. The \(Y\)-translation of each story is dominated by 2.48 Hz alone; the \(X\)-translation of story B peaks at 3.68 Hz; and the rotation of story B is dominated by 5.04 and 4.24 Hz, peaks that were minor in every raw \(Y\) channel. (That \(a_2 - a_1\) and \(a_4 - a_3\) have identical spectral shapes is a check that the floor is behaving as a rigid body: both are the same rotation, scaled by the edge lengths.) Which mode a peak belongs to has become a question the geometry can answer.
combos = {"story A, $Y$ translation $(a_1 + a_2)/2$": (seis[:, 0] + seis[:, 1])/2,
"story A, $X$ translation $(a_3 + a_4)/2$": (seis[:, 2] + seis[:, 3])/2,
"story A, rotation $a_2 - a_1$": seis[:, 1] - seis[:, 0],
"story B, $Y$ translation $(a_5 + a_6)/2$": (seis[:, 4] + seis[:, 5])/2,
"story B, $X$ translation $(a_7 + a_8)/2$": (seis[:, 6] + seis[:, 7])/2,
"story B, rotation $a_6 - a_5$": seis[:, 5] - seis[:, 4]}
fig, axs = plt.subplots(3, 2, figsize=(7.4, 5.2), sharex=True)
for (name, x), ax in zip(combos.items(), axs.T.flat):
A = amp_spec(x); ax.plot(f_case, A, lw=0.8)
ax.axvline(f_case[np.argmax(A)], color="C3", lw=0.8, ls=":")
ax.set_title(f"{name}: peak {f_case[np.argmax(A)]:.2f} Hz", fontsize=8); ax.set_xlim(0, 8); ax.set_yticks([])
for ax in axs[-1]: ax.set_xlabel("frequency (Hz)")
plt.tight_layout(); plt.show()
What we have not done is prove that every peak is a mode. Output-only spectra are the product of the structure’s response and the excitation’s spectrum (Equation 7.23), and a shaker with a strong component at some frequency would leave a peak there too; a recorded base channel would settle it by division. Nor have we extracted damping or mode shapes, which are the business of a structural dynamics course. What the FFT has delivered, in four lines of code, is the list of frequencies at which this structure prefers to move and a first assignment of each to a story and a direction — which is what a sensing engineer hands to the structural engineer, and what the next lecture teaches you to attach an uncertainty to.
Cover Figure 7.14 and look only at the layout in Figure 7.12. Which channels should share a dominant peak (same direction, same story)? Which pairs should show identical spectral shapes if the floors are rigid? Which combination of channels would you form to see the rotation of story A on its own? Write your predictions down, then uncover the figures and check each one.
The lecture opened with a summary that could not tell a signal from its shuffled copy and closed with a spectrum that names a building’s modes, and every step between was one more assumption made explicit about a finite record. The Fourier series decomposes a periodic signal into harmonics; the transform is its limit for signals that never repeat; the DFT is the transform we can compute, at the price of pretending the record repeats; windows pay that price gracefully; the PSD, if you took the optional detour, extends all of it to random signals; and the convolution theorem turns every LTI block in the chain — sensor, filter, structure — into a multiplication that can be read off a Bode plot. Each of the three questions of Section 7.2 now has its place in the spectrum: the mean is bin zero, the RMS is the square root of the area, and how the signal moves is the shape of the curve in between.
This is the last purely signal-side tool in the SENSORS third. Lecture 7 asks how much of what the spectrum shows is the structure and how much is the measurement — noise, quantization, and the uncertainty that every number in Figure 7.14 deserves to carry. After that the DATA third takes custody of the record, and the spectrum will reappear when time-series data need to be stored, compressed, and modeled.
A pump’s output pressure is a full-wave-rectified 60 Hz sine, \(p(t) = |P_0 \sin 120\pi t|\).
A 60 Hz hum is sampled at 400 Hz for 64 samples.
Band variances are \(6 \times 10^{-4}\) and \(3 \times 10^{-4}\ (\mathrm{m/s^2})^2\), so the band RMS values are \(0.0245\) and \(0.0173\ \mathrm{m/s^2}\). Their sum, \(0.042\), exceeds the 1–10 Hz RMS of \(0.030\) because it is variances that add across disjoint bands (Equation 7.21), not RMS values: \(\sqrt{6 \times 10^{-4} + 3 \times 10^{-4}} = 0.030\).
A first-order low-pass filter with a 50 Hz corner precedes the ADC. A modal peak sits at 30 Hz.
Using channel \(a_7\) of the case-study data, compute the Hann-windowed amplitude spectrum of:
Report \(\delta f\) and the dominant peak of each, and explain what became of the peaks at 2.96, 3.20, and 3.68 Hz.
numpy.fft: https://numpy.org/doc/stable/reference/routines.fft.htmlscipy.signal.welch: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html