View Raw SPL
/*****************************************************************************
* *
* MEDFREQ.SPL Copyright (C) 2024 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Randy Race *
* *
* Synopsis: Compute the median frequency of a power spectrum *
* *
* Revisions: 26 Oct 2024 RRR Creation *
* *
*****************************************************************************/
#if @HELP_MEDFREQ
MEDFREQ
Purpose: Computes the median frequency of the power spectrum of a series.
Syntax: MEDFREQ(x, fs, frange)
(f, p) = MEDFREQ(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: MEDFREQ(Pxx, f, frange, "psd")
(f, p) = MEDFREQ(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 median frequency of the power spectrum
over the frequency range.
(f, p) = MEDFREQ(x, fs, frange) returns the median frequency
and power over the frequency range.
Example:
W1: gsin(1000, 1/1000, 40) + gsin(1000, 1/1000, 160)
W2: {medfreq(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 median frequency and converts the scalar
result to a one point series. The median frequency is 40.5 Hz.
Example:
W1: gsin(1000, 1/1000, 40) + gsin(1000, 1/1000, 160)
W2: (f, p) = medfreq(w1);{f, p};table
W2 == {40.5, 1}
Same as above, except W2 contains both the median frequency,
40.5 Hz and the power, 1.0.
Example:
W1: gsin(1000, 1/1000, 40) + gsin(1000, 1/1000, 160);setvunits("V")
W2: {medfreq(w1)}
W3: 10*log10(powspec(w1));linedraw(w0, red, w2[1], -inf, w2[1], inf)
Same as the first example except W3 computes the power spectrum
of the series in W1 and a vertical line is drawn at the median
frequency.
Example:
W1: gsin(10000, 1/1000, 100)^3
W2: {medfreq(w1)}
W3: psd(w1)
W4: {medfreq(w3, "psd")}
W2 == 100.005556
W4 == 100.005556
W2 computes the median frequency from the time domain series
in W1.
W4 computes the median frequency from the PSD in W3. The value
is identical to W2.
Example:
W1: gsin(10000, 1/1000, 100)^3
W2: {medfreq(w1)}
W3: powspec(w1)
W4: {medfreq(w3, "powspec")}
W2 == 100.00556
W4 == 100.00556
Same as above except W4 computes the median frequency from the
power spectrum in W3.
Example:
W1: gsin(1000, 1/1000, 20)^2
W2: gsin(10000, 1/10000, 100)^2
W3: {medfreq(W1, {30, 50})}
W4: {medfreq(W2, {100, 300})}
W5: ravel(w1, w2)
W6: {medfreq(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: {medfreq(w1, {150, 300})}
W4: {medfreq(w2, {150, 300}, "psd")}
W3 == 200.014561
W4 == 200.014561
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: {medfreq(w1, rate(w1), {150, 300})}
W4: {medfreq(w2, xvals(w2), {150, 300}, "powspec")}
W3 == 200.014561
W4 == 200.014561
Same as above except the sample rate and frequencies are
explicit and the input to W4 is the power spectrum in W2.
Remarks:
MEDFREQ computes the median frequency of the energy
distribution of the series spectrum. The median frequency
is the frequency that divides the spectrum into two halves
of equal power.
The returned frequency is linearly interpolated within the
frequency band that is determined to contain the median
frequency.
The input is considered a frequency domain PSD if PSDFLAG
is "psd". If PSDFLAG is "powspec, the input is considered
a power spectrum as computed by POWSPEC.
See MEANFREQ to compute the mean frequency.
See Also:
Bandpower
Hamming
Meanfreq
Powspec
Psd
Xylookup
#endif
/* find the median frequency of the power spectrum */
ITERATE medfreq(s, f, frange, psdflag)
{
local ps, cp, minx, maxx, cf, P1, P2, hp, P, mf;
/* parse and process input args to get the power spectrum */
(ps, frange) = freqband_process_args(__FUNC__, s, f, frange, psdflag);
/* cumulative power */
cp = cumsum(ps);
/* frequency range */
(minx, maxx) = xminmax(cp);
/* bin width shift for narrowband signals */
cp = delay(cp, 1);
/* shift frequencies by 1/2 binwidth */
cf = xvals(cp) - deltax(ps) / 2;
cf[1] = 0;
cf[end] = maxx;
if (length(frange) > 1)
{
/* lookup and linearly interpolate half power */
P1 = xylookup(cf, cp, frange[1]);
P2 = xylookup(cf, cp, frange[2]);
/* half power threshold within the band */
hp = (P1 + P2) / 2;
/* full power in band */
P = P2 - P1;
}
else
{
/* full power */
P = max(cp);
/* half power = total power / 2 */
hp = P / 2;
}
/* lookup and linearly interpolate half power frequency */
mf = xylookup(cp, cf, hp);
return(mf, P);
}