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

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

30.11.2015, 12:13. Просмотров 4034. Ответов 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
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
04.02.2016, 11:19  [ТС] 41
Выброс устранен, шкала дБ. Цветовая шкала - увеличил кол-во "контуров".
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
04.02.2016, 11:35 42
Средняя картинка замечательная! С правой я не понял - почему там масштаб по частоте "поплыл"?
Рекомендую расширить частотный диапазон до частоты Найквиста.
Шкалу рекомендую 0; -120 дБ ограничить (и весь цвет на этот диапазон растянуть)

Добавлено через 7 минут
Причём растянуть нелинейно 0,-3,-6,-12,-24,-48 ....дБ
У каждого дБ - свой цвет
0
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
04.02.2016, 12:09  [ТС] 43
Вот такая картинка получается.
0
Миниатюры
Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
04.02.2016, 12:32 44
Ну и отлично! Надеюсь появление составляющих с отрицательным df/dt не смущает?
Тогда можно переходить к анализу своего сигнала.
Мне очень интересно, что за Х-хромосома там вырисовывается
0
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
04.02.2016, 13:33  [ТС] 45
Надеюсь появление составляющих с отрицательным df/dt не смущает?
Появление не смущает. Смущает, что я не знаю как его объяснить.
Прежде чем перейти к анализу исходного сигнала, хотелось бы узнать следующее:
каким образом Вы заранее "рассчитываете" ширину окна, шаг окна и выбираете оконную функцию?
0
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
04.02.2016, 13:54 46
Начну отвечать с конца.
1. Оконная функция. Выбор то не очень велик.
Fortran
1
2
3
4
5
    er=1                                    !Прямоугольное окно
    er=0.5*(1-cos(ex))                      !Окно Ханна
    er=0.42 -0.5*cos(ex)+0.08*cos(2*ex)     !Окно Блэкмана
    er=0.35875 -0.48829*cos(ex)+0.14128*cos(2*ex)-0.01168*cos(3*ex) !Окно Блэкмана-Харриса
    er=0.25 -0.4825*cos(ex)+0.3225*cos(2*ex)-0.097*cos(3*ex)- 0.008*cos(4*ex) !Flat-top
Можно все перебрать и выбрать лучшее.
Мне больше нравится Окно Блэкмана-Харриса и Flat-top.
2. Шаг окна - любой, пока не надоест ждать окончания работы программы. Уменьшение шага с какого-то момента перестаёт влиять на разрешение картинки.
3. Ширина окна - вот единственный интересный параметр.
Сам поиграй с ним и увидишь появление "тянучек" по вертикали или по горизонтали. Можешь это в своём тестовом примере посмотреть, чтобы принцип понять. Но при переходе к данным с другой дискретизацией и длительностью этот параметр всё равно станет другим!
Цитата Сообщение от Farmazon Посмотреть сообщение
Появление не смущает. Смущает, что я не знаю как его объяснить.
Постараюсь твою картинку подробнее рассмотреть с интерпретацией, но чуть позже.
1
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
04.02.2016, 16:42 47
Давайте подробно рассмотрим картинку, представленную Farmazon
С первой гармоникой всё понятно (100-800 Гц)
Вторая гармоника тоже будет ЛЧМ, но 200-1600 Гц
Третья 300-2400 Гц
А вот с четвёртой гармоникой (и более высокими) гораздо интересней:
4-я гармоника 400-3200 Гц
Как видим она "перепрыгивает" частоту Найквиста FN=2500 Гц
На самом деле из-за эффектов дискретизации (см. теорему Котельникова) составлюющие с
частотами выше Найквиста "отражаются" обратно в область 0-FN.
Мало того, частоты выше частоты дискретизации Fd=5000 Гц отражаются и от
линии F=0. Это видно на гармониках 7 (700-5600 Гц) и выше.

Точно такую же картинку мы бы получили, посылая аналоговый ограниченный сигнал на АЦП без
антиаляйсингового фильтра (фильтра подавляющего все частоты выше Найквиста).

Но у нас тестовый сигнал. Мы его создали искусственно и уже с заданным шагом дискретизации.
То есть просто так мы его уже отфильтровать не сможем.
Технология создания таких тестовых сигналов - в создании сигнала на частоте дискретизации в несколько раз большей требуемой, фильтрации на частоту FN (в нашем случае 2500 Гц) и последующей децимации.
Но еще раз подчеркнуть - это касается создания широкополосных тестовых сигналов.
Если сигнал - просто синус, или ЛЧМ без ограничения амплитуды, то такие усложнения конечно ни к чему

