R语言spec.pgram/spectrum功率谱缩放规则及信号轴正确缩放咨询

Let’s break down exactly how to set up the axes correctly when using spec.pgram() or spectrum() in R, given your signal parameters:

• Signal duration: 100ms (0.1s)
• Sampling rate (fs): 1000Hz (so sampling interval Δt = 0.001s)
• Total samples: N = 0.1s * 1000Hz = 100 samples

Key Background on R’s Scaling Rule #

R’s spectral functions follow the S-PLUS definition, using a scaling factor of 1/frequency(x) to ensure the power spectral density (PSD) behaves as a density over the frequency range (-fs/2, fs/2]. Here, frequency(x) refers to your signal’s sampling rate (1000Hz), so the scaling factor becomes 1/1000.

1. Frequency Axis (X-Axis) Setup #

R’s spectrum()/spec.pgram() will handle this automatically if you specify the sampling rate correctly in the function call:

# Assume your signal is stored in variable 'x'
spec <- spectrum(x, frequency = 1000, plot = FALSE)

• The spec$freq vector will contain actual frequency values in Hz, ranging from 0 up to the Nyquist frequency (500Hz, since fs/2 = 1000/2).
• If you forget to set frequency = 1000, the function will default to normalized frequencies (0 to 0.5). To convert these to Hz, multiply each value by 1000.
• When you plot the spectrum with plot(spec), the X-axis will automatically label itself in Hz (thanks to the frequency parameter).

2. Power Spectral Density Axis (Y-Axis) Setup #

The scaling factor 1/fs ensures the PSD has units of (signal units)²/Hz (e.g., V²/Hz if your signal is voltage). This scaling guarantees that:

• The integral of the PSD over the full frequency range (-500Hz to 500Hz) equals the variance of your original signal.
• Since R defaults to a one-sided spectrum (only positive frequencies), you can verify the scaling by checking that sum(spec$spec * diff(spec$freq)[1]) * 2 is approximately equal to var(x) (the factor of 2 accounts for the mirrored negative frequency energy).

Quick Example Workflow #

# Generate a test signal matching your parameters
set.seed(123)
fs <- 1000
duration <- 0.1
n_samples <- fs * duration
x <- rnorm(n_samples)  # White noise signal

# Compute spectrum with correct scaling
spec <- spectrum(x, frequency = fs, plot = TRUE)

# Verify variance match
psd_integral <- sum(spec$spec * diff(spec$freq)[1]) * 2
cat("Signal variance:", var(x), "\n")
cat("PSD integral (2x one-sided):", psd_integral, "\n")

You’ll see the values are very close, confirming the scaling is correct.

内容的提问来源于stack exchange，提问作者NicolasBourbaki