Acoustic Signal Processing: Welch’s Method and SPL vs PSD Scaling
This article teaches signal processing for acoustics using Python.
It provides an interactive example using a realistic pressure signal from a propeller simulation and demonstrates producing Sound Pressure Level (SPL) and Power Spectral Density (PSD) spectra.
Phase-averaging is introduced to decompose the signal into deterministic (tonal) and nondeterministic (broadband) components.
Then, the “industry standard” Welch’s method is proposed as the best tool in practice.
It is explained using the equivalent case to the basic Fast Fourier Transform, then further advantages and nuances about scalings for tonal vs broadband noise are covered.
Open this Jupyter notebook in Colab to follow along and modify, or save for future use:
For a refresher on the fundamentals of signals and Fast Fourier Transforms, refer to the first article: FFT Signal Analysis using Python. This article will assume understanding of the topics presented there.
# We will use standard numpy and matplotlib, as well as welch from the scipy.signal library.
import numpy as np
from scipy.signal import welch
import matplotlib.pyplot as plt
1. The Raw Signal, and Tonal and Broadband noise
First the example pressure time-series data is loaded and plotted. It is the signal in Pascals (Pa) produced by a 5-bladed propeller measured at a location ~$5.9\times R_{propeller}$ from the center at an angle ~$60^{\circ}$ below the disk. The propeller is spinning at a constant 3806 RPM, so it is convenient to nondimensionalize time by the period of rotation. The signal includes a strong 5-per-revolution signal corresponding to the passage of each blade, called the Blade Passage Frequency (BPF), as well as broadband noise occurring across all frequencies. The repeated component is deterministic and can be called tonal noise. The random component is nondeterministic and is caused by random processes like turbulence. It is often referred to as broadband noise, the most common example being white noise. The human ear can perceive frequencies between 20 — 20,000 Hz, but is most sensitive to frequencies around 3,000 Hz, where human speech lies.
# Load time and pressure data
time, pressure = np.loadtxt("yaw70_U10_mic244_pressure.txt").T
rpm = 3806
period = 60 / rpm # 60 [s/min] * (1/rpm) [min/rev] = period [s/rev]
num_revs = 10
blade_pass_freq = rpm / 60 * 5 # (rpm) [rev/min] * (1/60) [min/s] * 5 [blade passages/rev] = BPF [Hz]
View plotting code
# Plot the original signal for a few revolutions
plt.figure(figsize=(10, 4))
plt.plot(time / period, pressure, color='royalblue', linewidth=1.5)
plt.xlim([0, num_revs])
plt.ylim([-0.2,0.3])
plt.xlabel("Revolutions", fontsize=FONT_SIZE_AXIS)
plt.ylabel("Pressure [Pa]", fontsize=FONT_SIZE_AXIS)
plt.xticks(np.arange(11), fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.grid(True)
plt.title("Original Pressure Signal", fontsize=FONT_SIZE_TITLE)
2. Phase-Averaging the Deterministic Component
The first step is phase-averaging to separate deterministic and nondeterministic noise. This is common for rotating machinery like propellers and aircraft turbofan engines. Time $t$ can be divided into strides of the period $T$ as:
\[t = \tau + n T\]where $n$ is the index of the period ($0,1, 2, … N-1$) and $\tau$ is the time within a single revolution ($0 \le \tau < T$). This allows the pressure $p(t)$ to be separated into the phase-averaged mean $\bar{p}(\tau)$ and residual $p’(t)$:
\[p(\tau + nT) = \bar{p}(\tau) + p'(\tau + nT)\]The phase-averaged signal $\bar{p}(\tau)$ is computed across all $N$ revolutions as:
\[\bar{p}(\tau) = \frac{1}{N} \sum_{n=0}^{N-1} p(\tau + nT)\]In Python, this is efficiently calculated by reshaping the 1D signal array into a 2D array of shape (num_revs, samples_per_rev) and taking the mean across the first axis.
The random broadband noise averages out, leaving only the tonal component that is consistently present across all revolutions.
The following plot shows the phase-averaged mean superimposed over each period of the original signal, where the revolution-to-revolution variance is the broadband noise.
# There are 10 revolutions of data with 720 samples each.
samples_per_rev = 720
num_revs = 10
# Reshape so each period is in a separate row
pressure_per_rev = pressure.reshape((num_revs, samples_per_rev))
# Phase average by taking the mean across the rows (revolutions)
pressure_phase_avg = np.mean(pressure_per_rev, axis=0)
View plotting code
# Plot all revolutions superimposed, followed by the phase average
plt.figure(figsize=(10, 4))
time_one_rev = time[:samples_per_rev] / period
# Plot individual revolutions with a colormap
colors = sbn.color_palette("colorblind", num_revs)
lines = []
for i in range(num_revs):
line, = plt.plot(time_one_rev, pressure_per_rev[i, :], color=colors[i])
lines.append(line)
# Plot phase-averaged signal
line_avg, = plt.plot(time_one_rev, pressure_phase_avg, color='black', linewidth=2.5)
plt.xlim([0, 1])
plt.ylim([-0.2,0.3])
plt.xlabel("Revolutions", fontsize=FONT_SIZE_AXIS)
plt.ylabel("Pressure [Pa]", fontsize=FONT_SIZE_AXIS)
plt.xticks(fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.grid(True)
# Temporarily thicken lines so the legend picks it up
for line in lines:
line.set_linewidth(3.0)
plt.legend([tuple(lines), line_avg], ['Individual Revs', 'Phase Average'],
handler_map={tuple: HandlerTuple(ndivide=None)}, loc="upper right",
fontsize=FONT_SIZE_LEGEND, ncols=2)
# Reset lines to their original thin width
for line in lines:
line.set_linewidth(1.0)
plt.title("Phase-Averaging of Pressure Signal", fontsize=FONT_SIZE_TITLE)
3. Deterministic and Nondeterministic Decomposition
The nondeterministic component can be obtained by subtracting the phase-averaged deterministic component from the total original signal.
\[p'(t) = p(\tau + nT) - \bar{p}(\tau)\]This is done in code by using numpy.tile to repeat the phase averaged pressure by the number of revs to match the original signal length, and subtracting it from the original signal. The equivalent MATLAB function is repmat. The following figure shows one period of the original signal and the separated deterministic and nondeterministic components.
# Subtract per-rev average out of original pressure to get broadband component
pressure_broadband = pressure - np.tile(pressure_phase_avg, num_revs)
View plotting code
plt.figure(figsize=(10, 4))
plt.plot(time_one_rev, pressure[:samples_per_rev], color='black', linewidth=2.5, label='Total Signal')
plt.plot(time_one_rev, pressure_phase_avg, color='crimson', linestyle='dashed', linewidth=2, label="Deterministic (Tonal)")
plt.plot(time_one_rev, pressure_broadband[:samples_per_rev], color='royalblue', linestyle='dotted', linewidth=2, label="Nondeterministic (Broadband)")
plt.xlim([0, 1])
plt.ylim([-0.2,0.3])
plt.xlabel("Revolutions", fontsize=FONT_SIZE_AXIS)
plt.ylabel("Pressure [Pa]", fontsize=FONT_SIZE_AXIS)
plt.xticks(fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.grid(True)
plt.legend(loc="upper right", fontsize=FONT_SIZE_LEGEND, ncols=2)
plt.title("Signal Decomposition (One Revolution)", fontsize=FONT_SIZE_TITLE)
4. Frequency Domain: FFT vs Welch’s Method
Noise analysis ultimately requires examination in the frequency domain, and in practice, engineers usually use Welch’s method. This approach can be explained by introducing the simple case in which Welch’s method is equivalent to a standard Fast Fourier Transform (FFT). (For a refresher on FFT fundamentals, refer to the article: FFT Signal Analysis using Python)
First, a one-sided FFT is computed manually to obtain the spectrum of the total signal:
num_samples = len(pressure)
dt = time[1] - time[0]
# Manual FFT
fft_spectrum = np.fft.fft(pressure)
fft_freqs = np.fft.fftfreq(num_samples, d=dt)
# The one-sided spectrum includes DC (0 Hz) and the Nyquist frequency (if N is even).
# To convert two-sided spectrum to one-sided, first slice away negative frequencies and take out the scaling factor
end_idx_including_nyquist = num_samples // 2 + 1
fft_onesided = np.abs(fft_spectrum[:end_idx_including_nyquist]) / num_samples
fft_freqs = fft_freqs[:end_idx_including_nyquist]
# Next, multiply the positive frequencies by 2 (excluding DC (0) and Nyquist frequency).
# Even-length signals include the Nyquist frequency, which should not be doubled. Odd-length signals do not.
# Using ceiling function captures the correct stop index to double, excluding Nyquist if it exists.
end_idx_excluding_nyquist = int(np.ceil(num_samples / 2))
fft_onesided[1:end_idx_excluding_nyquist] *= 2
Welch’s method is an approach used to estimate spectra of random data. For rigorous derivation of these formulas, refer to Bendat and Piersol’s Random Data: Analysis and Measurement Procedures. The variance of a signal of length $N$ can be reduced by splitting it into $K$ overlapping segments of length $M$ and averaging the spectrum of each segment. Additionally, a window $w(n)$ can be multiplied to each segment to reduce spectral leakage. Another crucial detail is that Welch’s method directly computes the mean-square spectrum, also called the autospectrum. The one-sided autospectrum using FFT via Welch’s method is:
\[P_{xx}(f) = \frac{1}{2 K} \sum_{k=1}^{K} \left| \frac{2 U}{M} \sum_{n=0}^{M-1} p_k(n) w(n) e^{-j 2\pi f n} \right| ^2 \tag{1}\]where $p_k(n)$ is the time-domain pressure signal of segment $K$ and $U$ is a normalization factor which compensates for the amplitude reduction caused by applying the window. This value will be discussed in the next section. Using a ‘boxcar’ window ($w(n)=1$ so $U=1$) and a single segment ($M=N$ so $K=1$) simplifies Welch’s equation to:
\[P_{xx}(f) = \frac{1}{2} |X(f)|^2 = \overline{p^2}\]where $X(f)$ represents the normalized and unwindowed one-sided FFT of the entire signal (exactly the fft_onesided array computed manually above).
Because this FFT algorithm produces a spectrum of the amplitude of sine waves, the factor of one half is necessary to calculate the root-mean-square.
# FFT output is amplitude of sine waves. The RMS of a sine wave is A / sqrt(2)
# Squaring it yields the Mean Square
fft_rms_sq = 0.5 * fft_onesided**2
# Mathematically equivalent operation using Welch's method (scaling='spectrum' for tonal amplitudes)
welch_freqs, welch_spectrum = welch(pressure, fs=1/dt, window='boxcar', nperseg=num_samples, scaling='spectrum')
Finally, the spectra are converted to Sound Pressure Level (SPL) in decibels. This is a logarithmic scale contextualized to human perception of noise. It uses a specific reference pressure $p_{ref} = 20 \mu\text{Pa}$ corresponding to the absolute threshold of human hearing, thus any noise with SPL below $0$ dB is too quiet to be heard by humans.
\[SPL = 10 \log_{10} \left( \frac{\overline{p^2}}{p_{ref}^2} \right) \quad [\text{dB}]\]When these steps are taken, identical SPL spectra are produced by the manual FFT and Welch’s algorithm, as demonstrated by the figure below.
# Convert to decibels (Sound Pressure Level)
fft_spl = 10 * np.log10(fft_rms_sq / (2e-5)**2)
welch_spl = 10 * np.log10(welch_spectrum / (2e-5)**2)
View plotting code
plt.figure(figsize=(10, 5))
plt.semilogx(fft_freqs, fft_spl, color='black', linewidth=2.5, label="Manual FFT")
plt.plot(welch_freqs, welch_spl, color='crimson', linestyle='dotted', linewidth=2, label="Welch Method")
plt.ylim(0, 75)
plt.xlim(100, 20000)
plt.xlabel("Frequency [Hz]", fontsize=FONT_SIZE_AXIS)
plt.ylabel("SPL [dB]", fontsize=FONT_SIZE_AXIS)
plt.xticks(fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.grid(True)
plt.legend(loc="upper right", fontsize=FONT_SIZE_LEGEND, ncols=2)
plt.title("Equivalence between FFT and Welch's Method", fontsize=FONT_SIZE_TITLE)
5. Window Correction Factors (Amplitude vs Energy)
From the Welch equation in Section 4 (Equation 1), the normalization factor $U$ compensates for the energy reduction caused by the window function. However, the correction depends on whether the noise is tonal or broadband. Analysis of tones that occur at exact frequencies requires the preservation of peak amplitudes, while analysis of continuous broadband spectra requires preservation of total energy.
SciPy’s welch function encodes this distinction in the scaling argument. This corresponds directly to MATLAB’s spectrumType argument for the pwelch function.
Python (SciPy) welch |
MATLAB pwelch |
|
|---|---|---|
| Argument | scaling |
spectrumType |
| Power Spectral Density | "density" (default) |
"psd" (default) |
| Power Spectrum | "spectrum" |
"power" |
Correction in SciPy Welch Function
From scipy/signal/_spectral_py.py (v1.18.0):
if scaling == 'density':
scale = 1.0 / (fs * (win*win).sum())
elif scaling == 'spectrum':
scale = 1.0 / win.sum()**2
In the SciPy implementation, the scale variable also includes the normalization by the segment length $M$ and is moved outside of the squared term as written in Equation 1.
The divisor $\mathtt{fs} = M / T$ used in density scaling introduces the frequency bin width $1/T$ to achieve a density spectrum.
For density scaling:
and for spectrum scaling:
The correction factors for each scaling method are derived below using a Hanning window for example.
Energy Correction — scaling='density' (Broadband Noise)
For broadband noise, energy is spread across all frequencies. The goal is to preserve the total energy, or variance of the signal. When frequency contributions are independent, their energies add incoherently, so the total energy across all bins scales with $\sum_n w(n)^2$. For a Hanning window, this sum evaluates to $3M/8$. Therefore, $U_{density} = \sqrt{\frac{M}{3M/8}} = \sqrt{8/3} \approx \mathbf{1.63}$.
Amplitude Correction — scaling='spectrum' (Tonal Noise)
For tonal noise, energy is concentrated at discrete frequencies. The goal is to recover the true peak amplitude at each tone. A tone oscillates at a single frequency $f_0$, so all $M$ windowed samples add coherently in the Fourier-transformed sum. The resulting magnitude $A$ at a single bin is proportional to the sum of the window weights: \(A \sum_{n=0}^{M-1} w(n)\). For a Hanning window, this sum evaluates to $M/2$, therefore $U_{spectrum} = \frac{M}{M/2} = \mathbf{2.0}$.
Summary
scaling |
When to use | What it preserves | Correction | Hanning Correction |
|---|---|---|---|---|
'density' |
Analyzing broadband/continuous noise | Total integrated energy | $\overline{w(n)^2}^{-1/2}$ | $\approx 1.63$ |
'spectrum' |
Analyzing discrete tones (e.g., BPF harmonics) | Peak amplitude of discrete tones | $\overline{w(n)}^{-1}$ | $2.0$ |
Because these correction factors differ ($2.0 \neq 1.63$), the two scaling options are not interchangeable.
Even if a density scaling output is multiplied by the bin width $\Delta f$ to recover SPL, the result will
still underpredict the peak amplitude of a discrete tone by $20\log_{10}(2.0/1.63) \approx 1.76\ \text{dB}$.
This is demonstrated in the figure below.
freqs_tonal, spec_tonal_no_window = welch(pressure_phase_avg, fs=1/dt, window='boxcar',
nperseg=samples_per_rev, scaling='spectrum')
_, spec_tonal1 = welch(pressure_phase_avg, fs=1/dt, window='hann',
nperseg=samples_per_rev, scaling='spectrum')
_, spec_tonal2 = welch(pressure_phase_avg, fs=1/dt, window='hann',
nperseg=samples_per_rev, scaling='density')
spl_tonal_no_window = 10 * np.log10(spec_tonal_no_window / (2e-5)**2)
spl_tonal_spectrum = 10 * np.log10(spec_tonal1 / (2e-5)**2)
# Multiply density (V^2/Hz) by df (Hz) to get power (V^2), then convert to SPL
df = freqs_tonal[1]
spl_tonal_density_integrated = 10 * np.log10(spec_tonal2 * df / (2e-5)**2)
View plotting code
plt.figure(figsize=(10, 5))
plt.plot(freqs_tonal / blade_pass_freq, spl_tonal_no_window, color='black', linestyle='dotted',
linewidth=2, label="No Window")
plt.plot(freqs_tonal / blade_pass_freq, spl_tonal_spectrum, color='crimson',
linewidth=2.5, label="Hann, scaling='spectrum'")
plt.plot(freqs_tonal / blade_pass_freq, spl_tonal_density_integrated, color='royalblue',
linewidth=2.5, linestyle='dashed', label="Hann, scaling='density' * df")
plt.xlim((0.5, 3.5))
plt.ylim((20, 75))
plt.xticks(np.arange(1, 4), fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.xlabel("Frequency / BPF", fontsize=FONT_SIZE_AXIS)
plt.ylabel("SPL [dB]", fontsize=FONT_SIZE_AXIS)
plt.title("Non-Equivalence from Correction in SciPy Welch's Spectrum vs Density Scaling", fontsize=FONT_SIZE_TITLE)
plt.grid(True)
plt.legend(loc="upper right", fontsize=FONT_SIZE_LEGEND)
plt.show()
6. Comparing Tonal Extract to Total Signal
To analyze rotating machinery it is common to look at multiples of the Blade Passage Frequency (BPF) $f \times N_{blades}$, which will be the frequencies of the tonal noise peaks. The Sound Pressure Level (SPL) metric is used since the amplitude of the discrete tones is of interest.
Below, the SPL of the total signal’s $p(t)$ is compared to the peaks of the SPL of the phase-averaged tonal noise $\bar{p}(\tau)$.
freqs_tonal, spec_tonal = welch(pressure_phase_avg, fs=1/dt, window='boxcar',
nperseg=samples_per_rev, scaling='spectrum')
spl_tonal = 10 * np.log10(spec_tonal / (2e-5)**2)
# Extract just the peaks at the BPF harmonics
# BPF is 5 times the fundamental rotation frequency. The frequency resolution df is exactly the rotation frequency.
# Therefore, the BPF harmonics are located at indices 5, 10, 15, 20, 25...
bpf_indices = np.arange(5, len(freqs_tonal), 5)
bpf_harmonics = freqs_tonal[bpf_indices]
spl_bpf_peaks = spl_tonal[bpf_indices]
View plotting code
plt.figure(figsize=(10, 5))
plt.plot(welch_freqs / blade_pass_freq, welch_spl, color='black', linestyle='dotted', linewidth=2, label='Total Signal')
plt.scatter(bpf_harmonics / blade_pass_freq, spl_bpf_peaks, c='None', edgecolors='crimson', s=80, linewidths=2, zorder=5,
label='Phase Average Peaks (Tonal)',marker='d')
plt.ylim(0, 80)
plt.xlim(0.5, 5.5)
plt.xticks(np.arange(1, 6), fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.xlabel("Frequency / BPF", fontsize=FONT_SIZE_AXIS)
plt.ylabel("SPL [dB]", fontsize=FONT_SIZE_AXIS)
plt.grid(True)
plt.legend(loc="upper right", fontsize=FONT_SIZE_LEGEND, ncols=2)
plt.title("Tonal Peaks vs Total Signal", fontsize=FONT_SIZE_TITLE)
7. Comparing Broadband Extract to Total Signal
Because broadband noise occurs at continuous frequencies across the spectrum, Power Spectral Density (PSD) metric is used, which describes the power contained within each frequency bin. If SPL is used, which does not normalize by the bin width, amplitudes will vary depending on the frequency resolution. The independence from frequency resolution makes PSD ideal for displaying continuous broadband noise.
\[PSD = 10 \log_{10} \left( \frac{\overline{p^2} }{\Delta f~p_{ref}^2} \right) \quad [\text{dB/Hz}]\]Welch’s Method Features that are Convenient for Broadband Noise Analysis
Unlike the Welch calculation for the tonal spectrum which used a single segment of the entire signal length, the broadband analysis uses the full capability of Welch’s method by splitting the signal into multiple overlapping segments. This was not necessary for the tonal noise since the averaging was performed in the time domain.
When the signal is split into $K$ segments, each segment produces an independent estimate of the PSD. Averaging these $K$ estimates reduces the variance of the result. SciPy defaults to 50% overlap between segments, which nearly doubles $K$ for a given record length. A Hanning window is also applied to each segment to reduce spectral leakage from the broadband energy. The figure below demonstrates the effect of averaging these spectra on reducing the variance of the broadband signal $p’$.
# Compute the PSD for every individual overlapping segment
freqs_seg, times_seg, psd_segments_raw = spectrogram(pressure_broadband, fs=1/dt, window='hann',
nperseg=samples_per_rev, noverlap=samples_per_rev//2,
scaling='density')
psd_segments = 10 * np.log10(psd_segments_raw / (2e-5)**2)
# The Welch estimate is simply the average of these individual segment PSDs
psd_welch_avg = 10 * np.log10(np.mean(psd_segments_raw, axis=-1) / (2e-5)**2)
K_segments = psd_segments.shape[1]
View plotting code
plt.figure(figsize=(10, 5))
# Plot all individual segments using a colormap
colors = sbn.color_palette("colorblind", K_segments)
lines = []
for i in range(K_segments):
line, = plt.semilogx(freqs_seg, psd_segments[:, i], color=colors[i])
lines.append(line)
# Plot the Welch averaged PSD
line_avg, = plt.semilogx(freqs_seg, psd_welch_avg, color='black', linewidth=2.5)
plt.ylim(0, 70)
plt.xlim(100, 20000)
plt.xlabel("Frequency [Hz]", fontsize=FONT_SIZE_AXIS)
plt.ylabel("PSD [dB/Hz]", fontsize=FONT_SIZE_AXIS)
plt.xticks(fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.grid(True)
# Temporarily thicken lines so the legend picks it up
for line in lines:
line.set_linewidth(3.0)
plt.legend([tuple(lines), line_avg], [f'Spectra of Individual Segments ($K={K_segments}$)', 'Average Spectrum'],
handler_map={tuple: HandlerTuple(ndivide=None)}, loc="upper right",
fontsize=FONT_SIZE_LEGEND)
# Reset lines to their original thin width
for line in lines:
line.set_linewidth(1.0)
plt.title("Variance Reduction by Welch's Method using Multiple Segments", fontsize=FONT_SIZE_TITLE)
Finally, the averaged PSD of broadband $p’$ using Welch’s method is compared to the total signal in the figure below. The peaks at the BPF harmonics have been removed.
# Welch's method (scaling='density' for Power Spectral Density)
welch_freqs, welch_psd_raw = welch(pressure, fs=1/dt, window='boxcar', nperseg=num_samples, scaling='density')
welch_psd = 10 * np.log10(welch_psd_raw / (2e-5)**2)
freqs_bb, spec_broadband = welch(pressure_broadband, fs=1/dt, window='hann', nperseg=samples_per_rev, scaling='density')
psd_broadband = 10 * np.log10(spec_broadband / (2e-5)**2)
View plotting code
plt.figure(figsize=(10, 5))
plt.semilogx(welch_freqs, welch_psd, color='black', linestyle='dotted', linewidth=2, label='Total Signal')
plt.semilogx(freqs_bb, psd_broadband, color='royalblue', linewidth=2.5, label='Broadband')
plt.ylim(0, 70)
plt.xlim(100, 20000)
plt.xlabel("Frequency [Hz]", fontsize=FONT_SIZE_AXIS)
plt.ylabel("PSD [dB/Hz]", fontsize=FONT_SIZE_AXIS)
plt.xticks(fontsize=FONT_SIZE_TICKS)
plt.yticks(fontsize=FONT_SIZE_TICKS)
plt.grid(True)
plt.legend(loc="upper right", fontsize=FONT_SIZE_LEGEND, ncols=2)
plt.title("Broadband vs Total Signal", fontsize=FONT_SIZE_TITLE)