Spectral Estimation and Time–Frequency Analysis

Introduction

The Fourier transform provides a powerful description of signal content in frequency, but practical digital signals are finite, noisy, sampled, and often nonstationary. Spectral estimation addresses the problem of inferring the distribution of signal power over frequency from observations. Time–frequency analysis extends this objective by describing how spectral content changes with time.

These methods are central to radar and sonar, vibration monitoring, biomedical instrumentation, speech processing, wireless communications, power-system analysis, and machine learning for signal classification. At postgraduate level, it is important not only to calculate a spectrum but also to understand estimator bias, variance, resolution, leakage, statistical confidence, and the assumptions behind stationarity.

Learning Objectives

After studying this chapter, the reader should be able to:

1. distinguish energy spectra, power spectral density, and finite-record spectral estimates;
2. derive and interpret the periodogram and averaged spectral estimators;
3. explain spectral leakage, windowing, resolution, bias, variance, and aliasing;
4. apply Welch, Bartlett, Blackman–Tukey, parametric, and multitaper methods;
5. formulate and interpret short-time Fourier transform and spectrograms;
6. compare wavelet and quadratic time–frequency representations;
7. select suitable methods for stationary and nonstationary signals;
8. evaluate results critically using quantitative and statistical criteria.

1. Spectral Concepts and Statistical Assumptions

For a deterministic finite-energy sequence x[n], its discrete-time Fourier transform (DTFT) is

where ω is angular frequency in radians per sample. The energy spectral density is

For a random, wide-sense stationary process x[n], the autocorrelation is

where E{·} denotes expectation and * denotes complex conjugation. The power spectral density (PSD) is the Fourier transform of the autocorrelation:

This is the discrete-time Wiener–Khinchin relation. The PSD is nonnegative and periodic with period 2π. For a real-valued process, it is even: S_x(e^{jω}) = S_x(e^{-jω}).

If the sampling frequency is f_s, physical frequency f and digital frequency are related by ω = 2πf/f_s. A one-sided PSD for a real signal is commonly obtained by doubling the positive-frequency values, except at direct current and the Nyquist frequency. This convention preserves total mean-square value:

The equality assumes wide-sense stationarity and finite average power. PSD units are signal-unit² per hertz when expressed against f, or signal-unit² per radian/sample when expressed against ω.

Terminology

• Stationarity: Statistical properties do not change with time; wide-sense stationarity requires constant mean and autocorrelation depending only on lag.
• Ergodicity: Time averages converge to ensemble averages under suitable conditions.
• Spectral line: Concentrated power associated with a sinusoidal or periodic component.
• Broadband component: Power spread over a frequency interval.
• Spectral estimate: A statistic computed from finite observations to approximate an unknown spectrum.
• Resolution: The ability to distinguish nearby spectral components.
• Leakage: Redistribution of energy into frequencies other than the component's actual frequency.

2. The Periodogram

Given N observations x[0], …, x[N−1], the finite Fourier transform is

The modified periodogram with window w[n] is

where

and f_s is the sampling frequency. This normalization produces PSD units per hertz for a windowed sequence. The periodogram is generally asymptotically unbiased under suitable regularity conditions, but it is not a consistent estimator: increasing N does not cause its variance at a fixed frequency to approach zero. Its variance is approximately comparable to the square of the underlying PSD, even for long records.

For computation, the discrete frequencies are f_k = kf_s/N, k = 0, …, N−1. Zero-padding increases the density of plotted frequency samples but does not increase true resolution or information.

Spectral Leakage and Windows

A finite record is equivalent to multiplying an infinite signal by a rectangular window. In frequency, multiplication becomes convolution with the window spectrum. A sinusoid whose frequency is not exactly aligned with a DFT bin therefore spreads over nearby bins.

Window selection is a trade-off:

The main-lobe width controls separation of nearby tones; sidelobes control contamination from strong components. No window simultaneously gives the narrowest main lobe and the lowest sidelobes.

A useful frequency-resolution approximation is Δf ≈ f_s/N for a rectangular-window DFT grid, but the effective resolution is determined by the window's main-lobe width. For a Hann window, the ability to distinguish tones is poorer than the bin spacing alone suggests.

