View Raw SPL
/*****************************************************************************
* *
* PWELCH.SPL Copyright (C) 2005-2025 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Welch method of PSD estimation *
* *
* Revisions: 25 Aug 2005 RRR Creation *
* 22 Jun 2024 RRR freq, detrend, zeropad and range options *
* *
*****************************************************************************/
#include
#if @HELP_PWELCH
PWELCH
Purpose: Calculates the PSD using the Welch method of segment averaging.
Syntax: PWELCH(s, win, olag, nfft, fs, detrend, zeropad, range)
s - A series or array, the input.
win - An integer or series, the segment size or window.
If a series, each input segment is multiplied by
WIN with the segment size set to length(WIN).
Defaults to a Hamming window of where the length
is 1/8 the length of the input series.
olap - Optional. An integer, the number of points to
overlap. Defaults to length(WIN)/2.
nfft - Optional. An integer or series. If NFFT is an
integer, NFFT is the FFT length. If NFFT is a
series, the series contains the frequencies used
to evaluate the PSD. Defaults to the integer
bestpow2(WIN).
fs - Optional. A real, the input series sample rate.
Defaults to rate(S).
detrend - Optional. A string, the segment detrend mode:
"none" : no detrending (default)
"constant" : remove the mean of the segment
"linear" : remove the linear trend of the
segment
zeropad - Optional. A string, the short segment zero pad mode:
"none" : no zero padding, ignore the last column
if it is shorter than the other columns
(default)
"zeropad" : add zeros to the last column if it is
shorter than the other columns to force
the same length for all columns
range - Optional. A string, the frequency range of the PSD.
"onesided" : one sided PSD from 0 to Fs/2 (default)
"twosided" : full PSD from 0 to Fs
"center" : centered PSD from -Fs/2 to Fs/2
Returns: A series or array.
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
W1 contains a 2000 Hz sinewave with an additive ramp samples at 10 kHz.
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.
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 / 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.
Remarks:
PWELCH estimates the PSD of the input series by first dividing
the series into equal sized, overlapping segments. Each segment
is multiplied by WIN and the average of the PSDs 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(S)
are wrapped to the range between 0 and rate(S).
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(S) 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:
{0, 1, 2, 3, 4, 5, -4, -3, -2, -1}
If RANGE is "center", the full PSD is displayed from
-rate(S)/2 to rate(S)/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:
{-4, -3, -2, -1, 0, 1, 2, 3, 4, 5}
The true PSD for the series is continuous. Because the
PSD function returns samples of the continuous power
spectral density, increasing the number of samples results
in the area providing a closer approximation of the mean
squared power.
See Also:
Cpsd
Hamming
Hanning
Kaiser
Mscohere
PSD
Tfestimate
#endif
/* psd using built-in welchseg function */
SERIES ITERATE pwelch(s, win={}, olap={}, nfft={}, fs={}, detrend="", zeropad="", range="")
{
local wpsd, seg, lseg, numsegs, dflag, welch_flag, f, hu, output = "mag";
/* parse input args */
if (argc < 1 || not(isarray(s)))
{
error(sprintf("%s - input series required", __FUNC__));
}
(win, olap, nfft, fs, dflag, range, output) = welch_parse_args(s, {}, win, olap, nfft, fs, detrend, zeropad, range, output, "", __FUNC__);
if (range != "onesided")
{
dflag += WELCH_F_FULL;
}
/* detrend option returned as dflag */
welch_flag = WELCH_F_PSD + dflag;
if (isarray(nfft))
{
/* explicit frequencies, use Goertzel algorithm */
welch_flag += WELCH_F_GOERTZEL;
if (min(nfft) < 0 || max(nfft) > rate(s) / 2)
{
range = "twosided";
welch_flag |= WELCH_F_FULL;
}
/* compute psd of each segment and average */
wpsd = welchseg(s, win, nfft, olap, length(nfft), fs, welch_flag);
if (range == "onesided")
{
/* onesided - check DC and Nyquist */
f = find(nfft == 0);
wpsd[f] /= 2;
f = find(nfft == rate(s) / 2);
wpsd[f] /= 2;
}
/* x units */
hu = gethunits(wpsd);
/* XY */
wpsd = xy(nfft, wpsd);
sethunits(wpsd, hu);
}
else
{
/* Normal Welch method usng FFT */
wpsd = welchseg(s, win, olap, nfft, fs, welch_flag);
}
/* normalize entire psd */
wpsd /= mean(win * win);
if (range == "centered")
{
wpsd = welchshift(wpsd);
}
if (outargc > 1)
{
return(yvals(wpsd), xvals(wpsd));
}
else
{
return(wpsd);
}
}