View Raw SPL
/*****************************************************************************
* *
* MEANFREQ.SPL Copyright (C) 2024 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Compute the mean frequency of a power spectrum *
* *
* Revisions: 26 Oct 2024 RRR Creation *
* *
*****************************************************************************/
#if @HELP_MEANFREQ
MEANFREQ
Purpose: Computes the mean frequency of the power spectrum of a series.
Syntax: MEANFREQ(x, fs, frange)
(f, p) = MEANFREQ(x, fs, frange)
x - 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: MEANFREQ(Pxx, f, frange, "psd")
(f, p) = MEANFREQ(Pxx, f, frange, "psd")
Pxx - A series or array, the input frequency domain
PSD, power spectral density.
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" - Optional. A string. The flag "PSD" can be used to
indicate the input is a frequency domain PSD. If the
flag is "POWSPEC", the input is a frequency domain
power spectrum as computed by POWSPEC.
Returns: A scalar or array, the mean frequency of the power spectrum
over the frequency range.
(f, p) = MEANFREQ(x, fs, frange) returns the mean frequency
and power over the frequency range.
Example:
W1: gsin(1000, 1/1000, 40) + gsin(1000, 1/1000, 160)
W2: {meanfreq(w1)}
W1 contains a 1000 samples series that is the sum of a 40 Hz
sinewave and a 160 Hz sinewave sampled at 1000 Hz.
W2 computes the mean frequency and converts the scalar result
to a one point series. The mean frequency is 100 Hz.
Example:
W1: gsin(1000, 1/1000, 40) + gsin(1000, 1/1000, 160)
W2: (f, p) = meanfreq(w1);{f, p};table
W2 == {100, 1}
Same as above, except W2 contains both the mean frequency,
100 Hz and the power, 1.0.
Example:
W1: gsin(1000, 1/1000, 40) + gsin(1000, 1/1000, 160);setvunits("V")
W2: {meanfreq(w1)}
W3: 10*log10(psd(w1));linedraw(w0, red, w2[1], -inf, w2[1], inf)
Same as the first example except W3 computes the PSD of the series
in W1 and a vertical line is drawn at the mean frequency.
Example:
W1: gsin(10000, 1/1000, 100)^3
W2: {meanfreq(w1)}
W3: psd(w1)
W4: {meanfreq(w3, "psd")}
W2 == 120
W4 == 120
W2 computes the mean frequency from the time domain series in W1.
W4 computes the mean frequency from the PSD in W3. The value
is identical to W2.
Example:
W1: gsin(10000, 1/1000, 100)^3
W2: {meanfreq(w1)}
W3: powspec(w1)
W4: {meanfreq(w3, "powspec")}
W2 == 120
W4 == 120
Same as above except W4 computes the mean frequency from the
power spectrum in W3.
Example:
W1: gsin(1000, 1/1000, 20)^2
W2: gsin(10000, 1/10000, 100)^2
W3: {meanfreq(W1, {30, 50})}
W4: {meanfreq(W2, {100, 300})}
W5: ravel(w1, w2)
W6: {meanfreq(W5, ravel({30, 50}, {100, 300}))}
W3 == 40
W4 == 200
W6 == {{40, 200}}
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(w1)
W3: {meanfreq(w1, {150, 300})}
W4: {meanfreq(w2, {150, 300}, "psd")}
W3 == 202.126447
W4 == 202.126447
W1 contains a squared 100 Hz sinewave sampled at 1 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 PSD of W1.
W3 computes the series power of W1 between 150 Hz and 300 Hz.
W4 computes the series power of W1 between 150 Hz and 300 Hz
directly from the PSD in W2.
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(w1)
W3: {meanfreq(w1, rate(w1), {150, 300})}
W4: {meanfreq(w2, xvals(w2), {150, 300}, "powspec")}
W3 == 202.126447
W4 == 202.126447
Same as above except the sample rate and frequencies are
explicit and the input to W4 is the power spectrum in W2.
Remarks:
MEANFREQ computes the average frequency of the energy
distribution of the series spectrum. It is the weighted
average of the frequencies in the spectrum where the weights
are proportional to the power at each frequency.
The input is considered a frequency domain power spectrum if
PSDFLAG is "psd". If PSDFLAG is "powspec, the input is considered
a power spectrum as computed by POWSPEC.
See MEDFREQ to compute the median frequency.
See Also:
Bandpower
Hamming
Medfreq
Powspec
Psd
#endif
ITERATE meanfreq(s, f, frange, psdflag)
{
local ps, P, mf;
/* parse and process input args to get the power spectrum */
ps = freqband_process_args(__FUNC__, s, f, frange, psdflag);
/* total power */
P = sum(ps);
/* weighted frequencies divided by total power */
mf = sum(ps * xvals(ps)) / P;
return(mf, P);
}