3. Averaged Nonparametric Estimators

Bartlett and Welch Methods

The Bartlett method divides a record into K non-overlapping segments, computes a periodogram for each, and averages them:

Averaging reduces variance but shortens the segments and therefore reduces frequency resolution.

Welch's method generalizes Bartlett's approach by applying a window and allowing overlap. If segments are sufficiently independent, averaging approximately K periodograms gives an effective increase in degrees of freedom. With overlap, adjacent estimates are correlated, so the improvement is less than that predicted by K independent segments.

Welch estimation is often the best default for engineering PSD measurement because it provides a practical compromise between leakage, variance, and computational cost. Segment length, overlap, window, detrending, and scaling must be reported; otherwise results are difficult to reproduce.

Blackman–Tukey Estimation

The sample autocorrelation is first estimated and then windowed:

where v[k] is a lag window and M is the retained correlation length. Truncating the autocorrelation improves variance but introduces frequency smoothing. Blackman–Tukey methods are useful when correlation structure is of primary interest.

Parametric Estimation

Parametric methods assume that the data arise from a model with a small number of parameters. An autoregressive AR(p) model is

where e[n] is zero-mean white innovation noise with variance σ_e². Its PSD is

The sign convention follows the model above. AR spectra can provide high resolution for short records, especially for sharp resonances, but an incorrect model order can create artificial peaks or suppress genuine structure. The Yule–Walker equations estimate AR coefficients from autocorrelation values:

with the corresponding prediction-error equation at m = 0. Model validation is essential.

The moving-average and ARMA models extend this framework. For ARMA(p,q),

where A(z) and B(z) are denominator and numerator polynomials. Parametric estimates are efficient when the model is correct but biased when it is not.

Multitaper Estimation

Multitaper methods use several orthogonal tapers, commonly discrete prolate spheroidal sequences, to produce multiple approximately independent spectral estimates. Their weighted average reduces variance while controlling leakage. A generic form is

where X_r is the tapered transform and μ_r are nonnegative weights. Multitaper analysis is valuable for short records, closely spaced components, and signals with strong spectral dynamic range. The time–bandwidth product and number of tapers determine the resolution–variance compromise.

4. Statistical Interpretation and Confidence

A periodogram ordinate should not be interpreted as an exact spectrum. For a large sample and under approximate Gaussian assumptions, an averaged PSD estimate can be represented approximately as a scaled chi-square random variable. If ν degrees of freedom apply,

where S(f) is the true PSD and χ²_ν denotes a chi-square variable with ν degrees of freedom. Thus a two-sided 100(1−α)% confidence interval is approximately

where χ²_{ν,p} is the p-quantile. The approximation is frequency dependent and can fail near strong deterministic lines, boundaries, or for strongly non-Gaussian data. Confidence intervals should therefore be accompanied by the estimator, effective degrees of freedom, and assumptions.

5. Cross-Spectral Analysis

For two jointly stationary processes x[n] and y[n], the cross-spectrum is

where r_xy[k] = E{x[n]y*[n−k]}. The magnitude-squared coherence is

and satisfies 0 ≤ C_xy(f) ≤ 1 when spectra are properly defined. High coherence indicates a stable linear relationship at a frequency, not necessarily causation. The cross-spectrum phase can estimate phase delay; for a pure delay τ, phase is approximately −2πfτ, subject to phase wrapping and model assumptions.

Cross-spectral methods support modal analysis, transfer-function estimation, noise-source identification, and system diagnostics.

6. Time–Frequency Analysis

A stationary spectrum averages information over the entire record. When frequency content evolves, this averaging can hide important events. Time–frequency methods represent signal energy jointly in time and frequency, necessarily accepting a localization trade-off: a narrow window improves frequency localization but worsens time localization, while a short window does the reverse.

Short-Time Fourier Transform

The STFT of x[n] is

where w[n] is an analysis window, m is the frame index, and R is the hop size in samples. The spectrogram is

possibly with a normalization factor for PSD units. A smaller window provides better time localization; a larger window provides narrower frequency resolution. The uncertainty principle prevents simultaneous arbitrary precision in both domains.

