medfreq

R2026b

Median frequency

Syntax

fMed = medfreq(x)

fMed = medfreq(x,Fs)

fMed = medfreq(x,Ts)

fMed = medfreq(x,t,Domain="time")

fMed = medfreq(pxx,f)

fMed = medfreq(sxx,f,rbw)

fMed = medfreq(___,freqRange)

fMed = medfreq(___,Domain=sigDomain)

[fMed,bandPow] = medfreq(___)

medfreq(___)

Description

Time-Domain Inputs

fMed = medfreq(x) estimates the median normalized frequency, fMed, of the power spectrum of a time-domain signal, x.

• Specify x as a vector, matrix, or a MATLAB^® timetable (since R2026b). If x is a matrix or a timetable, then the function computes the median frequency of each column or variable independently.

• To compute the power spectrum, medfreq uses the periodogram function with a rectangular window and number of DFT points equal to the length of x.

fMed = medfreq(x,Fs) estimates the median frequency of x (vector or matrix) with a sample rate Fs.

fMed = medfreq(x,Ts) estimates the median frequency of x (vector or matrix) with a sample time Ts. (since R2026b)

fMed = medfreq(x,t,Domain="time") estimates the median frequency of x (vector or matrix) with time instances t. (since R2026b)

Note

Specify Domain="time" to distinguish time-domain inputs (x,t) from frequency-domain inputs (pxx,f).

Frequency-Domain Inputs

fMed = medfreq(pxx,f) returns the median frequency of a power spectral density (PSD) estimate pxx. The frequencies f correspond to the estimates in pxx.

fMed = medfreq(sxx,f,rbw) returns the median frequency of a power spectrum estimate sxx with resolution bandwidth rbw. The frequencies, f, correspond to the estimates in sxx.

Additional Options

fMed = medfreq(___,freqRange) specifies the frequency interval freqRange over which to compute the median frequency.

• This syntax can include any combination of input arguments from previous syntaxes, as long as you specify freqRange before the Domain name-value argument.

• If you specify the second input as empty and x is not a timetable, then the function assumes normalized frequencies.

fMed = medfreq(___,Domain=sigDomain) specifies the signal domain of the input. (since R2026b)

[fMed,bandPow] = medfreq(___) also returns the band power of the spectrum, bandPow. If you specify freqRange, then bandPow contains the band power within freqRange.

medfreq(___) with no output arguments plots the PSD or power spectrum and annotates the median frequency.

Examples

Median Frequency of Chirps

Generate 1024 samples of a chirp sampled at 1024 kHz. Specify the chirp so that it has an initial frequency of 50 kHz and reaches 100 kHz at the end of the sampling. Add white Gaussian noise such that the signal-to-noise ratio is 40 dB.

nSamp = 1024;
Fs = 1024e3;
SNR = 40;

t = (0:nSamp-1)'/Fs;
x0 = chirp(t,50e3,nSamp/Fs,100e3);

x1 = x0 + randn(size(x0))*std(x0)/db2mag(SNR);

Estimate the median frequency of the chirp. Plot the power spectral density (PSD) and annotate the median frequency.

medfreq(x1,Fs);

Generate another chirp. Specify an initial frequency of 200 kHz, a final frequency of 300 kHz, and an amplitude that is twice that of the first signal. Add white Gaussian noise.

x2 = 2*chirp(t,200e3,nSamp/Fs,300e3);
x2 = x2 + randn(size(x2))*std(x2)/db2mag(SNR);

Concatenate the chirps to produce a two-channel signal. Estimate the median frequency of each channel.

y = medfreq([x1 x2],Fs)

y = 1×2
105 ×

    0.7500    2.4999

Plot the PSDs of the two channels and annotate their median frequencies.

medfreq([x1 x2],Fs);

Add the two channels to form a new signal. Plot the PSD and annotate the median frequency.

medfreq(x1+x2,Fs);

Median Frequency of Sinusoids

Generate 1024 samples of a 100.123 kHz sinusoid sampled at 1024 kHz. Add white Gaussian noise such that the signal-to-noise ratio is 40 dB.

nSamp = 1024;
Fs = 1024e3;
SNR = 40;

t = (0:nSamp-1)'/Fs;
x0 = sin(2*pi*t*100.123e3);

x1 = x0 + randn(size(x0))*std(x0)/db2mag(SNR);

Use the periodogram function to compute the power spectral density (PSD) of the signal. Specify a Kaiser window with the same length as the signal and a shape factor of 38. Estimate the median frequency of the signal and annotate it on a plot of the PSD.

[Pxx,f] = periodogram(x1,kaiser(nSamp,38),[],Fs);

medfreq(Pxx,f);

Generate another sinusoid, this one with a frequency of 257.321 kHz and an amplitude that is twice that of the first sinusoid. Add white noise.

x2 = 2*sin(2*pi*t*257.321e3);
x2 = x2 + randn(size(x2))*std(x2)/db2mag(SNR);

