Machining Sounds - Part B: Sound¶
Part B: Sound¶
A FAB26 workshop lesson, companion to Machining Dynamics: Jupyter Notebook Edition (DOI: forthcoming)
Authors: Tony L. Schmitz and Michael F. Gomez
"The profound study of nature is the most fertile source of mathematical discoveries." Joseph Fourier, The Analytical Theory of Heat (1822)
Part A ended without fully describing the procedure to calculate the vibration spectrum of sound signals and the FRFs from tap tests. In Part B, we define the procedure. Starting from the sound that is incident on a microphone, we follow the analysis chain that turns a continuous sound into a recording with its own spectrum. Specifically, the continuous sound becomes a list of numbers in the recording, the periodic signal contained in the recording becomes a sum of sinusoids, and the discrete Fourier transform connects the two. We'll understand the compute_spectrum and compute_spectrogram functions that produced the figures from Part A. These short, inspectable computations will then become trusted analysis tools.
The purposes of Part B are to:
- Describe sound as a traveling pressure wave and identify how a microphone responds to sound.
- Explain sampling, the sample rate, the Nyquist frequency, and aliasing.
- Develop the Fourier series and the idea that a periodic signal can be represented by a sum of scaled harmonics.
- Define the discrete Fourier transform, its fast algorithm, and its frequency resolution.
- Build the spectrogram from short-time spectra and examine its time-frequency tradeoff.
- Apply the complete toolkit to a recording of your own, choosing the sampling and analysis parameters.
- Connect a periodic force to a harmonic sound, the observation that carries us into Part C.
Each concept is demonstrated interactively: run the cells, listen, drag a slider, and run again. An understanding of concepts from Part A is assumed; no other knowledge is required.
Nomenclature¶
The following symbols are used consistently throughout the lesson:
| Symbol | Meaning | Units |
|---|---|---|
| $t$ | time | s |
| $f$ | frequency | Hz |
| $c$ | speed of sound in air | m/s |
| $f_s$ | sample rate | Hz |
| $f_N$ | Nyquist frequency, $f_s / 2$ | Hz |
| $N$ | number of samples in a record | samples |
| $T$ | record length, $N / f_s$ | s |
| $\Delta f$ | frequency resolution, $1 / T$ | Hz |
| $f_0$ | fundamental frequency of a periodic signal | Hz |
| $k$ | harmonic index | dimensionless |
| $A_k$, $\phi_k$ | amplitude and phase of the $k$-th harmonic | signal units, rad |
| $x_n$ | the $n$-th sample of a recording | signal units |
| $X_m$ | the $m$-th coefficient of the DFT | signal units |
In the code cells, $f_s$ appears as the constant AUDIO_RATE, matching Part A; other symbols keep their names (f0, N, df).
Setup: Install and Import Libraries¶
Run the install cell once. If this is the first install on your machine, restart the kernel afterward (Kernel menu, Restart), then continue from the import cell.
%pip install -q numpy matplotlib plotly
Note: you may need to restart the kernel to use updated packages.
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from matplotlib.textpath import TextPath
from matplotlib.font_manager import FontProperties
from IPython.display import Audio, display
AUDIO_RATE = 48000 # Hz, sample rate for all synthesized sound
print('Libraries loaded.')
Libraries loaded.
The utility functions used throughout the lesson are defined in the collapsed cell below and documented in the Appendix. They are the same functions used in Part A.
# --- Utility functions used throughout the lesson (documented in the Appendix) ---
def fade(x, ms=15, fs=AUDIO_RATE):
"""Apply a short raised-cosine fade-in and fade-out so audio clips start
and stop without clicks."""
n = int(ms * 1e-3 * fs)
ramp = 0.5 * (1 - np.cos(np.pi * np.arange(n) / n))
y = x.copy()
y[:n] *= ramp
y[-n:] *= ramp[::-1]
return y
def compute_spectrum(x, fs, fmax=4000):
"""Amplitude spectrum of signal x on a linear scale, normalized to a peak
value of one. Returns the frequency vector (Hz, limited to fmax) and the
normalized spectrum."""
x = np.asarray(x, dtype=float)
w = np.hanning(len(x))
X = np.fft.rfft((x - x.mean()) * w)
f = np.fft.rfftfreq(len(x), 1 / fs)
amp = 2 * np.abs(X) / np.sum(w)
amp = amp / amp.max()
sel = f <= fmax
return f[sel], amp[sel]
def compute_spectrogram(x, fs, nfft=2048, hop=512, fmax=4000):
"""Short-time amplitude spectrogram on a linear scale, normalized to a
peak value of one. Returns (t_frames, f, S) where S has one column per
time frame and one row per frequency."""
w = np.hanning(nfft)
n_frames = 1 + max(0, (len(x) - nfft) // hop)
f = np.fft.rfftfreq(nfft, 1 / fs)
sel = f <= fmax
S = np.empty((int(sel.sum()), n_frames))
for j in range(n_frames):
S[:, j] = np.abs(np.fft.rfft(x[j * hop:j * hop + nfft] * w))[sel]
S /= S.max()
t_frames = (np.arange(n_frames) * hop + nfft / 2) / fs
return t_frames, f[sel], S
def decimate_for_plot(t, x, max_points=None):
"""Pass plotting data through unchanged by default (max_points=None).
To thin a long record for display, set max_points to an integer: each
output bin then keeps its local minimum and maximum, preserving the
visual envelope without aliasing. Computations always use the
full-rate data."""
n = len(x)
if max_points is None or n <= max_points:
return np.asarray(t), np.asarray(x)
bins = max_points // 2
edge = (n // bins) * bins
xb = np.asarray(x)[:edge].reshape(bins, -1)
tb = np.asarray(t)[:edge].reshape(bins, -1)
tt = np.repeat(tb.mean(axis=1), 2)
xx = np.empty(2 * bins)
xx[0::2] = xb.min(axis=1)
xx[1::2] = xb.max(axis=1)
return tt, xx
B.1 Sound as a Pressure Wave¶
A vibrating surface pushes on the air ahead of it as it moves forward (squeezing) and leaves a slight rarefaction as it moves back (expanding). Each layer of air nudges the next, and the disturbance travels outward at the speed of sound, about 343 m/s in room temperature air. As Sect. A.1 recounted, Newton computed this speed from the air's pressure and density in 1687 and Laplace corrected the missing thermal effect in 1816. The traveling disturbance is a pressure wave: at a fixed point in the room, the air pressure rises and falls slightly about its resting value as the wave passes.
A microphone measures exactly that. Its diaphragm is deformed by the pressure fluctuation and the recording is a single number, pressure, reported over and over as time advances. Two properties of the wave map onto what we hear. The rate of the fluctuation sets the pitch: human hearing spans roughly 20 Hz to 20 kHz and faster fluctuation is heard as higher pitch. The size of the fluctuation sets the loudness.
Look closely at what a sound recording contains, though: pressure against time and nothing else. In the time domain, there is no frequency axis, no peaks, no spectrum. Every frequency domain figure in Part A was computed from a record of exactly this kind and the remainder of Part B describes that computation from beginning to end. In the first step, we need to get the continuous time domain microphone signal into the phone or computer. However, a computer is digital. It does not work with continuous signals, so the continuous data must first be discretized into a list of values. We call this sampling the signal.
▸ Key terms: pressure wave, speed of sound, pitch, loudness
B.2 Sampling¶
A phone converts the continuous pressure signal into a series of numbers by measuring it at regular intervals. This sampling occurs $f_s$ times per second, the sample rate. Each measurement is a sample and the recording becomes the list $x_0, x_1, x_2, \ldots$, with the true wave known only at the sampling time and lost between samples. The measurement is accurate as long as the wave does not vary significantly between samples. We quantify the required sample rate using the sampling theorem of Nyquist and Shannon [1, 2], which states that frequencies contained in a continuous signal are only accurately represented up to the Nyquist frequency
$$f_N = \frac{f_s}{2} \tag{B.1}$$
and no higher. Phones sample sound at $f_s = 48000$ samples per second, so $f_N = 24$ kHz sits comfortably above the limit of human hearing, which is why a phone recording can accurately represent any audible sound. The explorer of Fig. B.1 samples a 440 Hz tone at an adjustable rate. The dashed curve is the lowest frequency wave consistent with the samples, which is the wave any playback or analysis will report. At sampling rates which exceed twice the tone's frequency the continuous signal is accurately represented by the samples. However, drag the slider below 880 samples per second and the sampled (digital) signal is no longer accurate.
# --- Fig. B.1: a tone sampled at an adjustable rate ---
f_tone = 440.0
t_fine = np.arange(int(0.01 * AUDIO_RATE)) / AUDIO_RATE # 10 ms at 48 kHz
FS_STEPS = [500, 600, 700, 800, 1000, 1500, 3000, 8000]
START_FS = 6 # slider starts at 3,000 per second, above the Nyquist limit
fig = go.Figure()
fig.add_trace(go.Scatter(x=t_fine * 1e3, y=np.sin(2 * np.pi * f_tone * t_fine),
mode='lines', line=dict(color='#1f77b4', width=1),
name='pressure wave, 440 Hz'))
for i, fs_i in enumerate(FS_STEPS):
t_s = np.arange(int(0.01 * fs_i) + 1) / fs_i
d = f_tone - fs_i * round(f_tone / fs_i) # signed frequency implied by the samples
vis = (i == START_FS)
fig.add_trace(go.Scatter(x=t_s * 1e3, y=np.sin(2 * np.pi * f_tone * t_s),
mode='markers', marker=dict(color='#d62728', size=8),
name='samples', visible=vis))
fig.add_trace(go.Scatter(x=t_fine * 1e3, y=np.sin(2 * np.pi * d * t_fine),
mode='lines', line=dict(color='#ff7f0e', dash='dash'),
name=f'lowest-frequency wave through the samples, {abs(d):.0f} Hz',
visible=vis))
steps = []
for i, fs_i in enumerate(FS_STEPS):
vis = [True] + [False] * (2 * len(FS_STEPS))
vis[1 + 2 * i:3 + 2 * i] = [True, True]
steps.append(dict(method='update', label=f'{fs_i:,}', args=[{'visible': vis}]))
fig.update_layout(
sliders=[dict(active=START_FS, currentvalue=dict(prefix='samples per second = '),
steps=steps)],
title='Fig. B.1 — A 440 Hz tone sampled at an adjustable rate',
template='simple_white', height=430,
xaxis_title='t (ms)', yaxis_title='Pressure (normalized)',
legend=dict(x=0.98, y=0.98, xanchor='right'),
margin=dict(l=60, r=20, t=60, b=90))
fig.show()
Below the Nyquist limit of 880 Hz for this example, the samples are recorded, but they no longer represent the continuous wave that produced them. At 600 samples per second, the explorer shows that the 440 Hz tone appears as a 160 Hz wave. This sampling failure is called aliasing and the digital recording reports the lower frequency, aliased signal in place of the true one. The worst part is that we do not know that we have a measurement error. Fig. B.2 contains an example from the audible frequency range. A 3200 Hz tone sampled at 4000 samples per second produces an 800 Hz tone.
# --- Fig. B.2: aliasing of a tone above the Nyquist frequency ---
f_true = 3200.0
fs_slow = 4000
t_fine = np.arange(int(0.005 * AUDIO_RATE)) / AUDIO_RATE # 5 ms
t_samp = np.arange(int(0.005 * fs_slow)) / fs_slow
f_alias = abs(f_true - fs_slow * round(f_true / fs_slow)) # 800 Hz
fig = go.Figure()
fig.add_trace(go.Scatter(x=t_fine * 1e3, y=np.sin(2 * np.pi * f_true * t_fine),
mode='lines', line=dict(color='#1f77b4', width=1),
name=f'true tone, {f_true:.0f} Hz'))
fig.add_trace(go.Scatter(x=t_fine * 1e3, y=-np.sin(2 * np.pi * f_alias * t_fine),
mode='lines', line=dict(color='#ff7f0e', dash='dash'),
name=f'aliased signal, {f_alias:.0f} Hz'))
fig.add_trace(go.Scatter(x=t_samp * 1e3, y=np.sin(2 * np.pi * f_true * t_samp),
mode='markers', marker=dict(color='#d62728', size=8),
name='samples at 4,000 per second'))
fig.update_layout(title='Fig. B.2 — Aliasing of a 3,200 Hz tone sampled at 4,000 per second',
template='simple_white', height=380,
xaxis_title='t (ms)', yaxis_title='Pressure (normalized)',
legend=dict(x=0.98, y=0.02, xanchor='right', yanchor='bottom'),
margin=dict(l=60, r=20, t=60, b=50))
fig.show()
Both curves pass through every sample, but the sampled record contains no information to indicate that something has gone amiss. The first clip in the following is the true 3200 Hz tone. The second is what a system sampling at 4000 per second records. Although the starting signal is the same, the dominant tone in the digital signal is the aliased version at 800 Hz.
# --- Hearing aliasing ---
dur = 1.2
t = np.arange(int(dur * AUDIO_RATE)) / AUDIO_RATE
true_tone = 0.5 * np.sin(2 * np.pi * 3200 * t)
fs_slow = 4000
t_slow = np.arange(int(dur * fs_slow)) / fs_slow
sampled = 0.5 * np.sin(2 * np.pi * 3200 * t_slow)
crude = np.repeat(sampled, AUDIO_RATE // fs_slow) # each sample held until the next
print('The true tone, 3,200 Hz:')
display(Audio(fade(true_tone), rate=AUDIO_RATE))
print('The aliased recording (dominant aliased tone at 800 Hz, plus playback artifacts):')
display(Audio(fade(crude), rate=AUDIO_RATE))
The true tone, 3,200 Hz:
The aliased recording (dominant aliased tone at 800 Hz, plus playback artifacts):
Audio equipment settles this problem by its choice of sample rate. Phones use 48000 samples per second, and audio compact discs, the standard before streaming, use 44100 samples per second. These sample rates are selected so the Nyquist frequency is above the roughly 20 kHz limit of human hearing. Furthermore, the recording hardware filters out content above $f_N$ before sampling so aliasing is eliminated. Although Eq. B.1 requires only that $f_s$ exceed twice the highest frequency of interest, measurement practice is more conservative. For vibration measurements we typically sample at roughly 10 times the highest frequency we expect to observe, so the recorded waveform is smooth and accurate. We also filter the continuous signal at the Nyquist frequency before sampling. For machine listening, the sound sampling rates happen to be adequate. Tool tip natural frequencies, spindle rotation rates, and tooth passing frequencies rarely exceed a few kHz and a phone samples at 10 times that comfortably. Every recording in this lesson series uses $f_s = 48000$.
▸ Key terms: sample, sample rate, Nyquist frequency, aliasing
B.3 The Fourier Series¶
After sampling, the recording is now a list of numbers and the question from Sect. B.1 stands: where does frequency come from? Let's begin with a simple observation about the pitch of sounds. A sustained tone from a string, a pipe, or a humming machine produces a pressure signal that repeats a pattern, of whatever shape, every $T_0$ seconds. The repetition rate $f_0 = 1/T_0$ is the fundamental frequency and it is what our ears interpret as the pitch.
Fourier's assertion of 1822, with the precise conditions later supplied by Dirichlet, is that any such periodic signal can be written as a sum of sinusoids at integer multiples of the fundamental frequency [3]:
$$x(t) = A_0 + \sum_{k=1}^{\infty} A_k \sin\!\left(2\pi k f_0 t + \phi_k\right) \tag{B.2}$$
The component at $k f_0$ is the $k$-th harmonic, and the list of amplitudes $A_1, A_2, A_3, \ldots$ describes how strongly each harmonic is present. These concepts were referenced in Part A, where we informally noticed which harmonics were present and their magnitudes. Equation B.2 is a synthesis statement. It says a periodic signal can be represented as a sum of scaled harmonics. This assembly can be thought of as a recipe, and the assembly can be observed or heard one term at a time. The following explorer synthesizes a sawtooth-like wave whose recipe is $A_k = 1/k$; drag the slider to add terms (harmonics).
# --- Fig. B.3: a periodic signal assembled from its harmonics ---
f0 = 250.0
K_STEPS = [1, 2, 3, 4, 5, 8, 12, 20, 40]
START_K = 2 # slider starts at K = 3
t_two = np.arange(int(2 / f0 * AUDIO_RATE)) / AUDIO_RATE # two periods
fig = make_subplots(rows=1, cols=2, column_widths=[0.55, 0.45],
subplot_titles=('Waveform (two periods)', 'Recipe used so far'))
for i, K in enumerate(K_STEPS):
xK = np.zeros_like(t_two)
for k in range(1, K + 1):
xK += (1 / k) * np.sin(2 * np.pi * k * f0 * t_two)
vis = (i == START_K)
fig.add_trace(go.Scatter(x=t_two * 1e3, y=xK, mode='lines',
line=dict(color='#1f77b4'), showlegend=False,
visible=vis), row=1, col=1)
fig.add_trace(go.Bar(x=[k * f0 for k in range(1, K + 1)],
y=[1 / k for k in range(1, K + 1)],
width=30, marker_color='#1f77b4', showlegend=False,
visible=vis), row=1, col=2)
steps = []
for i, K in enumerate(K_STEPS):
vis = [False] * (2 * len(K_STEPS))
vis[2 * i:2 * i + 2] = [True, True]
steps.append(dict(method='update', label=str(K), args=[{'visible': vis}]))
fig.update_layout(
sliders=[dict(active=START_K, currentvalue=dict(prefix='harmonics K = '), steps=steps)],
title='Fig. B.3 — A periodic signal assembled from its harmonics',
template='simple_white', height=420,
margin=dict(l=60, r=20, t=70, b=90))
fig.update_xaxes(title_text='t (ms)', row=1, col=1)
fig.update_xaxes(title_text='Frequency (Hz)', range=[0, 41 * f0], row=1, col=2)
fig.update_yaxes(title_text='x', row=1, col=1)
fig.update_yaxes(title_text='Amplitude, A<sub>k</sub>', range=[0, 1.05], row=1, col=2)
fig.show()
The convergence to the actual sawtooth-like wave is audible as more harmonics are added. Each clip that follows represents a sum of harmonics: one harmonic (the fundamental frequency only) gives a pure tone, a few more harmonics add buzz, and by 40 harmonics the timbre has stopped changing to most ears, since the remaining harmonics are weak and high frequency. The pitch is identical in every clip, because the fundamental frequency does not change.
# --- Hearing the partial sums converge ---
dur = 1.2
t = np.arange(int(dur * AUDIO_RATE)) / AUDIO_RATE
for K in [1, 3, 10, 40]:
xK = np.zeros_like(t)
for k in range(1, K + 1):
xK += (1 / k) * np.sin(2 * np.pi * k * 250 * t)
xK *= 0.5 / np.max(np.abs(xK))
print(f'K = {K} harmonic{"s" if K > 1 else ""}:')
display(Audio(fade(xK), rate=AUDIO_RATE))
K = 1 harmonic:
K = 3 harmonics:
K = 10 harmonics:
K = 40 harmonics:
Timbre¶
Equation B.2 also explains why instruments sound different when playing the same note. Ohm proposed in 1843 that the ear decomposes a sound into its harmonics and Helmholtz established the case in 1863. The pitch is set by the fundamental frequency and the timbre, the recognizable character of the sound, is set by the recipe of harmonic amplitudes [4]. Two signals with the same $f_0$ and different recipes have the same pitch and different timbres. Fig. B.4 compares two: the $1/k$ recipe from Fig. B.3 with every harmonic present, and a recipe with the same amplitudes but with only the odd harmonics ($k = 1, 3, 5, \ldots$).
# --- Fig. B.4: two recipes at the same pitch ---
f0 = 250.0
dur = 1.2
t = np.arange(int(dur * AUDIO_RATE)) / AUDIO_RATE
t_two = np.arange(int(2 / f0 * AUDIO_RATE)) / AUDIO_RATE
recipes = [('every harmonic, 1/k', range(1, 9), '#1f77b4', 'solid'),
('odd harmonics, 1/k', range(1, 9, 2), '#ff7f0e', 'dash')]
fig = make_subplots(rows=1, cols=2, column_widths=[0.55, 0.45],
subplot_titles=('Waveform (two periods)', 'Spectrum'))
sounds = {}
for name, ks, color, dash in recipes:
x_sig = np.zeros_like(t)
x_two = np.zeros_like(t_two)
for k in ks:
x_sig += (1 / k) * np.sin(2 * np.pi * k * f0 * t)
x_two += (1 / k) * np.sin(2 * np.pi * k * f0 * t_two)
sounds[name] = 0.5 * x_sig / np.max(np.abs(x_sig))
fig.add_trace(go.Scatter(x=t_two * 1e3, y=x_two, mode='lines',
line=dict(color=color, dash=dash), name=name),
row=1, col=1)
f_sp, a_sp = compute_spectrum(x_sig, AUDIO_RATE, fmax=2500)
fig.add_trace(go.Scatter(x=f_sp, y=a_sp, mode='lines',
line=dict(color=color, dash=dash), showlegend=False),
row=1, col=2)
fig.update_layout(title='Fig. B.4 — Two recipes at the same pitch',
template='simple_white', height=400,
legend=dict(x=0.44, y=0.98, xanchor='right'),
margin=dict(l=60, r=20, t=70, b=50))
fig.update_xaxes(title_text='t (ms)', row=1, col=1)
fig.update_xaxes(title_text='Frequency (Hz)', row=1, col=2)
fig.update_yaxes(title_text='x', row=1, col=1)
fig.update_yaxes(title_text='Amplitude (normalized)', range=[0, 1.05], row=1, col=2)
fig.show()
for name in sounds:
print(name + ':')
display(Audio(fade(sounds[name]), rate=AUDIO_RATE))
every harmonic, 1/k:
odd harmonics, 1/k:
Both clips have a pitch of 250 Hz, but they are easy to tell apart. The odd-harmonic recipe has the hollow character of a clarinet-like tone, while the full recipe has the brighter buzz. This observation matters for machine listening as well. A milling machine has a voice for the same reason an instrument does: a repeating process with a particular recipe of harmonics, and changes in that recipe are audible even when the pitch is unchanged.
Drawing with Harmonics¶
Equation B.2 applies to anything periodic and a closed curve traced in the x-y plane is periodic in its own way. The curve travels around the shape outline and the path repeats. Writing the path as a complex number $z(u) = x(u) + i\,y(u)$ over one lap, the same synthesis sum reconstructs the outline from a recipe of rotating components. Truncating the recipe at $K$ harmonics gives a smoothed approximation that sharpens as $K$ grows. The following explorer applies this to the outlines of the letters and numbers that compose FAB26. Watch where the reconstruction struggles the longest. We see that the sharp corners demand a larger number of harmonics. This is similar to sounds, such as an impulse or a chirp, which are only accurately represented by many harmonics.
# --- Fig. B.5: the letters FAB26 traced by partial Fourier sums ---
tp = TextPath((0, 0), 'FAB26', size=100,
prop=FontProperties(family='DejaVu Sans', weight='bold'))
contours = [p for p in tp.to_polygons() if len(p) >= 4]
def resample_contour(poly, M=400):
d = np.linalg.norm(np.diff(poly, axis=0), axis=1)
s = np.concatenate([[0], np.cumsum(d)])
u = np.linspace(0, s[-1], M, endpoint=False)
return np.interp(u, s, poly[:, 0]) + 1j * np.interp(u, s, poly[:, 1])
def partial_sum(z, K):
M = len(z)
C = np.fft.fft(z) / M
ks = np.fft.fftfreq(M, 1 / M).astype(int)
n = np.arange(M)
zr = np.zeros(M, dtype=complex)
for i in np.where(np.abs(ks) <= K)[0]:
zr += C[i] * np.exp(2j * np.pi * ks[i] * n / M)
return zr
Z = [resample_contour(p) for p in contours]
K_IMG = [1, 2, 3, 5, 8, 12, 20, 40, 80]
START_IMG = 2 # slider starts at K = 3
fig = go.Figure()
# the target outline, always visible
gx, gy = [], []
for z in Z:
gx += list(np.append(z.real, z.real[0])) + [None]
gy += list(np.append(z.imag, z.imag[0])) + [None]
fig.add_trace(go.Scatter(x=gx, y=gy, mode='lines',
line=dict(color='#7f7f7f', width=1), opacity=0.35,
hoverinfo='skip', showlegend=False))
for i, K in enumerate(K_IMG):
rx, ry = [], []
for z in Z:
zr = partial_sum(z, K)
rx += list(np.append(zr.real, zr.real[0])) + [None]
ry += list(np.append(zr.imag, zr.imag[0])) + [None]
fig.add_trace(go.Scatter(x=rx, y=ry, mode='lines',
line=dict(color='#1f77b4', width=2),
visible=(i == START_IMG), showlegend=False))
steps = []
for i, K in enumerate(K_IMG):
vis = [True] + [False] * len(K_IMG)
vis[1 + i] = True
steps.append(dict(method='update', label=str(K), args=[{'visible': vis}]))
fig.update_layout(
sliders=[dict(active=START_IMG, currentvalue=dict(prefix='harmonics per contour K = '),
steps=steps)],
title='Fig. B.5 — The letters FAB26 traced by partial Fourier sums',
template='simple_white', height=420,
xaxis=dict(visible=False),
yaxis=dict(visible=False, scaleanchor='x', scaleratio=1),
margin=dict(l=20, r=20, t=60, b=90))
fig.show()
When only a few harmonics are added per contour, every letter is a rounded blob. The straight lines arrive first and the corners arrive last. The audio version of a corner is a click or an abrupt change. Any signal feature that changes quickly, in space or in time, requires a high number of harmonics to represent it accurately. This is one reason the hard hammer tip of Sect. A.7 excited a wider frequency range than the soft tip. The short, abrupt impulse applied by the hard tip is the temporal cousin of a spatial sharp corner.
▸ Key terms: periodic signal, fundamental frequency, harmonic, Fourier series, timbre
B.4 The DFT and the FFT¶
Equation B.2 synthesizes a signal using its recipe. Measurement and signal analysis proceeds in the opposite direction. In this case, we first record the signal and then we determine its recipe. For a sampled record of $N$ numbers, this inverse operation is the discrete Fourier transform (DFT):
$$X_m = \sum_{n=0}^{N-1} x_n \, e^{-i 2\pi m n / N} \tag{B.3}$$
The DFT takes the $N$ samples $x_n$ and returns one coefficient per frequency bin. Bin $m$ corresponds to the frequency
$$f_m = \frac{m f_s}{N} \tag{B.4}$$
and $|X_m|$ reports the size of that frequency in the record. The transform retains all the information in the recording, because the samples are recovered exactly by the inverse DFT:
$$x_n = \frac{1}{N} \sum_{m=0}^{N-1} X_m \, e^{+i 2\pi m n / N} \tag{B.5}$$
which is the discrete counterpart of the synthesis sum of Eq. B.2. Analysis and synthesis are the two directions of one relationship. From Eq. B.4, adjacent bins are separated by the frequency resolution, set by the record length $T = N / f_s$:
$$\Delta f = \frac{f_s}{N} = \frac{1}{T} \tag{B.6}$$
A one-second record gives a frequency spectrum with points separated by 1 Hz. In other words, the frequency resolution is 1 Hz. If the record length is only 0.1 s, the resolution is reduced to 10 Hz. Longer time records provide higher resolution and more frequency detail, but there is a tradeoff we will study in Sect. B.5. Finally, for a sinusoidal component landing in bin $m$, the amplitude follows from the coefficient as
$$A_m = \frac{2\,|X_m|}{\sum_n w_n} \tag{B.7}$$
where $w_n$ is the window discussed next. Without a window ($w_n = 1$) the denominator is simply $N$. Computed directly, Eq. B.3 costs about $N^2$ multiplications, which for a one-second phone recording is over two billion. The fast Fourier transform (FFT) of Cooley and Tukey computes the identical result in about $N \log_2 N$ operations [5], roughly a thousandfold savings at that size, and it is the reason the spectrum is a routine, real-time computation. In the code, np.fft.rfft is the FFT, specialized to real-valued inputs.
The compute_spectrum Function¶
The spectra of Part A, the middle column of Fig. A.5 and the force spectra of Fig. A.2, were produced by compute_spectrum:
def compute_spectrum(x, fs, fmax=4000):
x = np.asarray(x, dtype=float)
w = np.hanning(len(x)) # window: taper the record ends
X = np.fft.rfft((x - x.mean()) * w) # Eq. B.3, by the FFT
f = np.fft.rfftfreq(len(x), 1 / fs) # the bin frequencies of Eq. B.4
amp = 2 * np.abs(X) / np.sum(w) # amplitude by Eq. B.7
amp = amp / amp.max() # normalize the peak to one
sel = f <= fmax
return f[sel], amp[sel]
One line deserves an explanation. The window w tapers the record smoothly to zero at both ends before the transform is applied. The DFT treats the record as one period of a repeating signal and a record whose ends do not match creates an artificial jump at their interface. The tapering at the beginning and end removes the jump at the cost of slightly widening each spectral peak. Every spectrum in this lesson series uses this window.
With the algorithm available, we can test it against a known answer. The recipe of Fig. B.3 was $A_k = 1/k$; the following cell synthesizes that signal, computes its spectrum, and reads the amplitudes back at each harmonic.
# --- Reading a known recipe back from its spectrum ---
f0 = 250.0
dur = 1.0
t = np.arange(int(dur * AUDIO_RATE)) / AUDIO_RATE
x_sig = np.zeros_like(t)
for k in range(1, 9):
x_sig += (1 / k) * np.sin(2 * np.pi * k * f0 * t)
f_sp, a_sp = compute_spectrum(x_sig, AUDIO_RATE, fmax=2500)
df = AUDIO_RATE / len(t)
print(f'Record length T = {dur:.1f} s, so the resolution of Eq. B.6 is {df:.1f} Hz per bin.')
print()
print('harmonic frequency recipe A_k measured amplitude')
for k in range(1, 9):
idx = np.argmin(np.abs(f_sp - k * f0))
print(f' {k} {k * f0:6.0f} Hz {1 / k:.3f} {a_sp[idx]:.3f}')
Record length T = 1.0 s, so the resolution of Eq. B.6 is 1.0 Hz per bin. harmonic frequency recipe A_k measured amplitude 1 250 Hz 1.000 1.000 2 500 Hz 0.500 0.500 3 750 Hz 0.333 0.333 4 1000 Hz 0.250 0.250 5 1250 Hz 0.200 0.200 6 1500 Hz 0.167 0.167 7 1750 Hz 0.143 0.143 8 2000 Hz 0.125 0.125
The measured amplitudes reproduce the recipe, normalized so the fundamental value is one. The spectrum is analysis in exactly the sense that Eq. B.2 is synthesis and the frequency axes of Part A are the output of this computation and nothing else.
The speed of the FFT is worth seeing once. The following cell transforms a record with a quarter million samples and reports the time.
# --- The cost of a spectrum ---
import time
N = 2**18 # 262,144 samples, about 5.5 s of audio
x_big = np.sin(2 * np.pi * 440 * np.arange(N) / AUDIO_RATE)
t0 = time.perf_counter()
np.fft.rfft(x_big)
t1 = time.perf_counter()
print(f'FFT of N = {N:,} samples: {(t1 - t0) * 1e3:.1f} ms')
print(f'The direct sum of Eq. B.3 needs about N^2 = {N**2:.1e} multiplications;')
print(f'the FFT needs about N log2 N = {N * np.log2(N):.1e}, a factor of '
f'{N / np.log2(N):,.0f} fewer.')
FFT of N = 262,144 samples: 4.8 ms The direct sum of Eq. B.3 needs about N^2 = 6.9e+10 multiplications; the FFT needs about N log2 N = 4.7e+06, a factor of 14,564 fewer.
▸ Key terms: discrete Fourier transform (DFT), fast Fourier transform (FFT), frequency resolution, window
B.5 The Spectrogram¶
The DFT answers one question well: which frequencies does the record contain overall? Sect. A.8 identified its blind spot. The spectrum of the musical scale reported every note and nothing about their order. The spectrogram addresses this blind spot by dividing the record into short frames, computing a spectrum for each, and stacking the results side by side. Each frame is advanced from the last by a fixed hop. The compute_spectrogram utility performs these steps. It loops through frames with one windowed FFT per frame. It then assembles the columns of frequencies for each frame into the heatmap image.
Equation B.6 defines the spectrogram limitation. There is a price to pay for the short frames. Because each frame is short, each frame's resolution is low: a 2048-sample frame at 48000 samples per second spans 43 ms and offers a resolution of only about 23 Hz. Shrinking the frame time provides more insight into signals that vary over time, but it also blurs the frequencies; expanding the frame does the opposite. This time-frequency tradeoff is a property of the analysis and it depends on the frame length. The following explorer shows the same musical scale from Sect. A.8 analyzed using three frame lengths; drag the slider between them.
# --- Fig. B.6: one scale, three frame lengths ---
up = [261.63, 293.66, 329.63, 349.23, 392.00, 440.00, 493.88, 523.25]
notes = up + up[-2::-1]
dur_scale = 3.0
note_len = int(dur_scale / len(notes) * AUDIO_RATE)
scale_sig = np.zeros(len(notes) * note_len)
for j, f_note in enumerate(notes):
seg = np.arange(note_len) / AUDIO_RATE
scale_sig[j * note_len:(j + 1) * note_len] = fade(
0.6 * np.sin(2 * np.pi * f_note * seg), ms=10)
NFFTS = [512, 2048, 8192]
START_N = 1 # slider starts at 2,048
fig = go.Figure()
for i, nfft in enumerate(NFFTS):
t_g, f_g, S = compute_spectrogram(scale_sig, AUDIO_RATE,
nfft=nfft, hop=nfft // 4, fmax=700)
fig.add_trace(go.Heatmap(x=t_g, y=f_g, z=S, colorscale='Inferno',
showscale=False, visible=(i == START_N)))
steps = []
for i, nfft in enumerate(NFFTS):
vis = [False] * len(NFFTS)
vis[i] = True
frame_ms = nfft / AUDIO_RATE * 1e3
df_bin = AUDIO_RATE / nfft
steps.append(dict(method='update',
label=f'{nfft} ({frame_ms:.0f} ms, {df_bin:.1f} Hz)',
args=[{'visible': vis}]))
fig.update_layout(
sliders=[dict(active=START_N, currentvalue=dict(prefix='frame length = '), steps=steps)],
title='Fig. B.6 — One scale at three frame lengths',
template='simple_white', height=430,
xaxis_title='t (s)', yaxis_title='Frequency (Hz)',
margin=dict(l=60, r=20, t=60, b=90))
fig.show()
At 512 samples per frame the note changes are crisply timed, but each note is represented by a thick, uncertain frequency band. At 8192 samples, on the other hand, the notes are represented by thin, precise frequency lines, but the beginning and end of each note smears into its neighbors. The middle setting of 2048 samples is the one used throughout Part A. It was chosen because tones that are 10s of Hz apart and events 10s of ms apart both matter for machine sounds. The softened staircase edges visible in Fig. A.5 have the same origin: each frame that straddles a note change reports both notes in the local spectrum.
▸ Key terms: frame, hop, time-frequency tradeoff
B.6 Analyzing Your Own Recording¶
Everything in this Part now runs in one direction. We begin with a continuous signal, sample that signal, calculate the DFT, and present the three views of the information. In this section you will apply the toolkit to your own recording, with the analysis decisions in your hands. Record a few seconds of sound on your phone, such as a whistle, a sung vowel, a clap, or a running machine, export or convert it to a WAV file, and place it beside this notebook. If no file is found, the cell synthesizes an example recording containing a hummed tone, a clap, and a rising whistle, writes it as example_recording.wav, and analyzes that instead. The section therefore runs either way and switching to your own sound is a one-string edit.
Four parameters at the top of the cell are yours to adjust and each one represents a concept that we've introduced. The analysis window (t_start, duration) selects which portion of the recording to analyze, and its length $T$ sets the frequency resolution based on Eq. B.6. The resample rate (fs_new) reduces the sample rate by keeping every $n$-th sample, deliberately without the protective filter that real recording hardware applies, so content above the new Nyquist frequency of Eq. B.1 aliases as expected and you can produce the failure of Sect. B.2 on your own sound. The frame length (nfft) adjusts the spectrogram on the time-frequency tradeoff we discussed in Sect. B.5.
# --- Reading and synthesizing recordings ---
import wave as _wave
def load_wav(path):
"""Read a WAV file (PCM 8, 16, 24, or 32 bit, mono or stereo) and return
(sample_rate, mono_signal) with the signal scaled to the range -1 to 1."""
with _wave.open(path, 'rb') as wf:
fs = wf.getframerate()
n_ch = wf.getnchannels()
sw = wf.getsampwidth()
raw = wf.readframes(wf.getnframes())
if sw == 1:
x = (np.frombuffer(raw, dtype=np.uint8).astype(float) - 128.0) / 128.0
elif sw == 2:
x = np.frombuffer(raw, dtype='<i2').astype(float) / 32768.0
elif sw == 3:
b = np.frombuffer(raw, dtype=np.uint8).reshape(-1, 3).astype(np.int32)
x = b[:, 0] | (b[:, 1] << 8) | (b[:, 2] << 16)
x = np.where(x >= 2**23, x - 2**24, x).astype(float) / 2**23
elif sw == 4:
x = np.frombuffer(raw, dtype='<i4').astype(float) / 2**31
else:
raise ValueError(f'unsupported WAV sample width: {sw} bytes')
if n_ch > 1:
x = x.reshape(-1, n_ch).mean(axis=1)
return fs, x
def make_example_recording(path='example_recording.wav'):
"""Synthesize a 4 s example (hummed tone, clap, rising whistle), write it
as a 16-bit WAV beside the notebook, and return (sample_rate, signal)."""
fs = AUDIO_RATE
t = np.arange(int(4.0 * fs)) / fs
x = np.zeros_like(t)
seg = (t >= 0.2) & (t < 1.4) # hummed tone, 220 Hz recipe
for k in range(1, 7):
x[seg] += (0.5 / k) * np.sin(2 * np.pi * k * 220 * t[seg])
seg = (t >= 1.7) & (t < 1.9) # clap
rng = np.random.default_rng(26)
x[seg] += np.exp(-(t[seg] - 1.7) / 0.02) * rng.standard_normal(int(seg.sum())) * 0.8
seg = (t >= 2.2) & (t < 3.8) # rising whistle, 800 to 1600 Hz
ts = t[seg] - 2.2
f_w = 800 + 500 * ts
x[seg] += 0.5 * np.sin(2 * np.pi * np.cumsum(f_w) / fs)
x = fade(0.9 * x / np.max(np.abs(x)))
with _wave.open(path, 'wb') as wf:
wf.setnchannels(1)
wf.setsampwidth(2)
wf.setframerate(fs)
wf.writeframes((x * 32767).astype('<i2').tobytes())
return fs, x
# --- Analyze a recording. Adjust the parameters and rerun. ---
filename = 'my_recording.wav' # your WAV file; the example is used when absent
t_start = 0.0 # s, where the analysis window begins
duration = 4.0 # s, window length; sets the resolution of Eq. B.6
fs_new = None # Hz, resample rate, e.g. 4000; None keeps the original
nfft = 2048 # samples per spectrogram frame, Sect. B.5
fmax = 4000 # Hz, top of the frequency axes
import os
if os.path.exists(filename):
fs0, x_rec = load_wav(filename)
print(f'Loaded {filename}.')
else:
fs0, x_rec = make_example_recording()
print(f'{filename} not found; using the synthesized example '
'(written beside this notebook as example_recording.wav).')
i0 = int(t_start * fs0)
x_rec = x_rec[i0:i0 + int(duration * fs0)]
x_rec = x_rec / (np.max(np.abs(x_rec)) + 1e-12)
fs_rec = float(fs0)
if fs_new is not None and fs_new < fs0:
stride = int(round(fs0 / fs_new))
x_rec = x_rec[::stride] # deliberate, unguarded resampling
fs_rec = fs0 / stride
elif fs_new is not None:
print('fs_new is at or above the original rate; keeping the original.')
T_rec = len(x_rec) / fs_rec
f_top = min(fmax, fs_rec / 2)
print(f'Analysis window T = {T_rec:.2f} s at fs = {fs_rec:,.0f} per second: '
f'Nyquist {fs_rec / 2:,.0f} Hz (Eq. B.1), resolution {1 / T_rec:.2f} Hz (Eq. B.6).')
display(Audio(fade(0.8 * x_rec), rate=int(fs_rec)))
tt, xx = decimate_for_plot(np.arange(len(x_rec)) / fs_rec, x_rec,
max_points=6000) # display only; see Appendix
f_sp, a_sp = compute_spectrum(x_rec, fs_rec, fmax=f_top)
t_g, f_g, S = compute_spectrogram(x_rec, fs_rec, nfft=nfft, hop=nfft // 4, fmax=f_top)
fig = make_subplots(rows=1, cols=3,
column_titles=('time domain', 'spectrum', 'spectrogram'))
fig.add_trace(go.Scatter(x=tt, y=xx, mode='lines',
line=dict(color='#1f77b4', width=0.7),
showlegend=False), row=1, col=1)
fig.add_trace(go.Scatter(x=f_sp, y=a_sp, mode='lines',
line=dict(color='#1f77b4'), showlegend=False), row=1, col=2)
fig.add_trace(go.Heatmap(x=t_g, y=f_g, z=S, colorscale='Inferno',
showscale=False), row=1, col=3)
fig.update_xaxes(title_text='t (s)', row=1, col=1)
fig.update_xaxes(title_text='Frequency (Hz)', row=1, col=2)
fig.update_xaxes(title_text='t (s)', row=1, col=3)
fig.update_yaxes(title_text='Amplitude', row=1, col=1)
fig.update_yaxes(title_text='Amplitude (normalized)', range=[0, 1.05], row=1, col=2)
fig.update_yaxes(title_text='Frequency (Hz)', range=[0, f_top], row=1, col=3)
fig.update_layout(title='Fig. B.7 — A recording of your choice in three views',
template='simple_white', height=380,
margin=dict(l=60, r=20, t=70, b=50))
fig.show()
my_recording.wav not found; using the synthesized example (written beside this notebook as example_recording.wav). Analysis window T = 4.00 s at fs = 48,000 per second: Nyquist 24,000 Hz (Eq. B.1), resolution 0.25 Hz (Eq. B.6).
Use the example, or your own recording, to review the concepts described in this chapter. The hummed tone shows multiple harmonics in the spectrum, which represents the recipe of Sect. B.3. The clap is observed as a brief vertical stripe on the spectrogram with broad content in the spectrum. This is the temporal cousin of the spatial sharp corners from the Corners and Clicks box. The whistle with an increasing frequency draws a rising line that the whole-record spectrum smears into a block, which was the lesson of Sect. A.8.
Next, vary the parameters. Shorten duration to 0.5 s and watch every spectral peak widen as the resolution is reduced according to Eq. B.6. Set fs_new = 4000 and listen: any content above the new 2000 Hz Nyquist frequency appears at a lower aliased frequency, and the whistle's line bends where it crosses the Nyquist frequency boundary. Adjust nfft between 512 and 8192 and update your spectrogram based on the time-frequency tradeoff. Each adjustment is a decision a measurement engineer makes before trusting a spectrum and Eqs. B.1 and B.6 predict the consequence of these decisions.
Phone voice memo apps commonly export M4A, rather than WAV, format files. Most phones and computers convert with a share or export option and the free tool ffmpeg converts anything, including the soundtrack of a video, with one command:
ffmpeg -i input.mp4 -ac 1 output.wav. The same analysis chain, microphone to spectrum and spectrogram, runs continuously for live listening. Part C examines it for a machining process.
B.7 The Sound of a Periodic Force¶
One observation now transitions us from Part B into Part C. Milling is interrupted cutting: each tooth on the rotating cutter enters the workpiece, cuts material away in the form of a chip, and exits the workpiece. The next tooth repeats the cycle. The cutting force is therefore a train of similar pulses arriving at a fixed rate. A pulse train is a periodic signal whose fundamental frequency equals that rate. By Eq. B.2, its content is composed of harmonics at integer multiples of the fundamental frequency, so the sound of a stable milling cut is represented by evenly spaced spectral lines. Fig. B.8 shows the principle using a pulse train with a fundamental frequency of 250 Hz.
# --- Fig. B.8: a repeating pulse and its spectrum ---
f_rep = 250.0
dur = 1.0
t = np.arange(int(dur * AUDIO_RATE)) / AUDIO_RATE
width = 0.0004 # each pulse lasts 0.4 ms
phase = (t * f_rep) % 1.0
pulses = np.where(phase < width * f_rep, 1.0, 0.0)
fig = make_subplots(rows=1, cols=2, column_widths=[0.45, 0.55],
subplot_titles=('Pulse train (first 20 ms)', 'Spectrum'))
n_show = int(0.02 * AUDIO_RATE)
fig.add_trace(go.Scatter(x=t[:n_show] * 1e3, y=pulses[:n_show], mode='lines',
line=dict(color='#1f77b4'), showlegend=False), row=1, col=1)
f_sp, a_sp = compute_spectrum(pulses, AUDIO_RATE, fmax=4000)
fig.add_trace(go.Scatter(x=f_sp, y=a_sp, mode='lines',
line=dict(color='#1f77b4'), showlegend=False), row=1, col=2)
fig.update_layout(title='Fig. B.8 — A repeating pulse and its spectrum',
template='simple_white', height=380,
margin=dict(l=60, r=20, t=70, b=50))
fig.update_xaxes(title_text='t (ms)', row=1, col=1)
fig.update_xaxes(title_text='Frequency (Hz)', row=1, col=2)
fig.update_yaxes(title_text='Force (normalized)', row=1, col=1)
fig.update_yaxes(title_text='Amplitude (normalized)', range=[0, 1.05], row=1, col=2)
fig.show()
print('The pulse train, as sound:')
display(Audio(fade(0.5 * (pulses - pulses.mean()) / np.max(np.abs(pulses))),
rate=AUDIO_RATE))
The pulse train, as sound:
The spectrum consists of multiple harmonics. We observe these spectral lines at the 250 Hz fundamental and its integer multiples (harmonics), with heights that decrease gradually as their frequency increases. The sharper the pulse (i.e., the shorter its duration), the larger the spectral lines at higher harmonics appear. In a milling cut, the fundamental frequency is the tooth passing frequency, set by the product of the spindle speed and the number of teeth, and the harmonics it produces are the healthy hum we previewed in Sect. A.2. The frequency domain treatment of milling in Sect. 4.3.3 of Schmitz and Smith [6] expands the periodic cutting force in exactly this way, as a Fourier series.
Part C applies the concepts we learned in Parts A and B. We give these harmonics their machining names, simulate the cutting process in the time domain following the straight tooth model from Sect. 4.4 of Schmitz and Smith [6], and listen to what happens when the cut stops being periodic. This occurs with chatter, a self-excited vibration with its own frequency that is not (in general) a multiple of the fundamental tooth passing frequency. The spectrum in Part C becomes the tooth passing family of frequencies for a cut. The spectrogram becomes the view in which chatter appears between the family frequency lines at the time that it occurs. Finally, we apply the tools examined in Part B to hear the sound of chatter in Part C.
▸ Key terms: pulse train
Summary¶
- Sound is a traveling pressure wave; a microphone responds to this pressure over time, and pitch and loudness correspond to the rate and size of the fluctuation.
- Sampling converts the continuous wave into a digital signal by sampling $f_s$ samples per second. Frequencies up to the Nyquist frequency $f_N = f_s/2$ (Eq. B.1) are represented accurately and content above it is aliased to a lower frequency unless an analog filter is applied before sampling occurs.
- A periodic signal can be written as a Fourier series (Eq. B.2), a sum of scaled harmonics at multiples of the fundamental frequency. The amplitude recipe is the signal's timbre and the same synthesis traces closed spatial curves as readily as temporal waveforms.
- The discrete Fourier transform (Eq. B.3) recovers the recipe from $N$ samples with resolution $\Delta f = 1/T$ (Eq. B.6), the FFT computes it quickly, and
compute_spectrumprovides this computation with a window and normalization. - The spectrogram stacks windowed short-time spectra frame by frame; the frame length sets a time-frequency tradeoff, with shorter frames more accurate in time and longer frames more accurate in frequency.
- The complete chain can be applied to any recording: the analysis window, sample rate, and frame length are the analyst's decisions, and Eqs. B.1 and B.6 predict the consequences of each decision.
- A repeating milling force is periodic, so its sound contains harmonics at integer multiples of its fundamental frequency. We analyze a stable milling cut's sound in Part C.
Exercises¶
1. A data acquisition system samples at $f_s = 8000$ per second. (a) What is its Nyquist frequency? (b) Tones at 1500 Hz, 5000 Hz, and 7200 Hz are recorded. Which alias, and at what frequency does each aliased signal appear? The aliased frequency is $|f - f_s \cdot \mathrm{round}(f / f_s)|$.
2. Using the explorer of Fig. B.3, find the smallest $K$ at which you can no longer hear a difference from the next step, and compare with a neighbor. (a) At $f_0 = 250$ Hz, what frequency is the highest harmonic you kept? (b) Human hearing ends near 20 kHz. What value of $K$ would place the highest harmonic at that limit, and what does that suggest about the audibility of the terms beyond your answer to (a)?
3. In the cell of Fig. B.4, change the recipe amplitudes from $1/k$ to $1/k^2$ for both recipes and rerun. (a) Predict, before running, how the sound will change. (b) Confirm against the spectrum, and describe the change in timbre in one sentence.
4. Two tones 20 Hz apart must appear as separate peaks in a spectrum. (a) Using Eq. B.6, what minimum record length is required? (b) For a spectrogram at $f_s = 48000$, what minimum frame length in samples is required, and what does that frame length cost in timing precision?
5. In the cell of Fig. B.5, change the text string to your own initials and rerun. (a) Which features of your letters resolve last as $K$ grows? (b) Using the discussion in the Corners and Clicks box, explain the pattern in one sentence.
6. Record a whistle that slides from low to high, load it in Sect. B.6, and analyze it.
(a) From the spectrogram, what frequency range does your whistle cover?
(b) Choose fs_new so that the top of your whistle lies above the new Nyquist frequency but the bottom lies below it. Predict, using the aliased-frequency formula from Exercise 1, where the folded portion of the line will appear, then rerun and check.
Appendix: Utility Functions¶
The code cells of this lesson rely on the helper functions, defined once in the collapsed cell at the top of the notebook (they must be defined before their first use so that Restart and Run All succeeds). They are identical to the functions used in Part A. For reference:
fade(x, ms, fs)applies a short raised-cosine fade-in and fade-out so audio clips start and stop without clicks.compute_spectrum(x, fs, fmax)returns the amplitude spectrum of a signal on a linear scale, normalized to a peak of one; Sect. B.4 walks through it line by line.compute_spectrogram(x, fs, nfft, hop, fmax)returns the short-time spectrogram on the same linear, normalized scale; Sect. B.5 describes its frame-by-frame operation.decimate_for_plot(t, x, max_points)passes plotting data through unchanged by default; setmax_pointsto an integer to thin a long record for display, keeping each bin's minimum and maximum so the visual envelope is preserved without aliasing. Computations always use the full-rate data.
References¶
[1] Nyquist, H. (1928). Certain topics in telegraph transmission theory. Transactions of the AIEE, 47(2), 617 to 644.
[2] Shannon, C. E. (1949). Communication in the presence of noise. Proceedings of the IRE, 37(1), 10 to 21.
[3] Fourier, J. B. J. (1822). Theorie analytique de la chaleur. Paris: Firmin Didot. English translation: The Analytical Theory of Heat, Cambridge University Press, 1878.
[4] Helmholtz, H. (1863). Die Lehre von den Tonempfindungen. English translation by A. J. Ellis from the 4th German edition: On the Sensations of Tone as a Physiological Basis for the Theory of Music (1877). Reprinted by Dover Publications, 1954.
[5] Cooley, J. W., and Tukey, J. W. (1965). An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19(90), 297 to 301. DOI: 10.1090/S0025-5718-1965-0178586-1
[6] Schmitz, T. L., and Smith, K. S. (2019). Machining Dynamics: Frequency Response to Improved Productivity, 2nd ed. Springer. DOI: 10.1007/978-3-319-93707-6
[7] Schmitz, T. L., and Gomez, M. F. Machining Dynamics: Jupyter Notebook Edition (in preparation).