Определите частоту из данных сигнала в MATLAB
У меня есть данные с датчика, и мне нужно найти его частоту. Это выглядит как fft() Кажется, это правильный путь, но документы MATLAB показывают только, как получить график частот, я не знаю, что оттуда делать.
Вот как выглядят мои данные:
задан 07 марта ’13, 00:03
2 ответы
Один из способов — действительно использовать fft. Так как fft дает вам частотное представление сигнала, вы хотите найти максимум, а поскольку fft является комплексным сигналом, вам нужно сначала взять абсолютное значение. Индекс будет соответствовать нормированной частоте с максимальной энергией. Наконец, если ваш сигнал имеет смещение, как в случае с тем, который вы показываете, вы хотите избавиться от этого смещения, прежде чем принимать fft, чтобы вы не получали максимум в начале координат, представляющий компонент постоянного тока.
Все, что я описал, помещается в одну строку:
где indexMax — это индекс, по которому можно найти максимальное значение fft.
Примечание: чтобы перейти от indexMax к фактической интересующей частоте, вам необходимо знать длину L fft (такую же, как длина вашего сигнала) и частоту дискретизации Fs. Тогда частота сигнала будет:
В качестве альтернативы, быстрее и довольно хорошо работающей в зависимости от вашего сигнала, возьмите автокорреляцию вашего сигнала:
и найти первый максимум после центральной точки автокорреляции. (Автокорреляция будет симметричной с максимумом посередине.) Найдя этот максимум, вы найдете первое место, где смещенный сигнал выглядит более или менее похожим на себя. Т.е. вы находите период вашего сигнала. Так как сигнал, сдвинутый на кратное своему периоду, всегда будет выглядеть сам по себе, нужно убедиться, что найденный вами максимум действительно соответствует периоду сигнала, а не одному из его кратных.
Из-за шума в вашем сигнале абсолютный максимум вполне может иметь место, кратное вашему периоду, а не самому периоду. Таким образом, чтобы учесть этот шум, вы должны взять абсолютный максимум автокорреляции (автокорреляция(длина(автокорреляция)/2+1), а затем найти, где автокорреляция больше, чем, скажем, 95% от этого максимального значения для первого время во второй половине сигнала 95%, 99% или другое число будет зависеть от того, насколько шум искажает ваш сигнал.
ОБНОВЛЕНИЕ: я понимаю, что я предположил, что вы имели в виду под «частотой» вашего сигнала тональность или базовую гармонику или частоту с наибольшей энергией, как бы вы ни хотели на это смотреть. Если под частотой вы имели в виду частотное представление вашего сигнала, то в первом приближении вы просто хотите построить абс БПФ, чтобы получить представление о том, где находится энергия:
Если вы хотите понять, почему существует абс или какую важную информацию вы теряете, не представляя фазу БПФ, вы можете прочитать немного больше о преобразовании ДПФ, чтобы точно понять, что вы получаете.
Determine frequency from signal data in MATLAB
I have data from a sensor and I need to find the frequency of it. It looks like fft() seems to be the way to go, but the MATLAB docs only show how to get a graph of the frequencies, I don’t know what to do from there.
Here’s what my data looks like:
2 Answers 2
One way to go is indeed to use an fft. Since the fft gives you the frequency representation of the signal, you want to look for the maximum, and since the fft is a complex signal, you will want to take the absolute value first. The index will correspond to the normalized frequency with maximum energy. Last, if your signal has an offset, as is the case with the one you show, you want to get rid of that offset before taking the fft so that you do not get a max at the origin representing the DC component.
Everything I described put in one line would be:
where indexMax is the index where the max fft value can be found.
Note: to get from indexMax to the actual frequency of interest, you will need to know the length L of the fft (same as the length of your signal), and the sampling frequency Fs. The signal frequency will then be:
Alternatively, faster and working fairly well too depending on the signal you have, take the autocorrelation of your signal:
and find the first maximum occurring after the center point of the autocorrelation. (The autocorrelation will be symmetric with its maximum in the middle.) By finding that maximum, you find the first place where the shifted signal looks more or less like itself. I.e. you find the period of your signal. Since the signal shifted by a multiple of its period will always look like itself, you need to make sure that the maximum you find indeed corresponds to the period of the signal and not one of its multiples.
Because of the noise in your signal, the absolute maximum could very well occur at a multiple of your period instead of the period itself. So to account for that noise, you would take the absolute max of the autocorrelation (autocorrelation(length(autocorrelation)/2+1), and then find where the autocorrelation is larger than, say, 95% of that maximum value for the first time in the second half of the signal. 95%, 99%, or some other number would depend on how much noise corrupts your signal.
UPDATE: I realize that I assumed you meant by «frequency» of your signal the pitch or base harmonic or frequency with the most energy, however you want to look at it. If by frequency you meant the frequency representation of your signal, then to a first approximation, you just want to plot the abs of the FFT to get an idea of where the energy is:
If you want to understand why there is an abs, or what relevant info you are losing by not representing the phase of the fft, you may want to read a bit more about the DFT transform to understand exactly what you get.
Определить частоту по данным сигнала в MATLAB
У меня есть данные от датчика, и мне нужно найти частоту его. Похоже на fft() кажется, это путь, но документы MATLAB показывают только, как получить график частот, я не знаю, что делать оттуда.
вот как выглядят мои данные:
2 ответов
один из способов пойти действительно использовать БПФ. Поскольку БПФ дает вам частотное представление сигнала, вы хотите искать максимум, и поскольку БПФ является сложным сигналом, вы захотите сначала принять абсолютное значение. Индекс будет соответствовать нормированной частоте с максимальной энергией. Наконец, если ваш сигнал имеет смещение, как в случае с тем, который вы показываете, вы хотите избавиться от этого смещения перед принятием fft, чтобы вы не получили max в начале координат представление компонента DC.
все, что я описал в одной строке, будет:
где indexMax-это индекс, в котором можно найти максимальное значение fft.
Примечание: чтобы получить от indexMax до фактической частоты интереса, вам нужно будет знать длину L БПФ (такую же, как длина вашего сигнала) и частоту дискретизации Fs. Частота сигнала будет:
альтернативно, быстрее и работает довольно хорошо тоже в зависимости от сигнала взять автокорреляция сигнала:
и найти первый максимум, происходящий после центральной точки автокорреляции. (Автокорреляция будет симметричной с максимумом посередине.) Найдя этот максимум, вы найдете первое место, где смещенный сигнал выглядит более или менее похожим на себя. Т. е. вы находите период вашего сигнала. Поскольку сигнал, сдвинутый на кратное его периоду, всегда будет выглядеть как сам, вы нужно убедиться, что найденный вами максимум действительно соответствует периоду сигнала, а не одному из его кратных.
из-за шума в вашем сигнале абсолютный максимум может очень хорошо произойти в несколько вашего периода вместо самого периода. Поэтому для учета этого шума вы возьмете абсолютный максимум автокорреляции (автокорреляция (длина (автокорреляция)/2+1) , а затем найдете, где автокорреляция больше, чем, скажем, 95% от этого максимального значения впервые во второй половине сигнала. 95%, 99% или какое-либо другое число будет зависеть от того, сколько шума искажает ваш сигнал.
UPDATE: я понимаю, что я предположил, что вы подразумеваете под «частотой» вашего сигнала высоту тона или базовую гармонику или частоту с наибольшей энергией, однако вы хотите посмотреть на это. Если под частотой вы подразумевали частотное представление вашего сигнала, то в первом приближении вы просто хотите построить abs БПФ, чтобы получить представление о том, где энергия:
если вы хотите понять, почему существует abs, или какую соответствующую информацию вы теряете, не представляя фазу БПФ, вы можете прочитать немного больше о преобразовании DFT, чтобы точно понять, что вы получаете.
MATLAB 8.0 (R2012b): создание, обработка и фильтрация сигналов, Signal Processing Toolbox
Окно справки по пакету расширения Signal Processing Toolbox системы MATLAB 8.0 представлено на рис. 1 на фоне рабочего окна самой системы с открытой вкладкой каталога пакетов расширения APPS. Одна из кнопок в панели инструментов Signal Analysis дает доступ к браузеру сигналов, фильтров и спектров, показанному в правой части окна справки. В левой части этого окна указаны наименования разделов пакета Signal Processing Toolbox.
Рис. 1. Окно справки по пакету расширения Signal Processing Toolbox
Как видно на рис. 1, пакет Signal Processing Toolbox состоит из следующих разделов:
- Waveforms — создание сигналов с различной формой и разными законами модуляции;
- Convolution and Correlation — свертка и корреляция сигналов ;
- Transform — преобразование сигналов ;
- Analog and Digital Filters — аналоговые и цифровые фильтры ;
- Spectral Analysis — спектральный анализ.
Создание сигналов
Многие сигналы представлены как функции времени s(t), параметры которой можно изменять с помощью модуляции того или иного вида. Модуляцией называют процесс изменения какого-либо параметра (амплитуды, частоты, фазы и т. д.) по определенному закону, в результате чего сигнал становится переносчиком информации.
MATLAB 8.0 со своими встроенными средствами позволяет создавать множество сигналов. Например, простейшим является синусоидальный сигнал:
s = A sin(2π ft +j),
где A — амплитуда; f — частота; j — фаза сигнала. Задав эти параметры ( t — как вектор отсчетов сигнала, число элементов которого определяет число отсчетов сигнала):
можно легко построить график сигнала s ( t ):
Он показан на рис. 2 в графическом окне системы MATLAB. Такое построение осуществляется в командной строке, на что указывает приглашение к вводу >>. В командном окне также показано окно с краткими данными о системе MATLAB 8.0 (R2012b). Согласно этим данным система выпущена на рынок в 2012 году.
Рис. 2. Окно командного режима
В дальнейшем, учитывая множество сигналов, мы будем представлять их по четыре в каждом графическом окне, используя для этого функцию его разбиения на четыре под-окна subplot.
Рис. 3. Окно с подокнами меандра, апериодического треугольного импульса, пилообразного периодического импульса и периодического треугольного импульса
Представленная программа (она задается в редакторе) дает примеры создания и графической визуализации четырех типов простых сигналов (меандра, одиночного треугольного импульса, пилообразного импульса и симметричного треугольного импульса) (рис. 3):
Рис. 4. Временные диаграммы четырех сигналов, относящихся к функциям Гаусса
Важное значение имеет функция sinc (или sin(π t )/π t при t ≠ 0 и 1 при t = 0). Функция sinc( t ) представляет обратное преобразование Фурье для прямоугольного импульса с высотой 1 и шириной 2π:
Кроме того, эту функцию можно использовать как базисную для восстановления любого сигнала g ( t ) по его отсчетам, если спектр сигнала ограничен условием –p < w < p:
Рис. 5. Получение непрерывной кривой, проходящей через точки (отсчеты) произвольного сигнала
Это положение, вытекающее из известной теоремы Котельникова, иллюстрирует приведенный ниже пример для десяти случайных точек сигнала (рис. 5):
Подобный метод восстановления аналогового сигнала из его цифровых отсчетов ныне применяется во всех цифровых осциллографах.
Рис. 6. Сложные модулированные сигналы
Сложные модулированные сигналы показаны на рис. 6, здесь chirp — период частотно-модулированного сигнала, sawtooth — сигнал с частотной модуляцией по треугольному закону, pulstrain — скачок с линейным спадом, pulstran gauspuls — последовательность импульсов Гаусса с убывающей амплитудой. Для создания этих четырех модулированных сигналов служит программа:
формирует выборку (дискретные значения) косинусоидального сигнала с частотой от f 0 в начальный момент времени t до f 1 в конечный момент времени t 1. Звук такого сигнала напоминает визг, откуда и его название (chirp). По умолчанию t = 0, f 0 = 0 и f 1 = 100. Необязательный параметр phi (по умолчанию 0) задает начальную фазу сигнала. Другой необязательный параметр — method — задает закон изменения частоты: linear — линейный закон (по умолчанию), quadratic — квадратичный и logarithmic — логарифмический.
Функция strips позволяет детально (по частям с длительностью 0,25 с) рассмотреть первые два сложных сигнала, например частотно-модулированного от функции vco. В разделе модуляция/демодуляция есть пример с GUI (рис. 7), иллюстрирующий такие сигналы с различными видами модуляции.
Рис. 7. Визуализация модулированных сигналов
Свертка и корреляция сигналов
Пусть имеется две последовательности, представленные векторами a и b . Сверткой называют одномерный массив, вычисляемый следующим образом:
При записи этого выражения учтено, что нумерация индексов массивов в MATLAB идет с единицы (1). Свертка реализуется функцией conv(a,b):
Операцию свертки часто используют для вычисления сигнала на выходе линейной системы y по сигналу на входе x при известной импульсной характеристике системы h :
Эту операцию можно использовать для осуществления простейшей фильтрации сигнала:
Для двух векторов x и y с длиной m и n определена операция свертки:
Рис. 8. Линейная и циклическая свертка
Обратная свертке функция — это [q,r] = deconv(z,x). Она фактически определяет импульсную характеристику фильтра (рис. 8):
Рис. 9. Очистка от шума линейно-нарастающего сигнала
Свертку часто применяют для очистки сигналов от шума (рис. 9):
Для двумерных массивов также существует функция свертки: Z = conv2(X,Y) и Z = conv2(X,Y,‘option’). Возможна и многомерная свертка — функция convn. Новая функция mscohere строит график зависимости квадрата модуля функции когерентности от частоты (рис. 10):
Рис. 10. Зависимость квадрата модуля функции когерентности от частоты
Следующая программа строит отсчеты двух сигналов с задержкой на три отсчета — треугольного сигнала и шума (рис. 11):
Рис. 11. Два коррелированных сигнала — треугольный и шума
Кросс-корреляцию этих сигналов обеспечивает программа (рис. 12):
Рис. 12. Кросс-корреляция сигналов
Преобразование сигналов
В MATLAB всегда много внимания уделялось различным преобразованиям сигнала [2–5], в частности спектральным. Спектр дискретного сигнала является периодическим, и прямое дискретное преобразование Фурье (ДПФ или Discrete Fourier Transform, DFT) определяется выражением:
Для предотвращения растекания (размазывания) спектра дискретных сигналов часто используются окна. Для этого достаточно в формуле прямого ДПФ под знаком суммы ввести еще один множитель — W ( k ). Соответственно, обратное дискретное преобразование Фурье задается выражением:
ДПФ легко обеспечивает восстановление непрерывных периодических сигналов с ограниченным спектром. Для этого нужно номер отсчета k поменять на нормированное время t / T . Тогда формула восстановления при четном числе отсчетов будет иметь вид:
Для получения полосы частот сигнала от 0 до π/ T приходится смещать нумерацию отсчетов. При нечетном числе отсчетов суммирование ведется при n , меняющемся от –( N –1)/2 до ( N –1)/2. Коэффициенты X ·( n ) с отрицательными номерами вычисляют из соотношения симметрии.
Частотным спектром случайного процесса является преобразование Фурье от корреляционной функции случайного процесса Rx :
В радиоэлектронике особый интерес представляет спектральная оценка сильно зашумленных сигналов. Для таких сигналов применяются два подхода: непараметрический — использующий только информацию, извлеченную из сигнала (реализован в методах периодограмм и Уэлча), и параметрический — предполагающий наличие некоторой статистической модели сигналов, параметры которой подлежат определению. Реализовано восемь классов алгоритмовспектрального анализа : Periodogram, Welch, MTM (Thomson multitaper method), Burg, Covariance, Modified Covariance, Yule-Walker, MUSIC (Multiple Signal Classification) и Eigenvector. Их подробное описание дано в [4–6].
Периодограммы, спектрограммы и их применение
Обычный спектр строится методом быстрого преобразования Фурье (БПФ, FFT) часто с применением временного окна, предотвращающего разрывы сигнала на концах интервала анализа спектра. Например, так реализованы периодограммы. Для построения спектрограмм используется разбивка интервала анализа на короткие окна. Короткое окно пробегает общий интервал анализа, и в каждом частичном интервале строится свой спектр. Их наложение дает спектрограмму, определенную в пространстве «уровень – частота – время» (на плоскости уровень представляется цветом), тогда как обычный спектр определен в плоскости «уровень – частота».
Рис. 13. Построение спектрограмм сигналов, модулированных по различным законам: логарифмическому, линейному, треугольному и синусоидальному
На рис. 13 показано построение четырех спектрограмм для сигналов с различными законами частотной модуляции. Нетрудно увидеть, что во всех случаях закон модуляции отчетливо распознается и позволяет судить об области частот, в которой действует сигнал (chirp или vco). Этот рисунок строит следующая программа:
Рис. 14. Сравнение периодограммы и спектрограммы сигнала sinc(t)
Из сказанного может сложиться неверное представление о явных преимуществах спектрограмм по сравнению с периодограммами (функция psd). То, что это далеко не так, показывает программа для сигнала sinc (рис. 14):
Эта программа строит периодограмму и спектрограмму функции sinc(t). Теперь беспомощной оказывается спектрограмма, по которой ничего нельзя сказать о спектре сигнала и области его частот. А периодограмма более информативна: она указывает на вид спектра и область занимаемых им частот. В частности, хорошо видно постоянство спектра в начальной области частот, присущее этой функции. Правда, вид спектрограммы сильно зависит от типа короткого окна: в данном случае задано окно Блэкмана (Blackman). Вид периодограммы также зависит от выбора окна, но глобального.
Чем сложнее сигнал, тем более детальной и эффектной оказывается спектрограмма. На рис. 15 показан пример модуляции/демодуляции с построением спектрограммы сложного звукового сигнала. Кстати, оригинальный сигнал и сигнал, прошедший модуляцию/демодуляцию, можно воспроизвести на компьютере, оборудованном звуковой картой и акустической системой. Это можно сделать при различных видах модуляции.
Рис. 15. Спектрограмма сложного акустического сигнала с амплитудной модуляцией
Похожие на спектрограммы картинки дают вейвлетограммы и скайлеграммы, получаемые при вейвлет-анализе сигналов [7]. Порою они более информативны. Но вейвлеты в Signal Ptocessing Toolbox не реализованы: они описаны и используются в отдельном пакете расширения Wavelet Toolbox.
Для выполнения дискретного преобразования Фурье различными методами служит GUI-окно Discrete Fourier Transform (рис. 16).
Рис. 16. Спектр пилообразного сигнала в окне дискретного Фурье-преобразования
Оно позволяет задать один из трех видов сигнала (синусоида, меандр и пилообразный) и его периодограмму (спектр) при одном из шести видов окон. Как сигнал, так и окно можно загружать извне из файла или рабочего пространства MATLAB. Нетрудно убедиться в большом влиянии на вид спектра выбранного окна.
Оконные функции и браузер окон
Учитывая важную роль окон, в Signal Processing Toolbox входит 21 N ‑точечное окно ( N — целое число), называемое по фамилии предложивших окно ученых, например Hamming, Blackman, Bartlett, Chebyshev, Taylor, Kaiser и т. д. Для просмотра временных и амплитудно-частотных характеристик всех окон есть соответствующие функции, но удобно пользоваться GUI-браузером окон wvtool.
позволяет строить четыре типа 64‑точечных окон — прямоугольное, Хемминга, Ханна и Гаусса. Можно задать окна и отдельными командами.
Рис. 17. Сравнение в браузере трех окон Кайзера с параметром бетта 1,5, Блэкмана и Блэкмана-Харисса
Это, а также сравнительное построение трех других типов окон обеспечивает следующая программа (рис. 17):
Рис. 18. Конструктор окон
Существует также конструктор-анализатор окон с GUI-интерфейсом, окно которого (рис. 18) открывается командой:
По умолчанию в нем открывается окно объекта sigwin, но можно открыть и другие окна (в том числе окно пользователя) из списка Current Window Information. Открытые окна появляются в третьем нижнем окне. Окна можно скопировать, добавить в список, установить в рабочее пространство MATLAB или стереть.
Изменение числа отсчетов и интерполяция сигналов
Изменение числа отсчетов широко используется в технике цифровой обработки сигналов. Для уменьшения числа отсчетов применяются операция децимации и функция decimate (рис. 19):
Рис. 19. Уменьшение частоты дискретизации сигнала (децимация)
Для увеличения числа отсчетов исходный сигнал интерполируют (функция interp), а затем нужное число отсчетов сигнала берут из кривой интерполированного сигнала (рис. 20):
Рис. 20. Интерполяция сигнала и увеличение частоты его дискретизации
Аналоговые и цифровые фильтры
Фильтры имеют особое значение при обработке сигналов. С их помощью осуществляется очистка сигналов от шума или реализуются избирательные свойства систем. Будем считать, что читатель знаком с теорией фильтров.
Фильтрующие цепи обычно задаются своей операторной передаточной характеристикой:
h ( s ) = a ( s )/ b ( s ).
Имея векторы коэффициентов полиномов a ( s ) и b ( s ), с помощью функции freqs можно построить АЧХ и ФЧХ фильтрующей цепи в логарифмическом масштабе (рис. 21):
Рис. 21. Построение логарифмической АЧХ и ФЧХ системы по ее операторной передаточной характеристике
В Signal Processing Toolboox входит множество функций по расчету и проектированию различных фильтров — нижних, верхних частот и полосовых. Порою для получения важных характеристик фильтров достаточно задать нужную строку программного кода. Например, построение характеристик аналогового фильтра Бесселя (рис. 22) реализуется следующим программным фрагментом:
Для просмотра практически всех характеристик фильтров можно использовать визуализатор фильтров fvtool. Приведенная ниже программа дает характеристики фильтра Бесселя нижних частот (рис. 23):
Рис. 23. Характеристики аналогового эллиптического фильтра нижних частот
Открытое меню Analysis дает представление обо всех доступных видах анализа фильтров. Среди них построение АЧХ и ФЧХ фильтра, импульсная и переходная характеристики, групповая задержка и т. д.
Конструирование двух цифровых фильтров НЧ (FIR и Баттерворта) со сравнительным построением их характеристик в окне fvtool (рис. 24) обеспечивает следующая программа:
Рис. 24. АЧХ двух цифровых фильтров для сравнения
Программа строит и другие характеристики фильтров, например переходные и импульсные. Напомним, что переходная характеристика является реакцией системы (фильтра) на единичный скачок, а импульсная — на импульс единичной площади с длительностью, стремящейся к нулю (рис. 25).
Рис. 25. Импульсные характеристики двух цифровых фильтров
Интерактивный конструктор-анализатор фильтров fdatool
В Signal Proccesing Toolbox входит интерактивный конструктор-анализатор фильтров с GUI-интерфейсом. Он позволяет без какого-либо программирования анализировать и проектировать семь основных типов фильтров (таблица).