Concatenate the sinusoids to produce a two-channel signal. Estimate the PSD of each channel and use the result to determine the median frequency.

[Pyy,f] = periodogram([x1 x2],kaiser(nSamp,38),[],Fs);

y = medfreq(Pyy,f)

y = 1×2
105 ×

    1.0012    2.5731

Annotate the median frequencies of the two channels on a plot of the PSDs.

medfreq(Pyy,f);

Add the two channels to form a new signal. Estimate the PSD and annotate the median frequency.

[Pzz,f] = periodogram(x1+x2,kaiser(nSamp,38),[],Fs);

medfreq(Pzz,f);

Median Frequency of Bandlimited Signals

Generate a signal whose PSD resembles the frequency response of an 88th-order bandpass FIR filter with normalized cutoff frequencies $$0.25\pi$$ rad/sample and $$0.45\pi$$ rad/sample.

d = fir1(88,[0.25 0.45]);

Compute the median frequency of the signal between $$0.3\pi$$ rad/sample and $$0.6\pi$$ rad/sample. Plot the PSD and annotate the median frequency and measurement interval.

medfreq(d,[],[0.3 0.6]*pi);

Output the median frequency and the band power of the measurement interval. Specifying a sample rate of $$2\pi$$ is equivalent to leaving the rate unset.

[mdf,power] = medfreq(d,2*pi,[0.3 0.6]*pi);

fprintf("Median = %.3f*pi, power = %.1f%% of total \n", ...
    mdf/pi,power/bandpower(d)*100)

Median = 0.371*pi, power = 77.4% of total 

Add a second channel with normalized cutoff frequencies $$0.5\pi$$ rad/sample and $$0.8\pi$$ rad/sample and an amplitude that is one-tenth that of the first channel.

d = [d;fir1(88,[0.5 0.8])/10]';

Compute the median frequency of the signal between $$0.3\pi$$ rad/sample and $$0.9\pi$$ rad/sample. Plot the PSD and annotate the median frequency of each channel and the measurement interval.

medfreq(d,[],[0.3 0.9]*pi);

Output the median frequency of each channel. Divide by $$\pi$$.

mdf = medfreq(d,[],[0.3 0.9]*pi)/pi

mdf = 1×2

    0.3706    0.6500

Median Frequency of Quadratic Chirp Timetable with Missing Samples

Since R2026b

Generate a timetable that comprises a three-second quadratic swept-frequency cosine signal with a sample rate of 1 kHz. The frequency sweeping ranges from 150 Hz to 350 Hz.

fs = 1000;
t = (1/fs:1/fs:3)';
f1 = 150;
f2 = 350;
sg = chirp(t,f1,t(end),f2,"quadratic") + 1e-3*randn(size(t));

sgTT = timetable(sg,SampleRate=fs);

Estimate the median frequency of the signal timetable. Plot the PSD and annotate the median frequency.

medfreq(sgTT);

You can estimate the median frequency of the input signal timetable. Remove one sample from the signal timetable every 0.1 seconds.

sgTT(mod(seconds(sgTT.Time),0.1)<1/(2*fs), :) = [];

Estimate the median frequency of the modified timetable. Plot the PSD and annotate the median frequency. The estimated median-frequency value approximates well with the median-frequency estimation from the unmodified timetable.

medfreq(sgTT);

Input Arguments

Time Domain

x — Input signal

vector | matrix | timetable

Input signal, specified as a real-valued vector, matrix, or timetable.

• If x is a vector, then the function treats x as a single-channel signal.

• If x is a matrix or timetable, then medfreq computes the median frequency of each column or variable independently.

• You can specify x as a nonuniformly sampled timetable. (since R2026b)

• x must be finite-valued.

Note

If you specify x as a nonuniformly sampled signal, then medfreq estimates an effective sample rate and performs uniform resampling before transforming the signal into the frequency domain. You can specify x as a nonuniformly sampled signal if the associated times tv satisfy these conditions:

• median(diff(tv))/mean(diff(tv)) < 100

• mean(diff(tv))/median(diff(tv)) < 100

The medfreq function obtains tv either from x.Time property if x is a timetable or from t if x is a vector or matrix.

Example: cos(pi/4*(0:159))+randn(1,160) is a single-channel row-vector signal.

Example: cos(pi./[4;2]*(0:159))'+randn(160,2) is a two-channel signal.

Data Types: single | double

Fs — Sample rate

positive real scalar

Sample rate, specified as [] (empty) or a positive scalar.

• [] — medfreq uses normalized frequencies in π rad/sample.

• Positive scalar — Use a sample rate, which is the number of samples per unit time. If the time is measured in seconds, then the sample rate is in hertz.

This argument does not apply if x is a timetable.

Data Types: single | double

Ts — Sample time

duration scalar

Since R2026b

Sample time, specified as a duration scalar.

This argument does not apply if x is a timetable.

Data Types: duration

t — Times associated with input signal

vector of nonnegative scalars | duration vector | datetime vector