For a sampling rate f_s, window length L samples corresponds approximately to a time span L/f_s. If the window's effective bandwidth is B_w, then frequency detail is on the order of B_w, not merely f_s/L. Overlap reduces visual discontinuity and improves time sampling but does not independently create frequency information.

Reassignment and Synchrosqueezing

Spectrogram energy is often smeared because each window has finite support. Reassignment estimates local instantaneous frequency and relocates energy toward sharper coordinates. Synchrosqueezing similarly concentrates selected components while retaining an invertible or approximately invertible representation under suitable conditions. These methods are useful for separating ridges, but they are sensitive to noise, component crossings, and parameter choices.

Wavelet Analysis

The continuous wavelet transform is

where ψ(t) is the mother wavelet, a is scale, b is time shift, and * denotes conjugation. Small scales correspond broadly to high frequencies and large scales to low frequencies. Wavelets provide multiresolution analysis: short effective windows at high frequencies and long windows at low frequencies. This is advantageous for transients, discontinuities, and signals with scale-dependent structure.

Discrete wavelet transforms use filter banks. At each level, a low-pass filter produces approximation coefficients and a high-pass filter produces detail coefficients, followed by downsampling. Perfect reconstruction requires compatible analysis and synthesis filters and appropriate boundary treatment.

Wigner–Ville and Quadratic Distributions

The discrete Wigner–Ville distribution is conceptually based on

with suitable discrete implementation conventions. It has excellent concentration for single components but produces cross-terms between multiple components. Smoothed pseudo-Wigner–Ville distributions reduce cross-terms at the cost of resolution. Quadratic methods should therefore be used when component interaction and interpretation are understood.

7. Practical Method Selection

A defensible workflow is: inspect and detrend the data; verify sampling and anti-aliasing; identify stationarity; select window and segment length; calculate an estimate with documented scaling; inspect confidence or repeatability; compare with an alternative method; and validate conclusions in the time domain or against physical knowledge.

8. Worked Example: Welch PSD

Consider a vibration signal sampled at f_s = 1,000 Hz with N = 8,000 samples. Use Hann-windowed Welch estimation with segment length L = 1,000 samples and 50% overlap. The hop is R = 500 samples. The number of complete segments is

The DFT frequency spacing for each segment is

This is the frequency sampling interval, not necessarily the ability to resolve two tones. A Hann window has a broader main lobe than a rectangular window, so practical separation requires more than 1 Hz in many situations. Averaging 15 modified periodograms substantially reduces random fluctuations. Because overlapping Hann segments are correlated, the effective degrees of freedom are less than 30, the value that would result from 15 independent real-valued periodograms. A 20-Hz peak recurring across segments is therefore credible as a persistent vibration component, whereas a single isolated bin should not be overinterpreted.

9. Applied Case: Bearing-Fault Monitoring

Background

An industrial motor is monitored using an accelerometer sampled at 12.8 kHz. Bearing defects generate periodic impulsive excitations, often modulated by shaft rotation and resonance bands.

Problem

A conventional whole-record FFT shows a broad resonance but does not reveal whether impulses occur continuously or only during particular operating conditions.

Analysis and Approach

The signal is first high-pass filtered to remove slow structural drift, with filter transients excluded. A Welch PSD identifies the resonance band and estimates its stable power. An STFT using a window selected to contain several expected fault cycles tracks changes over time. The squared-envelope spectrum is then calculated after band-pass filtering around the resonance. Because envelope analysis can amplify noise and be affected by filter choice, coherence with a tachometer channel is also examined when available.

If the machine speed varies substantially, fixed-frequency interpretation is inadequate. Order tracking, which maps frequency to rotational order, can reduce smearing. Time–frequency ridges are compared with known bearing characteristic frequencies, but diagnosis is not based on one peak alone; load, speed, sensor location, and baseline measurements are considered.

Outcome

A persistent resonance with regularly recurring envelope impulses supports a bearing-related hypothesis. A transient ridge appearing only during load changes instead suggests a different mechanical event. Maintenance action is scheduled only after confirming the pattern across repeated runs and checking mechanical inspection data.

