Форум программистов, компьютерный форум, киберфорум
Наши страницы
Цифровая обработка сигналов
Войти
Регистрация
Восстановить пароль
 
 
Рейтинг 4.62/21: Рейтинг темы: голосов - 21, средняя оценка - 4.62
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
1

Определение частоты квазипериодического сигнала

30.11.2015, 12:13. Просмотров 4047. Ответов 98
Метки нет (Все метки)

Здравствуйте, уважаемые форумчане.
Необходимо экспертное мнение касательно следующей задачи: есть сигнал (задан дискретно, весь сигнал в файле data.txt, первый столбец "тики"*), он квазипериодический, т.е. частота сигнала меняется во времени (примем freq(t) = p0 + p1*t). Эту частоту, а точнее закон по которому она меняется необходимо определить.
Что было сделано:
1. Сигнал делился на куски, "живой" пример Рис. 1;
2. К данному куску сигнала применялось окно Блэкмана-Наталла, Рис. 2;
3. После применения оконной функции производилось БПФ - Рис. 3 - маркером отмечен максимум, с которым в последствии и работаем, ось Y - дБ;
4. Произведено фитирование параболой полученного максимума, Рис. 4, ось Y - дБ;
5. После получения коэффициентов параболы, вычисляется положение центра и откладывается на отдельном графике, далее линейный фит, Рис. 5;
*Оси - X представлены в виде безразмерных "тиков", каждый по 10 нс.
Что хотелось бы услышать:
Я не являюсь специалистом в области гармонического анализа, поэтому мне нужен совет: улучшение имеющегося алгоритма или вариант применения любого другого. Заранее спасибо!
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  
Лучшие ответы (1)
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
30.11.2015, 12:13
Ответы с готовыми решениями:

Различие частоты сигнала и частоты его АКФ
Здравствуйте. Я думал, что если взять обычный синусоидальный сигнал, найти его...

Извлечение сигнала заданной частоты из более сложного сигнала
Нужен совет. Предположим, у меня есть некий сложный сигнал. Я знаю, что в него...

БПФ бинарного сигнала, поиск частоты и фазы
Вляпался по неосторожности в ЦОС, совсем не мою область. Есть светодиод. Он...

Осциллограф показывает две разные частоты одного сигнала
Осциллограф GW INSTEK GDS-840C, на вход поданы прямоугольные импульсы с...

Выделение сигнала полученного с SDR в стороне от центральной частоты
Исходные условия: Имеется SDR приёмник на базе USB TV-тюнера. На нём...

98
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
30.11.2015, 13:20 2
Могу посоветовать непрерывный вейвлет-анализ
На первой картинке вейвлет-анализ отфильтрованного следящими режекторными фильтрами сигнала.
На последней - исходный ЛЧМ-сигнал.
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
19.01.2016, 15:34  [ТС] 3
И снова здравствуйте!
Взялся за вейвлет преобразования. Рассматриваю самый простой случай:
сигнал
f(t) = sin(2*pi/w*t) - задан дискретно с шагом 1 на интервале 0..100;
вейвлет "Mexican Hat"
mh(t) = (1-t*t)*exp(-t*t/2).
Результат вейвлет преобразования обозначим wtr(s, tau), где s - scale. Вот тут начинаются вопросы:
1. Мой scale должен пробегать какой диапазон? Т.к. w=50 в синусе мне известна, то и scale должен "болтаться" вокруг этого значения?
2. Прицепил gif'ку, синяя кривулина -> вейвлет, черная -> сигнал, зеленая -> 1/s*f(t)*mh((t-tau)/s). Хотя бы немного на правду похоже?
3. Если я правильно понимаю алгоритм, то далее мне необходимо посчитать интеграл зеленой кривулины для каждого значения tau?
Шаг по s большой исключительно для создания анимации.
0
Изображения
 
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
19.01.2016, 18:09  [ТС] 4
На данный момент получен следующий результат для scale = 1..30, t = 0..50, количество точек 75.
Вопросы остались такие же!
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
20.01.2016, 12:24 5
Советую для начала взять "длинную" синусоиду - периодов 100, лучше 1000
Параметр scale придётся подбирать экспериментально, т.к. невозможно одновременно получить точное время и точную частоту.
Всё-таки возьмите ЛЧМ сигнал - на нём будет значительно проще работать.