P.S. При переходе к диапазону 0-FN необходимо было подкорректировать ширину окна. Если внимательно посмотреть на картинку при малых временах - линия явно имеет "тянучку" по оси t
1
Миниатюры
Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
04.02.2016, 19:13 48
P.S. Литература по окнам
0
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
05.02.2016, 09:23 49
Продолжаем рассматривать исходный файл.
Сравнение 3-х методов определения f(t)
Вейвлет-анализ проводился на частотах 30-й гармоники целевого сигнала.
0
Миниатюры
Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
05.02.2016, 13:51 50
То есть мой расчёт практически совпадает с Quinn's estimator.
Честно говоря, я с большим недоверием отношусь к подобным способам вычисления "точной" частоты. Всё, в основном будет зависеть от уровня (и спектра) помех.На сигналах с приличным SNR - это ещё может прокатить, а на 30 гармонике вряд ли такой расчёт возможен вообще...
Интересная картинка на 14 МГц.
Мы видим, в принципе, два процесса (не зависимых) - ЛЧМ (это, кстати, 14 гармоника сигнала 1.038МГц) и ЧМ синусом (грубо говоря).
ЧМ синусом, мне больше напоминает работу системы АПЧ... Надо ещё уточнить, что это за гармоника...Но это я могу ошибаться, т.к. по сигналу бывает затруднительно определить аппаратурные особенности
Ну вот практически и всё, что я мог бы сказать по данному вопросу...
Успехов!
1
Миниатюры
Определение частоты квазипериодического сигнала  
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
15.02.2016, 18:25  [ТС] 51
Во-первых очень хочется сказать ОГРОМНОЕ спасибо за помощь в освоении нового алгоритма уважаемому A_Santik'у. Ваш алгоритм был взят на вооружение.
Во-вторых хотелось бы задать следующий вопрос:
1. у меня есть закон изменения частоты во времени - Рис.1, код 1;
C++
1
2
3
4
5
6
7
8
double test2(double x) {
  /*задаются все константы, использующиеся ниже*/
  double d = (x*NuclotronRho*c*particleCharge)/particleMass;
  double beta = d/sqrt(1+pow(d, 2));
  double f = NuclotronKf*c/NuclotronCircumference * beta;
 
  return sin(2*TMath::Pi()*f*x);
}
2. результат FFT - Рис.2(ось Х - Гц);
3. если я принимаю частоту равной 100кГц, то алгоритм отрабатывает замечательно - Рис.3(ось Х - Гц), код 2;
C++
1
2
3
4
double test2(double x) {
  ...
  return sin(2*TMath::Pi()*100000*x);
}
4. Помню свою ошибку с ЛЧМ сигналом и фазой, но тут вроде бы все логично, код 3.
C++
1
2
3
4
5
6
7
8
9
double sgnlBeg = 0.0, sgnlEnd = 0.005;//seconds
int Fd = 100000000;//sample rate
int n = (sgnlEnd-sgnlBeg)*Fd;
//fill histogram with values
  double x = 0;
  for (int i = 0; i <= n; i++) {
    x = (double(i)/n)*(sgnlEnd-sgnlBeg)+sgnlBeg;
    hf->SetBinContent(i+1, test2(x));
  }
Сам вопрос: какого черта у меня появляются частоты после 100кГц?
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
15.02.2016, 21:44 52
Рис.2 - если это сигнал 0-100 кГц, то картинка не верная.
Во- первых частоты больше 100 кГц, непонятный выброс на 200 кГц...
Вообще спектр ЛЧМ 0-100 кГц должен быть симметричен относительно 50 кГц

Добавлено через 2 часа 32 минуты
Цитата Сообщение от Farmazon Посмотреть сообщение
return sin(2*TMath::Pi()*f*x);
Вот это выражение скорее всего неправильное.
U(t)=sin(Ф(t)) вот выражение для сигнала.
http://www.cyberforum.ru/cgi-bin/latex.cgi?\phi (t)=2\pi \int f(t)dt
Т.е. перед тем подставить в аргумент синуса , f(t) - надо проинтегрировать!

Добавлено через 18 минут
Пусть частота постоянна во времени f(t) =F0
http://www.cyberforum.ru/cgi-bin/latex.cgi?\phi (t)=2\pi \int F_0 dt= 2\pi F_0 t
Ведь что такое частота? По определению частота
http://www.cyberforum.ru/cgi-bin/latex.cgi?f(t)=\frac{1}{2\pi} \frac{d\phi (t)}{dt}
1
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
16.02.2016, 14:22  [ТС] 53
Все равно "разрыв шаблона" у меня в голове.
У меня есть промежуток времени t1...t2, пусть на этом промежутке существует n точек. Зная f(t) я получаю массив f[n] значений частоты в каждой точке. И значение времени в каждой точке я тоже знаю, т.е. у меня есть массив t[n]. Далее я просто вычисляю n значений sin(2*Pi*f[i]*t[i]), i = 0..n-1, при этом f[i] однозначно соответствует t[i]. В итоге я получаю функцию(в виде гистограммы), представленную на Рис.1.

Получается, что частота полученной функции неправильна, хмм... Но с другой стороны, что интегрировать понятно, но с каким шагом? У меня и так информация для каждой точки, т.е. точнее быть не может. Уполз размышлять.
0
Миниатюры
Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
16.02.2016, 20:41 54
Цитата Сообщение от Farmazon Посмотреть сообщение
Далее я просто вычисляю n значений sin(2*Pi*f[i]*t[i]), i = 0..n-1, при этом f[i] однозначно соответствует t[i].
И это неправильно!
http://www.cyberforum.ru/cgi-bin/latex.cgi?sin(2\pi \int f(t)dt)
Рекомендую посмотреть, чем ФМ (фазовая модуляция) отличается от ЧМ (частотной модуляции)