Since R2026b

Times associated with the input signal, specified as a vector of nonnegative scalars, a duration vector, or a datetime vector. The function uses the value specified in this argument as the times associated with the vector or matrix x.

• If x is a vector, then x and t must have the same number of elements.

• If x is a matrix, then the number of rows in x and the number of elements in t must coincide.

• You can specify t with nonuniformly spaced elements. In this case, t represents the times associated with a nonuniformly sampled signal x.

This argument does not apply if x is a timetable

Data Types: single | double | duration | datetime

Frequency Domain

pxx — Power spectral density

vector | matrix

Power spectral density (PSD), specified as a vector or matrix. If pxx is a matrix, then medfreq computes the median frequency of each column of pxx independently.

The power spectral density must be expressed in linear units, not decibels. Use db2pow to convert decibel values to power values.

Example: [pxx,f] = periodogram(cos(pi./[4;2]*(0:159))'+randn(160,2)) specifies the periodogram PSD estimate of a noisy two-channel sinusoid sampled at 2π Hz and the frequencies at which it is computed.

Data Types: single | double

f — Frequencies

vector

Frequencies, specified as a vector.

Data Types: single | double

sxx — Power spectrum estimate

vector | matrix

Power spectrum estimate, specified as a vector or matrix. If sxx is a matrix, then medfreq computes the median frequency of each column of sxx independently.

The power spectrum must be expressed in linear units, not decibels. Use db2pow to convert decibel values to power values.

Example: [sxx,w] = periodogram(cos(pi./[4;2]*(0:159))'+randn(160,2),"power") specifies the periodogram power spectrum estimate of a two-channel sinusoid embedded in white Gaussian noise and the normalized frequencies at which it is computed.

Data Types: single | double

rbw — Resolution bandwidth

positive scalar

Resolution bandwidth, specified as a positive scalar. The resolution bandwidth is the product of two values: the frequency resolution of the discrete Fourier transform and the equivalent noise bandwidth of the window used to compute the PSD.

Data Types: single | double

Additional Options

freqRange — Frequency range

two-element vector

Frequency range, specified as a two-element vector of real values. If you do not specify freqRange, then medfreq uses the entire bandwidth of the input signal.

Data Types: single | double

sigDomain — Signal domain

"frequency" (default) | "time"

Since R2026b

Signal domain, specified as either of these values:

• "frequency" — The function computes the median frequency from the PSD estimate (pxx) or power spectrum estimate (sxx).

• "time" — The function computes the median frequency from the time-domain input signal x.

  Note

  Specify Domain="time" to distinguish time-domain inputs (x,t) from frequency-domain inputs (pxx,f).

Data Types: char | string

Output Arguments

fMed — Median frequency

scalar | vector

Median frequency, specified as a scalar or vector.

• If you specify a sample rate, then fMed has the same units as Fs.

• If you do not specify a sample rate, then fMed has units of rad/sample.

bandPow — Band power

scalar | vector

Band power, returned as a scalar or vector.

Algorithms

To determine the median frequency, medfreq computes a periodogram power spectrum estimate using a rectangular window. You can obtain the same value of median frequency, fMed, from a signal x at a sample rate Fs in these three ways.

Directly from the signal	fMed = medfreq(x,Fs)
From the periodogram of the signal	[P,F] = periodogram(x,[],length(x),Fs); fMed = medfreq(P,F)
From the power spectral estimate (Welch's PSD) of the signal	[P,F] = pwelch(x,rectwin(length(x)),[],length(x),Fs); fMed = medfreq(P,F)

Note

Because medfreq uses an intermediary representation to transform the input signal from the time domain to frequency domain, the returned median frequency might vary, depending on the signal transformation method, number of DFT points, and window size.

References

[1] Phinyomark, Angkoon, Sirinee Thongpanja, Huosheng Hu, Pornchai Phukpattaranont, and Chusak Limsakul. "The Usefulness of Mean and Median Frequencies in Electromyography Analysis." In Computational Intelligence in Electromyography Analysis – A Perspective on Current Applications and Future Challenges, edited by Ganesh R. Naik. London: IntechOpen, 2012. https://doi.org/10.5772/50639.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

GPU Code Generation
Generate CUDA® code for NVIDIA® GPUs using GPU Coder™.

Thread-Based Environment

Run code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool.

The medfreq function supports thread-based environments with these usage notes and limitations:

• The syntax with no output arguments is not supported.

For more information, see Run MATLAB Functions in Thread-Based Environment.

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2015a

R2026b: Use timetable, sample time, and nonuniform signal inputs

The medfreq function supports new ways to specify time-domain signal inputs:

• Uniformly sampled signals, either as a timetable or as a signal vector or matrix and a sample time.

• Nonuniformly sampled signals, either as a timetable or as a signal vector or matrix paired with a nonuniform time vector.

R2023a: Use gpuArray objects

The medfreq function supports gpuArray objects. You must have Parallel Computing Toolbox™ to use this functionality.