Добавлено через 1 час 52 минуты
Я вот так делаю:
1. Делаем FFT исходного сигнала.
2. Обнуляем спектр на отрицательных частотах (преобразование Гильберта)
3. Для каждой частоты fi :
а. умножаем спектр на оконную функцию (с центром на частоте fi) и шириной df (подбираем по максимальной разрешённости картинки)
б. Делаем обратное преобразование Фурье

Добавлено через 2 часа 52 минуты
Отмечу, что при таком подходе выражение для вейвлета вообще не используется. Хотя если перейти из частотной во временную область - получим что-то очень похожее на вейвлет Морле. Окно предпочтительнее выбирать "посложнее": Ханна, Блэкмана или Блэкмана-Харриса.
Основной недостаток метода - очень большое время расчёта. Да и выходной файл может гигантским получиться.
2
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
25.01.2016, 17:51  [ТС] 6
Прежде всего - ОГРОМНОЕ спасибо за советы!
Для начала стал разбираться с преобразованием Гильберта:
1. взял радиоимпульс s(n), n = 0...N-1, N = 256 - Рис. 1;
2. если я правильно понимаю, мне необходимо получить аналитическую функцию моего сигнала, т.е. произвести преобразование Гильберта сигнала s(n):
2.1 беру FFT этого сигнала Рис. 2;
2.2 обнуляем вторую половину спектра и умножаем первую половину на 2;
2.3 беру обратное FFT, которое даст мне комплексное z(n) = s(n)+j*s_ort(n), n = 0...N-1, N = 256;
2.4 строю действительную часть z(n), она совпадает с исходным сигналом s(n), и мнимую часть s_ort(n) - ортогональное дополнение моего исходного сигнала s(n) - Рис. 3, красная линия - исходный сигнал, синяя - ортогональное дополнение.