Добавлено через 4 часа 26 минут
P.S. Как правильно задать вопрос...
Можно так ( в курилке). "Иваныч, а чем ФМ отличается от ЧМ ?"
Но лучше так (в той же курилке): "Иваныч, у меня здесь, совершенно случайно, завалялись 2 бутылочки коньяка...."

Добавлено через 1 час 33 минуты
P.S. Это я к тому, что вопрос не такой простой, как кажется на первый взгляд ...
0
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
29.02.2016, 18:17  [ТС] 55
Продолжаю свои изыскания.
С моим предыдущим вопросом проблем больше нет - начитался и разобрался.
Хотелось бы поделиться некоторыми наблюдениями относительно 2-х методов: 1) ОПФ+QuinnEstimator(2nd) и 2) Метод, предложенный уважаемым A_Santik'ом.
Наблюдения эти сводятся к определению относительной точности 1 и 2-го методов. Сигнал выглядит как в моем последнем сообщении + "белый шум" с RMS = 10^(-3). Сигнал брался из соображений: максимально приближен к реальности в плане зависимости изменений частоты от времени.
Результаты на Рис.1, 2.
Несложно видеть, что относительная ошибка в случае метода 1 категорически большая по сравнению с методом 2, их разделяют целых 4 порядка! Хотелось бы задать очередной вопрос: а есть ли еще более точные или сравнимые по точности методы определения частоты?

На обоих графиках точки по времени совпадают.
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
29.02.2016, 19:20 56
Цитата Сообщение от Farmazon Посмотреть сообщение
Несложно видеть, что относительная ошибка в случае метода 1 категорически большая по сравнению с методом 2, их разделяют целых 4 порядка!
Мне кажется здесь какая-то ошибка
4 порядка - это слишком круто!
Картинка по методу 2 вообще какая-то странная...
Вопросы:
1. С каким шагом по времени менялось окно в методе 1?
Рекомендую поставить шаг = 1 дискрет
2. Как влияет уровень шума на относительную ошибку?
Желательно снять эти характеристики при разных уровнях шума.

Цитата Сообщение от Farmazon Посмотреть сообщение
...а есть ли еще более точные или сравнимые по точности методы определения частоты?
Такие методы есть.
0
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
01.03.2016, 15:54  [ТС] 57
Поясняющие картинки.
Зависимость частоты от времени - линейна, т.е. f(t) = a*t+b - Рис.1. Сигнал: sin(2*Pi*(a*t^2/2+b*t)). Результат при использовании алгоритма 2 - Рис.2. Сигнал рассматривается на отрезке 0.03-0.031 с, он разбит на 100000 точек, вычисляются значения частоты с 35000-65000 точек, шаг 100 точек для обоих алгоритмов, ширина окна для алгоритма 2 равна 5000 точек. Сигнал после наложения окна Гаусса - Рис.3. Результат Quinn'а на Рис.4. Красная и синяя линии - это диапазон частот, вычисленный с помощью непрерывной функции f(t), за f0 берется средняя частота на данном диапазоне, зеленая линия - это результат Quinn'а. Как-то совсем грустно.
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

Определение частоты квазипериодического сигнала  
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
01.03.2016, 16:07  [ТС] 58
Меняю "белый шум":
Рис.1 - RMS 10^(-1);
Рис.2 - RMS 10^(-2);
Рис.3 - RMS 10^(-3);
Рис.4 - сигнал без шума.
0
Миниатюры
Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала   Определение частоты квазипериодического сигнала  

Определение частоты квазипериодического сигнала  
A_Santik
148 / 129 / 18
Регистрация: 29.04.2015
Сообщений: 626
01.03.2016, 16:20 59
Да, как-то не очень понятно всё...
Чтобы подтвердить/опровергнуть эти данные, мне нужно будет моделировать Quinn'а
Может, всё-таки пока для прояснения обстановки выполнить мои пункты 1 и 2?ъ
И в 4 картинке укажи реальную частоту f для справки (или там треугольниками помечено 580,600,620...кГц?)
Странно, что-то мне кажется даже парабола будет точнее Quinn'а ...
0
Farmazon
9 / 1 / 1
Регистрация: 27.11.2015
Сообщений: 51
01.03.2016, 17:15  [ТС] 60
На 4-й картинке реальная частота обозначена f0 = 636106 Гц. На всех картинках именно так, я - голова садовая забыл об этом написать. Ваше предположение, что парабола будет точнее - тоже пришло в голову, тут сигнал "синус" и оконная функция Гауссиан, все как бы намекает на параболу. Пункты 1 и 2 - в процессе.
0
01.03.2016, 17:15
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
01.03.2016, 17:15

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

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

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


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

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

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