pwelch() R function from [gsignal]

Welch’s power spectral density estimate

Compute power spectral density (PSD) using Welch's method.


pwelch(
  x,
  window = nextpow2(sqrt(NROW(x))),
  overlap = 0.5,
  nfft = if (isScalar(window)) window else length(window),
  fs = 1,
  detrend = c("long-mean", "short-mean", "long-linear", "short-linear", "none"),
  range = if (is.numeric(x)) "half" else "whole"
)

## S3 method for class 'pwelch'
plot(
  x,
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  plot.type = c("spectrum", "cross-spectrum", "phase", "coherence", "transfer"),
  yscale = c("linear", "log", "dB"),
  ...
)

## S3 method for class 'pwelch'
print(
  x,
  plot.type = c("spectrum", "cross-spectrum", "phase", "coherence", "transfer"),
  yscale = c("linear", "log", "dB"),
  xlab = NULL,
  ylab = NULL,
  main = NULL,
  ...
)

Arguments

x: input data, specified as a numeric vector or matrix. In case of a vector it represents a single signal; in case of a matrix each column is a signal.
window: If window is a vector, each segment has the same length as window and is multiplied by window before (optional) zero-padding and calculation of its periodogram. If window is a scalar, each segment has a length of window and a Hamming window is used. Default: nextpow2(sqrt(length(x))) (the square root of the length of x rounded up to the next power of two). The window length must be larger than 3.
overlap: segment overlap, specified as a numeric value expressed as a multiple of window or segment length. 0 <= overlap < 1. Default: 0.5.
nfft: Length of FFT, specified as an integer scalar. The default is the length of the window vector or has the same value as the scalar window argument. If nfft is larger than the segment length, (seg_len), the data segment is padded nfft - seg_len zeros. The default is no padding. Nfft values smaller than the length of the data segment (or window) are ignored. Note that the use of padding to increase the frequency resolution of the spectral estimate is controversial.
fs: sampling frequency (Hertz), specified as a positive scalar. Default: 1.
detrend: character string specifying detrending option; one of:
- long-mean: remove the mean from the data before splitting into segments (default)
- short-mean: remove the mean value of each segment
- long-linear: remove linear trend from the data before splitting into segments
- short-linear: remove linear trend from each segment
- none: no detrending
range: character string. one of:
- "half" or "onesided": frequency range of the spectrum is from zero up to but not including fs / 2. Power from negative frequencies is added to the positive side of the spectrum.
- "whole" or "twosided": frequency range of the spectrum is -fs / 2 to fs / 2, with negative frequencies stored in "wrap around order" after the positive frequencies; e.g. frequencies for a 10-point "twosided" spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2. -0.1.
- "shift" or "centerdc": same as "whole" but with the first half of the spectrum swapped with second half to put the zero-frequency value in the middle.
Default: If x are real, the default range is "half", otherwise the default range is "whole".
xlab, ylab, main: labels passed to plotting function. Default: NULL
plot.type: character string specifying which plot to produce; one of "spectrum", "cross-spectrum", "phase", "coherence", "transfer"
yscale: character string specifying scaling of Y-axis; one of "linear", "log", "dB"
...: additional arguments passed to functions

Returns

An object of class "pwelch", which is a list containing the following elements:

freq: vector of frequencies at which the spectral variables are estimated. If x is numeric, power from negative frequencies is added to the positive side of the spectrum, but not at zero or Nyquist (fs/2) frequencies. This keeps power equal in time and spectral domains. If x is complex, then the whole frequency range is returned.
spec: Vector (for univariate series) or matrix (for multivariate series) of estimates of the spectral density at frequencies corresponding to freq.
cross: NULL for univariate series. For multivariateseries, a matrix containing the cross-spectral density estimates between different series. Column $i + (j - 1) * (j - 2)/2$ of contains the cross-spectral estimates between columns $i$ and $j$ of $x$ , where $i < j$ .
phase: NULL for univariate series. For multivariate series, a matrix containing the cross-spectrum phase between different series. The format is the same as cross.
coh: NULL for univariate series. For multivariate series, a matrix containing the squared coherence between different series. The format is the same as cross.
trans: NULL for univariate series. For multivariate series, a matrix containing estimates of the transfer function between different series. The format is the same as cross.
x_len: The length of the input series.
seg_len: The length of each segment making up the averages.
psd_len: The number of frequencies. See freq
nseries: The number of series
series: The name of the series
snames: For multivariate input, the names of the individual series
window: The window used to compute the modified periodogram
fs: The sampling frequency
detrend: Character string specifying detrending option

Details

The Welch method [1] reduces the variance of the periodogram estimate to the PSD by splitting the signal into (usually) overlapping segments and windowing each segment, for instance by a Hamming window. The periodogram is then computed for each segment, and the squared magnitude is computed, which is then averaged for all segments. See also [2].

The spectral density is the mean of the modified periodograms, scaled so that area under the spectrum is the same as the mean square of the data. This equivalence is supposed to be exact, but in practice there is a mismatch of up to 0.5

the data.

In case of multivariate signals, Cross-spectral density, phase, and coherence are also returned. The input data can be demeaned or detrended, overall or for each segment separately.

Note

Unlike the 'Octave' function 'pwelch', the current implementation does not compute confidence intervals because they can be inaccurate in case of overlapping segments.

Examples


fs <- 256
secs <- 10
freq <- 30
ampl <- 1
t <- seq(0, secs, length.out = fs * secs)

x <- ampl * cos(freq * 2 * pi * t) + runif(length(t))
Pxx <- pwelch(x, fs = fs)              # no plot
pwelch(x, fs = fs)                     # plot

# 90 degrees phase shift with with respect to x
y <- ampl * sin(freq * 2 * pi * t) + runif(length(t))
Pxy <- pwelch(cbind(x, y), fs = fs)
plot(Pxy, yscale = "dB")
plot(Pxy, plot.type = "phase")
# note the phase shift around 30 Hz is pi/2
plot(Pxy, plot.type = "coherence")

# Transfer function estimate example
fs <- 1000                 # Sampling frequency
t <- (0:fs) / fs           # One second worth of samples
A <- c(1, 2)               # Sinusoid amplitudes
f <- c(150, 140)           # Sinusoid frequencies
xn <- A[1] * sin(2 * pi * f[1] * t) +
      A[2] * sin(2 * pi * f[2] * t) +  0.1 * runif(length(t))
h <- Ma(rep(1L, 10) / 10)      # Moving average filter
yn <- filter(h, xn)
atfm <- freqz(h, fs = fs)
etfm <- pwelch(cbind(xn, yn), fs = fs)
op <- par(mfrow = c(2, 1))
xl <- "Frequency (Hz)"; yl <- "Magnitude"
plot(atfm$w, abs(atfm$h), type = "l", main = "Actual", xlab = xl, ylab = yl)
plot(etfm$freq, abs(etfm$trans), type = "l", main = "Estimated",
     xlab = xl, ylab = yl)
par(op)

References

[1] Welch, P.D. (1967). The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, AU-15 (2): 70–73.

[2] https://en.wikipedia.org/wiki/Welch%27s_method