Радиоимпульс был взят исключительно для простоты, далее буду использовать ЛЧМ-сигнал.
В предположении, что сделал верно, дальше продвинуться не могу, т.к. не понимаю пункт 3 Вашего алгоритма "3. Для каждой частоты fi :...".
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
26.01.2016, 15:35 7
Ну Да, в принципе ты получил самостоятельно почти-что вейвлет Морле...
Но пока еще не подошёл к анализу своего (произвольного) сигнала...
В принципе, ты понял идеологию
Теперь обрати внимание:
Пункт 2 делается для того, чтобы не умножать на окно (в частотной области) ещё и в отрицательных частотах.
Это не занимает особо много времени - но зачем?
Поясню пункт 3.
Можно взять в качестве Fi - любую частоту- и это достоинство метода.
То есть даже если исходник с частотой дискретизации 44100 ГЦ - нет причин разбивать его равномерно на частоты 0- 22050 Гц.
Можно смело взять наиболее "интересный" частотный диапазон. (Это, кстати, очень экономит время
То есть - частоту fi - это ты сам задаёшь - и получаешь A(f,t) - 3D функцию.

Добавлено через 20 часов 58 минут
Farmazon, ждём твоего решения для тестового ЛЧМ- сигнала.
Надеюсь формулу для ЛЧМ найдешь без проблем?
Любую частоту - это я так, для красного словца сказал... В реальности же приходится всегда выбирать некий диапазон частот с оптимальным шагом
2
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
27.01.2016, 15:56  [ТС] 8
Лучший ответ Сообщение было отмечено Витальич как решение

Решение

Для решения еще рано, честно говоря, идеологию я тоже не понимаю. Видимо, мне сильно мешает "врожденная тупизна".
Как я понимаю алгоритм в данный момент:
1. есть сигнал LFM, n = 0...N-1, N = 256, диапазон 0...500 - Рис.1;
2. необходимо перейти в частотную область - используем FFT;
3. избавляемся от отрицательных частот - используем преобразование Гильберта - Рис.2 и Рис.3;
Тут стоит отметить, что полностью понятные шаги алгоритма закончились и начались глупые вопросы:
Получается, что мы перешли в пространство частот и работаем теперь только в нем?
Пункт 2 делается для того, чтобы не умножать на окно (в частотной области) ещё и в отрицательных частотах.
Т.е. окно я выбираю именно в частотной области? Пусть оно будет длительностью 16 отсчетов, кусок частотного спектра в этом окне обозначим для краткости [h], где h - номер шага, в моем случае шагов 113. Соответственно это окно с шагом 1 начинает ползти от 16/2 по моему частотному спектру с шагом 1 отсчет. Рассмотрим, что происходит на итерации номер 0:
у меня есть [0] и к нему я должен применить оконную функцию, к примеру, Ханна. И вот тут у меня тупик: во временной области применение оконной функции очевидно, а как быть в частотной? Плюс она еще и в частотном окне должна быть задана. Взять FFT от оконной функции и "отрезать" от нее отрицательные частоты? Но чтобы взять FFT оконной функции мне необходимо задать ее на каком-то временном диапазоне, на каком?
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
27.01.2016, 19:41 9
Цитата Сообщение от Farmazon Посмотреть сообщение
Получается, что мы перешли в пространство частот и работаем теперь только в нем?
Да!
Работаешь в пространстве частот точно также как и во временной области.
1. Умножаешь вектор FFT [256] на оконную функцию с центром на частоте 10 Гц
Делаешь обратное FFT. Получаешь вектор [256] со значениями А(10,t)
2. Умножаешь вектор FFT [256] на оконную функцию с центром на частоте 20 Гц
Делаешь обратное FFT. Получаешь вектор [256] со значениями А(20,t)
.....
И т.д. (шаг 10 Гц - просто для примера. 1 Гц можно взять )
Частоту ЛЧМ лучше 50-500 взять

Добавлено через 13 минут
... и сгладить начало и конец хотя бы функцией 1-cos^2

Добавлено через 48 минут
Цитата Сообщение от Farmazon Посмотреть сообщение
И вот тут у меня тупик: во временной области применение оконной функции очевидно, а как быть в частотной? Плюс она еще и в частотном окне должна быть задана. Взять FFT от оконной функции и "отрезать" от нее отрицательные частоты? Но чтобы взять FFT оконной функции мне необходимо задать ее на каком-то временном диапазоне, на каком?
Да всё проще гораздо!
Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
! В FX(i) - значения спектра с обнулёнными отрицательными частотами
! N- количество точек
N2=N/2
dff=12  !Ширина частотного окна  12 Гц 
do 100 k=10,110 !Цикл по частоте 10-110 Гц
Ft=k    !Центральная частота окна
do i=1,N2 
ff=(i-1)*Fd/N2
if((ff .GT. (Ft-dff)) .AND. (ff .LT. (Ft+dff))) then  !если текущая частота в интервале (Ft-dff,Ft+dff)
ex=Pi*(ff-Ft-dff)/dff   
!er=0.5*(1-cos(ex)) !Окно Ханна
er=0.35875 -0.48829*cos(ex)+0.14128*cos(2*ex)-0.01168*cos(3*ex) ! Окно Блэкмана-Харриса
X(i)=FX(i)*er  !Умножение спектра на "окно" 
else
X(i)=(0.0,0.0)
end if
end do
 
CALL DFFTCB(N,X,FX1) !Обратное FFT
.......
100 continue
Добавлено через 13 минут
Надеюсь перевод с Фортрана на любой другой язык не составит труда

Добавлено через 1 час 31 минуту
P.S. После обратного FFT надо не забыть обнулить массив Х(i).
2
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
29.01.2016, 17:47  [ТС] 10
Беру LFM t = 0...0.1 s, n = 0...512, freqMin = 50 Hz, freqMax = 800 Hz - Рис.1,2;
Следую алгоритму, ширина окна ~40 Hz, центральная частота окна ~40 Hz, результат 2-х итерация - Рис.3,4;
Амплитуда после применения оконной функции для итерации с Рис.3 - Рис.5;
На выходе после ОБПФ - хрень.
В каком месте я - дурак?
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
30.01.2016, 10:27 11
Немного не понял, а что ты называешь итерациями???
1. Ну надо с выбранным окном по всему диапазону интересному пройтись 50-800 с шагом 5-10 Гц
2. Смотреть после обратного FFT надо на модуль FX1

Добавлено через 8 минут
т.е . у тебя центральная частота окна должна быть 50-800 Гц с шагом 5 Гц

Добавлено через 15 часов 56 минут
Farmazon, можешь свой исходный файл data.txt выложить? Потом сравним результаты...
0
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
30.01.2016, 13:51 12
Вот как у меня выглядит тестовый сигнал
ЛЧМ 100 - 800 Гц, 1 секунда, частота дискретизации 5000Гц, сглаживание 100 мс окном Ханна (в начале и в конце).
Ширина частотного окна 2*dff=100 Гц (dff=50), шаг окна 10 Гц, центральная частота окна 10-1000.
1
Миниатюры
Определение частоты квазипериодического сигнала  
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
31.01.2016, 17:51  [ТС] 13
Спасибо! Опередили меня. Красота! Очень завидно. У меня все скромнее, технические моменты сильно тормозят.
Вот результат для моего "случая".
Попробую повторить для Вашего сигнала.
1
Миниатюры
Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
31.01.2016, 20:14 14
Не, всё нормально, я тоже с таких картинок начинал (и еще хуже были )
1. Ты не стал на сигнал "окно" накладывать. Поэтому у тебя (даже на спектре) огромные выбросы.(Эффект Гиббса, мать его!)
Но это тестовая задача - можно конечно и не делать... Но это именно эталонный сигнал! Можно и о его "красоте" позаботиться.
2. Странное расположение осей.... Обычно удобнее нах на оси х время держать.
3. Картинку лучше строить в 20*Ln(|FX(i)|) - некие интересные "артефакты" можно разглядеть.
4. Теперь самое интересное осталось - подбор ширины окна.
Там принцип простой - или теряешь в разрешении по времени, или по частоте!
1
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
31.01.2016, 22:00 15
P.S. В чём смысл правильного выбора осей? Потом же результат придётся сравнивать вот с этим графиком:
...только там будет А(f,t)
0
Миниатюры
Определение частоты квазипериодического сигнала  
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
31.01.2016, 23:49  [ТС] 16
Спасибо за замечания, возникло несколько вопросов, но нужно самому сначала "осмыслить".
Полагаю, что возиться мне в ближайшее время стоит именно с тестовым сигналом, потом уже браться за "живой".
Файл с "живыми" данными data.txt, которые использовались при создании темы и из-за которых собственно весь сыр-бор и случился.
Структура файла такая: 1-й столбец - время(каждый отсчет 10 нс), 2-й столбец амплитуда сигнала.
Насчет 10 нс засомневался, но уточнить смогу только завтра.
0
Вложения
Тип файла: zip data.zip (246.8 Кб, 3 просмотров)
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
01.02.2016, 08:33 17
А я наверное неправильно понял исходную задачу...
Пока точно не знаю частоту дискретизации данных все шкалы - относительные.
Тупо взял 5000 выборок, частоту дискретизации 5000Гц
Между пиками - частота не меняется практически! Или я не вижу изменений в таком масштабе?
Я правильно понял, надо искать закон изменения частоты сигнала между пиками?
Или частота прихода "пиков" изменяется по определённому закону?
1
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
01.02.2016, 10:01 18
Вот картинка в другом масштабе:
1
Миниатюры
Определение частоты квазипериодического сигнала  
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
01.02.2016, 10:03  [ТС] 19
Вы все правильно делаете, это задача сама по себе дурная немного.
Между пиками - частота не меняется практически! Или я не вижу изменений в таком масштабе?
Совершенно верно! Частота меняется очень медленно, насчет масштаба затрудняюсь ответить. Если Вы посмотрите на последнюю картинку в моем первом сообщении, там параметр p1 = 8.9 x 10-5. С помощью Quinn's estimator получилось, что p1 = 1 x 10-5, т.е. порядок 10-5 подтверждение нашел.
0
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
01.02.2016, 11:31 20
Тупо взял из исходного файла первые и последние 5000 выборок.
Пока разницы не заметно
Надо гораздо больше точек брать - разрешение по частоте должно существенно увеличиться.
Частота дискретизации точно 100МГц ?
0
Миниатюры
Определение частоты квазипериодического сигнала  
01.02.2016, 11:31
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
01.02.2016, 11:31

Определение доминантной несущей частоты в числовом ряде (БПФ, автокорреляция)
Добрый день, в зарубежной литературе встретил приведенный ниже код для...

Алгоритм поиска пиков (пик-детектор) и определение частоты на малой длительности аудиосигнала
Hallo Forum Benutzern :) Подскажите, какие сейчас есть современные алгоритмы /...

Определение частотного сдвига сигнала после QPSK модуляции
Здравствуйте! Нужна очень ваша помощь. Я начинаю только ознакомление с...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
20
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Рейтинг@Mail.ru