Power spectral density estimation in Signal Processing

Introduction

Power spectral density (PSD) estimation helps us understand how the power of a signal is spread across different frequencies. It shows which frequencies carry the most energy.

Analyzing sound signals to find dominant tones or noise.

Studying heart rate signals to detect irregular rhythms.

Checking vibration signals in machines to find faults.

Examining brain waves to identify different brain states.

Understanding radio signals to improve communication quality.

Syntax

Signal Processing

frequencies, psd_values = signal.welch(signal_data, fs=sampling_rate, nperseg=segment_length)

signal_data is your input signal array.

fs is the sampling rate (how many samples per second).

nperseg is the length of each segment used for estimation.

Examples

Estimate PSD of data sampled at 1000 Hz with default segment length.

Signal Processing

frequencies, psd = signal.welch(data, fs=1000)

Estimate PSD with 500 Hz sampling rate and segments of 256 samples.

Signal Processing

frequencies, psd = signal.welch(data, fs=500, nperseg=256)

Use a Hann window to reduce spectral leakage during PSD estimation.

Signal Processing

frequencies, psd = signal.welch(data, fs=2000, nperseg=512, window='hann')

Sample Program

This program creates a signal with two sine waves at 50 Hz and 120 Hz plus some noise. It then estimates the power spectral density using Welch's method and prints the first few frequency and PSD values. Finally, it plots the PSD on a logarithmic scale to show the power distribution across frequencies.

Signal Processing

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

# Create a sample signal: 2 sine waves + noise
fs = 1000  # Sampling frequency 1000 Hz
T = 1.0    # seconds
N = int(T * fs)  # Number of samples

t = np.linspace(0, T, N, endpoint=False)
sig = 0.5 * np.sin(2 * np.pi * 50 * t) + 0.3 * np.sin(2 * np.pi * 120 * t)
sig += 0.1 * np.random.normal(size=N)  # Add noise

# Estimate power spectral density
frequencies, psd = signal.welch(sig, fs=fs, nperseg=256)

# Print first 5 frequency and PSD values
for f, p in zip(frequencies[:5], psd[:5]):
    print(f"Frequency: {f:.3f} Hz, PSD: {p:.6f}")

# Plot PSD
plt.semilogy(frequencies, psd)
plt.title('Power Spectral Density Estimation')
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD (V**2/Hz)')
plt.grid(True)
plt.show()

Important Notes

Welch's method splits the signal into overlapping segments to get a smoother PSD estimate.

Choosing segment length (nperseg) affects frequency resolution and variance of the estimate.

Window functions like Hann reduce spectral leakage but change the PSD shape slightly.

Summary

Power spectral density shows how signal power spreads over frequencies.

Welch's method is a common way to estimate PSD by averaging segments.

Adjust sampling rate and segment length to get better frequency details.