pwelch {gsignal}

Welch’s power spectral density estimate

Description

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

Usage

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

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.

Value

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

Note

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

Author(s)

Peter V. Lanspeary [e-mail].
Conversion to R by Geert van Boxtel, [e-mail].

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

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)