Learning Points

1. Global PSD and local time–frequency analysis answer different questions.
2. Preprocessing and operating-speed variation can dominate interpretation.
3. Cross-channel evidence improves confidence.
4. A statistically visible component is not automatically a mechanically causal fault.

10. Limitations, Governance, and Contemporary Developments

Spectral estimates can be distorted by aliasing, missing samples, clipping, nonstationarity, leakage, boundary effects, and inappropriate normalization. Anti-alias filtering must precede sampling; no later algorithm can recover frequencies already aliased. In time–frequency analysis, edge regions are less reliable because the analysis window extends beyond available data. Padding can improve appearance but does not remove the underlying uncertainty.

In safety-critical monitoring, thresholds should be validated on representative operating conditions and documented with false-alarm and missed-detection rates. Data provenance, calibration records, timestamp integrity, and secure storage are important governance requirements. Automated or AI-assisted spectral classifiers should expose quality indicators and confidence measures rather than presenting unexplained labels. Training data should include variation in speed, load, sensor placement, and normal degradation. Privacy may matter when acoustic or biomedical signals contain identifying information.

Contemporary research includes adaptive and Bayesian spectral estimation, robust methods for impulsive noise, graph and spatial spectra for sensor networks, differentiable time–frequency layers for deep learning, and learned transforms. These methods can improve performance but do not eliminate the need for physical interpretation, statistical validation, and reproducibility.

Best-Practice Checklist

• State whether the signal is treated as stationary.
• Report sampling rate, record length, detrending, window, segment length, overlap, and scaling.
• Distinguish DFT bin spacing from true spectral resolution.
• Use PSD units and verify total integrated power against time-domain mean square.
• Avoid interpreting zero-padding as added resolution.
• Quantify uncertainty or repeatability where decisions depend on spectral peaks.
• Compare a nonparametric estimate with a physically justified model when appropriate.
• Treat edge effects, leakage, aliasing, and nonstationarity explicitly.
• Validate automated conclusions against independent measurements.

Chapter Summary

Spectral estimation infers power distribution from finite observations. The periodogram is simple but has high variance; windowing controls leakage at the cost of resolution. Bartlett and Welch averaging reduce variance, while Blackman–Tukey methods smooth estimated autocorrelation. Parametric AR and ARMA methods can provide high resolution under correct modeling assumptions, and multitaper methods offer strong leakage control with useful statistical stability.

A stationary PSD cannot show when a component occurs. The STFT and spectrogram provide joint time–frequency information through a window-dependent resolution trade-off. Wavelets provide multiresolution behavior, reassignment concentrates smeared energy, and quadratic distributions can provide sharp localization but may create cross-terms. Cross-spectral coherence supports system diagnosis but does not establish causation. Reliable practice requires correct sampling, documented parameters, uncertainty awareness, physical validation, and careful interpretation of automated results.

Key Terms

Spectral estimation: Statistical inference of a signal's frequency-power distribution.

Power spectral density: Average power per unit frequency of a stationary process.

Periodogram: Squared magnitude of a finite-record Fourier transform with scaling.

Spectral leakage: Spreading of spectral energy caused by finite observation and windowing.

Window function: Weight sequence applied before Fourier transformation.

Resolution: Ability to distinguish nearby spectral components.

Welch method: Averaged, windowed, usually overlapping periodogram estimator.

Blackman–Tukey method: Spectrum estimated from a windowed sample autocorrelation.

Parametric estimator: Spectrum derived from an assumed stochastic model.

Autoregressive model: Model in which the current sample depends on previous samples and innovation noise.

Multitaper method: Estimator averaging spectra obtained with orthogonal tapers.

Cross-spectrum: Frequency-domain relation between two jointly stationary signals.

Coherence: Normalized squared magnitude of the cross-spectrum.

STFT: Fourier transform calculated over successive localized windows.

Spectrogram: Squared magnitude of an STFT.

Wavelet transform: Scale- and shift-based multiresolution signal representation.

Reassignment: Relocation of time–frequency energy toward estimated local coordinates.

