Understanding Power Spectral Density (PSD) for Acceleration Data

Power Spectral Density (PSD) is a key tool for analyzing vibration data, particularly in structural health monitoring, earthquake engineering, and mechanical systems. This blog explores what PSD is, why it matters, and how to compute it from acceleration time series data using practical Python examples.

What is Power Spectral Density (PSD)?

PSD represents how the power of a signal is distributed across different frequencies. When analyzing acceleration signals, PSD helps identify dominant frequencies, structural resonances, and noise characteristics.

Mathematically, PSD can be estimated using the Discrete Fourier Transform (DFT):

$S_{aa} (f_k) = \frac{{|A[k]|}^2}{Nf_s}$

where $f_k$ are the frequency bins, $N$ is the number of samples, and $f_s$ is the sampling frequency.

Steps to Compute PSD from Acceleration Time Series

To compute the PSD of an acceleration time series, follow these steps:

1. Obtain the Acceleration Time Series

Collect data from an accelerometer at a known sampling frequency $f_s$ .
Ensure the data is recorded over a sufficiently long duration to capture relevant frequencies.

2. Preprocess the Data

Demean the Signal: Remove the mean to eliminate DC bias.
$a'[n] = a[n] - \frac{1}{N} \sum_{i=0}^{N-1}a[i]$
Detrending: Remove any linear trends that may affect spectral estimation.
Apply a Window Function: Use a window function (e.g., Hann, Hamming) to reduce spectral leakage.

3. Compute the Fourier Transform

Use the Fast Fourier Transform (FFT) to convert the time-domain signal to the frequency domain.
$A[k] = \sum_{n=0}^{N-1} a'[n] e^{-j2\pi kn/N}$
The result $A[k]$ is the frequency-domain representation of the signal.

4. Estimate the PSD

Compute the squared magnitude of the FFT and normalize it by the signal length and sampling frequency.
$PSD (f_k) = \frac{|A[k]|^2}{Nf_s}$
Alternatively, use Welch’s method for improved spectral estimation by dividing the signal into multiple overlapping segments, applying a window function to each segment to minimize spectral leakage, computing individual periodograms, and averaging them to obtain a more stable and less noisy estimate of the PSD.

5. Plot and Analyze the PSD

Convert PSD values to decibels (dB) if needed.
Identify dominant frequencies that provide insights into system behavior.

Methods for Estimating PSD

Different techniques can be used to estimate PSD:

1. Periodogram

Direct estimation using the FFT.
High variance, making it less reliable.

2. Welch’s Method (Recommended)

Splits the signal into overlapping segments.
Applies a window function to reduce spectral leakage.
Averages periodograms for better accuracy.

3. Multitaper Method

Uses multiple window functions.
Further reduces spectral leakage and variance.

Computing PSD in Python

Let’s compute PSD using Welch’s method in Python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch

# Generate synthetic acceleration data
fs = 100  # Sampling frequency (Hz)
T = 10    # Duration (seconds)
t = np.linspace(0, T, T * fs, endpoint=False)
acceleration = np.sin(2 * np.pi * 5 * t) + 0.5 * np.random.randn(len(t))  # 5 Hz signal + noise

# Compute PSD
frequencies, psd_values = welch(acceleration, fs=fs, nperseg=1024, scaling='density')

# Plot PSD
plt.figure(figsize=(8, 4))
plt.semilogy(frequencies, psd_values)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD ($m^2/s^3$/Hz)')
plt.title('Power Spectral Density of Acceleration')
plt.grid()
plt.show()

Why PSD Matters

1. Structural Health Monitoring

Identifies natural frequencies of buildings and bridges.
Detects structural damage over time.

2. Earthquake Engineering

Analyzes ground motion frequency content.
Supports seismic hazard assessment.

3. Machinery Vibration Analysis

Detects faults in rotating machinery.
Identifies wear and imbalance.