View Raw SPL
/*****************************************************************************
* *
* BANDPOWER.SPL Copyright (C) 2024 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Compute the power of a series within a frequency range *
* *
* Revisions: 5 Mar 2024 RRR Creation *
* *
*****************************************************************************/
#if @HELP_BANDPOWER
BANDPOWER
Purpose: Computes the mean squared power of a series within a frequency
range.
Syntax: BANDPOWER(s, fs, frange)
s - A series or array, the input time domain series.
fs - Optional. A scalar, the sample rate of the input
data. Defaults to RATE(x).
frange - Optional. A two element array of frequency values,
the frequency interval to compute the series power.
Defaults to the entire frequency range.
Alternate Syntax: BANDPOWER(Pxx, f, frange, "psd")
Pxx - A series or array, the input frequency domain
power spectrum.
f - Optional. A series, the input frequency values.
Defaults to XVALS(Pxx).
frange - Optional. A two element array of frequency values,
the frequency interval to compute the series power.
"psd" - A string. The flag "PSD" used to indicate the input
is a frequency domain PSD. If the flag is "POWSPEC",
the input is a power spectrum.
Returns: A scalar or array, the series power over the frequency range.
Example:
W1: 1..100
W2: {bandpower(w1)}
W3: {sum(powspec(w1))}
W4: {mean(w1 * w1)}
W2 == 3383.5
W3 == 3383.5
W4 == 3383.5
W2, W3 and W4 compute the mean squared power of the series in
W1.
W2 and W3 perform identical calculations in the frequency
domain.
W4 performs the mean squared time domain calculation.
This form of BANDPOWER computes the sum of the power spectrum
computed by POWSPEC. From Parseval's theorem, the sum of the
power spectrum terms equals the mean of the series squared,
i.e.:
sum(powspec(x)) == mean(x*x)
Example:
W1: 1..100
W2: {bandpower(w1, {0.0, 0.5})}
W3: {sum(powspec(hamming(w1, 3)))}
W2 == 2779.693386
W3 == 2779.693386
Computes the series power over the entire frequency range. In
this case, BANDPOWER multiplies the input series by a HAMMING
window corrected for the mean squared amplitude.
Example:
W1: 1..100
W2: {bandpower(w1, rate(w1), {0.0, 0.5})}
W3: {sum(powspec(hamming(w1, 3)))}
W2 == 2779.693386
W3 == 2779.693386
Same as above except the sample rate is explicit.
Example:
W1: 1..100
W2: {bandpower(w1, {0.1, 0.3})}
W2 == 0.624722
Computes the series power of W1 over the frequency range
0.1 Hz to 0.3 Hz.
Example:
W1: gsin(1000, 1/1000, 20)^2
W2: gsqrwave(10000, 1/10000, 100)
W3: {bandpower(W1, {30, 50})}
W4: {bandpower(W2, {1000, 2000})}
W5: ravel(w1, w2)
W6: {bandpower(W5, ravel({30, 50}, {1000, 2000}))}
W3 == 0.125000
W4 == 0.005386
W6 == {{0.12500, 0.005386}}
W1 and W2 contain two different series with different sample
rates.
W3 computes the series power of W1 between 30 Hz and 50 Hz.
W4 computes the series power of W2 between 1000 Hz and 2000 Hz.
W5 combines the data in W1 and W2 into a single two column
series.
W6 computes the series power of each column of W5, producing the
same result as the individual computations in W3 and W4.
Example:
W1: seedrand(1);gsin(1000, 0.001, 100)^2+gnorm(1000, 0.001)*0.2
W2: psd(hamming(w1, 3))
W3: {bandpower(w1, {150, 300})}
W4: {bandpower(w2, {150, 300}, "psd")}
W3 == 0.144707
W4 == 0.144707
W1 contains a squared 100 Hz sinewave sampled at 10 kHz with
additive normally distributed random noise. The random number
generated is set with a seed of 1 for consistent repeated
results.
W2 computes the frequency domain power spectrum of W1. The
data is multiplied by a Hamming window with mean squared
amplitude correction.
W3 computes the series power of W1 between 150 Hz and 300 Hz.
W4 computes the series power of W1 by summing the power spectrum
in W2 between 150 Hz and 300 Hz.
W3 and W4 produce the same result.
Example:
W1: seedrand(1);gsin(1000, 0.001, 100)^2+gnorm(1000, 0.001)*0.2
W2: powspec(hamming(w1, 3))
W3: {bandpower(w1, rate(w1), {150, 300})}
W4: {bandpower(w2, xvals(w2), {150, 300}, "powspec")}
W3 == 0.144707
W4 == 0.144707
Same as above except the sample rate and frequencies are
explicit and the input to W4 is the power spectrum in W2.
Remarks:
BANDPOWER computes the series power by summing the appropriate
bands of the PSD as calculated by PSD or the power spectrum
as calculated by POWSPEC.
If a frequency range is provided, the input series is multiplied
by a mean squared amplitude corrected Hamming window, otherwise
all frequencies of the un-windowed power spectrum are summed.
The input is considered a frequency domain PSD if PSDFLAG is
"psd" and a frequency domain power spectrum if PSDFLAG is
"powspec".
See Also:
Hamming
Meanfreq
Medfreq
Powspec
Psd
#endif
/* compute mean squared series power in the frequency domain */
ITERATE bandpower(s, f, frange, psdflag)
{
local ps, P;
/* parse and process input args to get the power spectrum */
ps = freqband_process_args(__FUNC__, s, f, frange, psdflag);
/* total power */
P = sum(ps);
return(P);
}