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, fs, frange, "psdflag")

(f, p) = MEDFREQ(Pxx, fs, frange, "psdflag")

Pxx

-

A series or array, the input frequency domain PSD (power spectral density) or power spectrum.

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.

"psdflag"

-

Optional. A string. The string "PSD" indicates the input is a frequency domain PSD. The string "POWSPEC"  indicates 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, range) 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 sine wave and a 160 Hz sine wave 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.

 

W3 computes the power spectral density of 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.005556

W4 == 100.005556

 

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 sine wave 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

MEANFREQ

POWSPEC

PSD

XYLOOKUP