Wigner–Ville distribution: Quadratic time–frequency representation with high concentration and possible cross-terms.

Aliasing: Irreversible frequency overlap caused by insufficient sampling.

Stationarity: Invariance of relevant statistical properties over time.

Review Questions

Short Answer Questions

1. Why is the periodogram not a consistent estimator of the PSD?
2. What is the purpose of a window in finite-record spectral analysis?
3. Distinguish DFT frequency spacing from actual spectral resolution.
4. Define magnitude-squared coherence and state its interpretation.
5. Why does the STFT involve a time–frequency resolution trade-off?

Long Answer Questions

1. Derive the PSD of an AR(p) process and discuss the consequences of model-order selection.
2. Compare Bartlett, Welch, Blackman–Tukey, and multitaper estimation in terms of bias, variance, and resolution.
3. Explain spectral leakage using the time-domain multiplication and frequency-domain convolution relationship.
4. Compare STFT, wavelet, and Wigner–Ville representations for nonstationary signals.
5. Develop a complete, reproducible workflow for vibration PSD estimation, including uncertainty and validation.

Analytical / Application Questions

1. A 2,048-sample record is sampled at 8 kHz. Determine the DFT frequency spacing and explain why it is not sufficient by itself to specify resolution.
2. A machine signal has a stable resonance but intermittent impulses. Propose a PSD and time–frequency analysis plan, including preprocessing and validation.
3. Two sensors show high coherence only near one frequency band. Explain what conclusions are justified, what conclusions are not justified, and what additional measurements would improve diagnosis.

Multiple Choice Questions with Answers

Question 1

Which statement best describes the periodogram for a stationary random process?

A. It has zero variance for every finite record.

B. It is generally consistent without averaging.

C. It can have substantial variance even for long records.

D. It is independent of the selected window.

Correct Answer: C. It can have substantial variance even for long records.

Question 2

For a record sampled at 1,000 Hz and divided into segments of 1,000 samples, the DFT frequency spacing is:

A. 0.1 Hz

B. 1 Hz

C. 10 Hz

D. 1,000 Hz

Correct Answer: B. 1 Hz

Question 3

What is the primary effect of zero-padding a finite record?

A. It increases the information content of the signal.

B. It eliminates spectral leakage.

C. It increases the density of plotted frequency samples.

D. It reduces sampling-frequency aliasing.

Correct Answer: C. It increases the density of plotted frequency samples.

Question 4

Welch's method primarily reduces spectral-estimate variance by:

A. Averaging modified periodograms from segments.

B. Increasing the signal's sampling rate.

C. Removing all window sidelobes.

D. Replacing the data with an AR model.

Correct Answer: A. Averaging modified periodograms from segments.

Question 5

In an AR model, the driving innovation e[n] is ordinarily assumed to be:

A. A deterministic sinusoid

B. White noise

C. A window function

D. A frequency response

Correct Answer: B. White noise

Question 6

Magnitude-squared coherence is bounded, under standard spectral definitions, by:

A. −1 and 1

B. 0 and 1

C. 0 and infinity

D. −infinity and infinity

Correct Answer: B. 0 and 1

Question 7

A shorter STFT window generally provides:

A. Better time localization and poorer frequency localization

B. Poorer time localization and better frequency localization

C. Perfect localization in both domains

D. No dependence on resolution

Correct Answer: A. Better time localization and poorer frequency localization

Question 8

Which method naturally provides shorter effective windows at high frequencies and longer windows at low frequencies?

A. Rectangular periodogram

B. Wavelet transform

C. Cross-spectrum

D. White-noise model

Correct Answer: B. Wavelet transform

Question 9

A major limitation of the Wigner–Ville distribution for multicomponent signals is:

A. It cannot represent frequency.

B. It always has zero time resolution.

C. It can produce cross-terms between components.

D. It requires no sampling.

Correct Answer: C. It can produce cross-terms between components.

Question 10

High coherence between two signals at a frequency establishes:

A. Definite physical causation

B. A stable linear association at that frequency

C. That both signals are identical in time

D. That aliasing is absent

Correct Answer: B. A stable linear association at that frequency