PWELCH

Purpose:

Estimates the PSD using the Welch method of segment averaging.

Syntax:

PWELCH(series, win, overlap, nfft, fs, detrend, zeropad, range)

Returns:

A series or array, the PSD estimate of the input.

Example:

W1: 1..1024

W2: pwelch(W1, 128)

W3: pwelch(W1, hamming(128), 64, 128)

W2 and W3 produce the same PSD estimate by dividing the series into segments of length 128 that overlap each other by 64 points. Each segment is multiplied by a Hamming window and a 128 point PSD is computed. The PSD's are averaged to produce the final PSD estimate.

Example:

W1: gsin(10000, 0.0001, 2000);W0 + gline(length, deltax, 100, 0);setvunits("V")

W2: pwelch(W1, hanning(128, "periodic"));semilogy

W3: pwelch(W1, hanning(128, "periodic"), 64, 1024);semilogy

W4: pwelch(W1, hanning(128, "periodic"), 64, 1024, "linear");overp(w3);semilogy

W2, W3 and W4 compute the PSD using the Welch method with a periodic Hanning window. The PSD results are displayed with log y scales.

W4 removes the linear trend from each segment before FFT processing. The result from W3 is overplotted for a visual comparison. Removing the linear trend accentuates the peak.

Example:

W1: 1..1024;setdeltax(1/length)

W2: pwelch(W1, hamming(128), 64, 128)

W3: pwelch(W1, hamming(128), 64, {0, 8, 16, 24})

W4: w2;overp(w3);semilogy(1)

W2 produces a standard PSD estimate by dividing the series into segments of length 128 that overlap each other by 64 points. The result is 65 points long where the frequencies range from 0 to 512 Hz with a spacing of 8 Hz.

W3 computes the same PSD estimate except the output frequencies are explicitly set to 0, 8, 16 and 24 Hz, i.e. the first 4 values of the PSD in W2.

W4 compares the 65 point PSD of W2 to the 4 point PSD of W4 on a log Y basis. The values in W3 are identical to the values in W2 within machine precision.

Example:

A := 10.0;

N := 100000;

f := 2000.13;

W1: A * gcos(N, 1/N, f);setvunits("V")

W2: psd(w1)

W3: pwelch(W1, rectwin(N), 0, N)

W4: pwelch(W1, rectwin(N), 0, {f})

W5: pwelch(W1, rectwin(N), 0, N*50)

W6: xy({maxloc(w3)}, {max(w3)});

W7: xy({maxloc(w4)}, {max(w4)});

W8: xy({maxloc(w5)}, {max(w5)});

W9: {(deltax(w1) * length(w1) * A^2) / 2}

W1 contains a 100000 point cosine with a sample rate of 100 kHz and a frequency of 2000.13 Hz.

W2 computes a PSD estimate using all the data in W1 with no windowing and no segmenting.

W3 produces a PSD estimate using a single segment that is the length of the input with no overlap. The result is identical to the PSD function in W2.

W4 computes a PSD estimate at a single target frequency f.

W5 computes a PSD estimate using the same form as W3 except the input data is zero padded to 50 times the original length.

W6 displays the frequency and magnitude of the detected cosine using the single periodogram.

W7 displays the frequency and magnitude of the detected cosine when the target frequency, f, is explicit.

W8 displays the frequency and magnitude of the detected cosine using the single periodogram with a 50x zero padded input series.

W9 displays the value 50.0, the analytic PSD maximum of the input. The magnitude at f is ideally T L A² / 2 where T = deltax(w1) and L = length(w1).

The single periodogram approaches of W2 and W3 produce an inaccurate result for the cosine magnitude because the input frequency does not lie at the center of a single FFT frequency bin.

The PWELCH method of W3 with an explicit target frequency produces a value that compares well with the ideal result. The single periodogram method of W4 with the zero padded extended input produces the same result as W3 but with a much longer running time and memory requirements.

W3 computes the same PSD estimate except the output frequencies are explicitly set to 0, 8, 16 and 24 Hz, i.e. the first 4 values of the PSD in W2.

W4 compares the 65 point PSD of W2 to the 4 point PSD of W4 on a log Y basis. The values in W3 are identical to the values in W2 within machine precision.

Remarks:

Welch's method of Power Spectral Density (PSD) estimation is a technique that reduces the variance of the basic periodogram by averaging multiple modified periodograms.

1. The input series x[n] of length N is divided into K overlapping segments. Each segment has length L, normally a power of 2 for efficient FFT evaluation. The overlap is typically 50% (L/2 samples overlap each segment). The i-th segment is:

    

   

    

   where:

    

   D = L − S is the hop size,

    

   S is the overlap length (e.g., S = L/2 for 50% overlap),

    

   i = 0, 1, …, K−1.

    

   and K is the number of segments:

    

   

    

2. Each segment is multiplied by a tapering window w[n] (e.g., HAMMING, HANNING) to reduce spectral leakage:

    

   

    

3. For each windowed segment, a modified periodogram, P_i( f ), is calculated by taking the Discrete Fourier Transform (DFT) of the windowed segment and computing the squared magnitude.

    

   

    

   where F_s is the sampling frequency and U is a window normalization factor to ensure the total power is preserved:

    

   

    

4. The final PSD estimate is the average of all the individual periodograms:

    

   

    

   This averaging step significantly reduces the variance of the estimate compared to a single, full-length periodogram.

    

PWELCH estimates the PSD of the input series by following the steps listed above: the input series is divided into equal sized, overlapping segments, each segment is multiplied by win and the normalized average of the periodograms of all the windowed segments is returned.

If nfft is an array of one or more target frequencies, PWELCH uses the GOERTZEL algorithm to compute the PSD magnitudes at each input frequency. The input frequencies can be any real values. Frequency values less than zero or greater than rate(series) are wrapped to the range between 0 and rate(series).

If detrend is "constant" or "linear", the mean or linear trend is computed and removed from each segment before FFT processing.

If zeropad is "zeropad", the last column is padded with zeros to make it the same length as the other columns. If zeropad is "none" or empty, the last column is ignored if it is shorter than the others.

If range is "twosided", the full PSD is displayed from 0 to rate(series) in wrap-around order. For example, a 10 point full PSD with values at frequencies:

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

is mapped to frequencies:

{-4, -3, -2, -1, 0, 1, 2, 3, 4, 5}

If range is "center", the full PSD is displayed from -rate(series)/2 to rate(series)/2 with the zero frequency in the middle. For example, a 10 point full PSD with values at frequencies:

{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}

is mapped to frequencies:

{0, 1, 2, 3, 4, 5, -4, -3, -2, -1}

If both inputs are real, range defaults to "onesided", else if any input is complex, range defaults to "twosided".

The PSD as estimated by PWELCH is also referred to as the auto-power spectral density.

See PSD to compute a single periodogram.

See CPSD to estimate the complex cross-power spectral density.

See MSCOHERE to estimate the magnitude-squared coherence function.

See TFESTIMATE to estimate the complex transfer function.

See DADiSP/FFTXL to optimize the underlying FFT computation.

See Also:

CPSD

FFT

GOERTZELDFT

HAMMING

HANNING

INVPSD

KAISER

MSCOHERE

PSD

SPECTRUM

TFESTIMATE

References: