medfreq
Медианная частота.
| Библиотека |
|
Синтаксис
Вызов функции
-
freq, power = medfreq(___,freqRange)— также использует частотный интервалfreqRange, на котором вычисляется медианная частота. Этот синтаксис может включать любую комбинацию входных аргументов из предыдущих синтаксисов, при условии, что второй входной аргумент равен либоFsлибоf. Если второй входной аргумент передан пустым, тоmedfreqиспользует нормализованную частоту. Значение по умолчанию дляfreqRange— вся полоса пропускания входного сигнала.
-
medfreq(___,out=:plot)— строит график спектральной плотности мощности или спектра мощности и аннотирует расчетную медианную частоту.
Аргументы
Входные аргументы
#
x —
входной сигнал
вектор | матрица
Details
Входной сигнал, заданный как вектор или матрица. Если x — вектор, то medfreq обрабатывает его как одиночный канал. Если x — матрица, то функция вычисляет медианную частоту независимо для каждого столбца x.
| Типы данных |
|
#
Fs —
частота дискретизации
положительный вещественный скаляр
Details
Частота дискретизации, заданная как положительный вещественный скаляр. Частота дискретизации — это количество отсчетов в единицу времени. Если единицей времени являются секунды, то частота дискретизации указывается в Гц.
| Типы данных |
|
#
pxx —
спектральная плотность мощности
вектор | матрица
Details
Спектральная плотность мощности, заданная как вектор или матрица. Если pxx — матрица, то medfreq вычисляет медианную частоту независимо для каждого столбца pxx.
Спектральная плотность мощности должна быть выражена в линейных единицах, а не в дБ.
| Типы данных |
|
#
f —
частоты
вектор
Details
Частоты, заданные как вектор.
| Типы данных |
|
#
sxx —
оценка спектральной мощности
вектор | матрица
Details
Оценка спектральной мощности, заданная как вектор или матрица. Если sxx — матрица, то medfreq вычисляет медианную частоту независимо для каждого столбца sxx.
Спектральная плотность мощности должна быть выражена в линейных единицах, а не в дБ.
| Типы данных |
|
#
rbw —
разрешающая способность по частоте
положительный скаляр
Details
Разрешающая способность по частоте, заданная как положительный скаляр. Разрешающая способность по частоте определяется как произведение двух величин: частотного разрешения дискретного преобразования Фурье и эквивалентной ширины полосы шума окна, используемого для вычисления спектральной плотности мощности.
| Типы данных |
|
#
freqRange —
диапазон частот
двухэлементный вектор
Details
Диапазон частот, заданный как двухэлементный вектор вещественных значений. Если freqRange не указан, то medfreq использует всю полосу пропускания входного сигнала.
| Типы данных |
|
Выходные аргументы
#
power —
спектральная мощность в полосе частот
скаляр | вектор
Details
Спектральная мощность в полосе частот, возвращаемая в виде скаляра или вектора.
Примеры
Медианная частота чирп сигнала
Details
Сгенерируем 1024 отсчета чирп сигнала с частотой дискретизации 1024 кГц. Зададим частоту чирпа так, чтобы она начиналась с 50 кГц и достигала 100 кГц в конце сигнала. Добавим белый гауссовский шум таким образом, чтобы отношение сигнал/шум составляло 40 дБ. Перезапустим генератор случайных чисел для получения воспроизводимых результатов.
import EngeeDSP.Functions: medfreq, chirp, std
using Random
nSamp = 1024
Fs = 1024e3
SNR = 40
Random.seed!(0)
t = (0:nSamp-1)./ Fs
x = chirp(t, 50e3, nSamp/Fs, 100e3)
x = x.+ randn(nSamp).* std(x)[1]./(10^(SNR/20))Оценим медианную частоту чирп сигнала. Построим график спектральной плотности мощности и определим медианную частоту.
medfreq(x,Fs,out=:plot)
Сгенерируем еще один чирп сигнал. Зададим начальную частоту 200 кГц, конечную частоту 300 кГц и амплитуду, вдвое превышающую амплитуду первого сигнала. Добавим белый гауссовский шум.
x2 = 2*chirp(t,200e3,nSamp/Fs,300e3)
x2 = x2.+ randn(nSamp).* std(x2)[1]./(10^(SNR/20))Объединим чирп сигналы для получения двухканального сигнала. Оценим медианную частоту каждого канала.
y = medfreq([x x2],Fs)([75008.20087306881, 249993.82312082333], [0.4996644951191848, 1.9996231821604598])Построим графики спектральной плотности мощности двух каналов.
y = medfreq([x x2],Fs,out=:plot)
Сложим два канала, чтобы сформировать новый сигнал. Построим график спектральной плотности мощности и определим медианную частоту.
y = medfreq(x+x2,Fs,out=:plot)
Медианная частота синусоиды
Details
Сгенерируем 1024 отсчета синусоидального сигнала с частотой 100.123 кГц, дискретизированного с частотой 1024 кГц. Добавим белый гауссовский шум таким образом, чтобы отношение сигнал/шум составляло 40 дБ. Перезапустим генератор случайных чисел для получения воспроизводимых результатов.
import EngeeDSP.Functions: medfreq, sin, std
using Random
nSamp = 1024
Fs = 1024e3
SNR = 40
Random.seed!(0)
t = (0:nSamp-1)./ Fs
x = sin.(2π * t * 100.123e3)
x = x.+ randn(nSamp).* std(x)[1]./(10^(SNR/20))Оценим медианную частоту синусоидального сигнала. Построим график спектральной плотности мощности и определим медианную частоту.
medfreq(x,Fs,out=:plot)
Сгенерируем еще один синусоидальный сигнал с частотой 257.321 кГц и амплитудой, в два раза превышающей амплитуду первой синусоиды. Добавим белый шум.
x2 = sin.(2π * t * 257.321e3)
x2 = x2.+ randn(nSamp).* std(x2)[1]./(10^(SNR/20))Объединим синусоиды для получения двухканального сигнала. Оценим медианную частоту каждого канала.
y = medfreq([x x2],Fs)([100005.51887841288, 257099.81517133984], [0.49964321557716784, 0.4998751983807168])Построим графики спектральной плотности мощности двух каналов.
y = medfreq([x x2],Fs,out=:plot)
Сложим два канала, чтобы сформировать новый сигнал. Построим график спектральной плотности мощности и определим медианную частоту.
y = medfreq(x+x2,Fs,out=:plot)
Медианная частота сигналов ограниченных по полосе
Details
Создадим сигнал, спектральная плотность мощности которого напоминает частотную характеристику полосового КИХ-фильтра 88-го порядка с нормированными частотами среза.
import EngeeDSP.Functions: fir1, medfreq, bandpower
d = fir1(88,[0.25 0.45])
medfreq(d,[],[0.3 0.6]*pi, out=:plot)
Выведем медианную частоту и мощность для интервала измерения. Указание частоты дискретизации 2π эквивалентно тому, что частота дискретизации не задана.
mf=medfreq(d,[],[0.3 0.6]*pi)
println("Median = $(round(mf[1]/pi, digits=3))*pi, power = $(round(mf[2]/bandpower(d)*100, digits=1))% of total")Median = 0.371*pi, power = 77.4% of totalДобавим второй канал с нормализованными частотами среза 0.5π рад/отсчет и 0.8π рад/отсчет и амплитудой, составляющей одну десятую от амплитуды первого канала.
d2=fir1(88,[0.5 0.8])/10
medfreq([d d2],[],[0.3 0.9]*pi, out=:plot)
Выведем медианную частоту каждого канала. Разделим на π.
mf=medfreq([d d2],[],[0.3 0.9]*pi)
println(mf[1]/pi)[0.370596807684655, 0.6500040234987472]Алгоритмы
Details
Для определения медианной частоты medfreq вычисляет оценку спектра мощности периодограммы с помощью прямоугольного окна.
Одно и то же значение медианной частоты medFrq можно получить из сигнала x с частотой дискретизации Fs тремя способами.
| Непосредственно из сигнала |
|
| По периодограмме сигнала |
|
| На основе оценки спектральной мощности сигнала (спектральной плотности мощности Уэлча) |
|
Поскольку medfreq использует промежуточное представление для преобразования входного сигнала из временной области в частотную, возвращаемая медианная частота может меняться в зависимости от метода преобразования сигнала, количества точек ДПФ и размера окна.
|
Литература
-
Phinyomark, Angkoon, Sirinee Thongpanja, Huosheng Hu, Pornchai Phukpattaranont, and Chusak Limsakul, The Usefulness of Mean and Median Frequencies in Electromyography Analysis, 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.