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

@Farmazon · Регистрация: 27.11.2015

Студворк — интернет-сервис помощи студентам

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

@A_Santik · 01.02.2016, 13:01

Ещё один вариант: рассматривать не сам сигнал, а его гармоники. Если изменение частоты присутствует в исходном сигнале, то будет присутствовать и его гармониках, но df/dt будет гораздо больше.
К сожалению на 64k выборок - у меня проблема с выводом результатов. Я Surfer8 использую... 10-15 мин прорисовка картинки занимает

@Farmazon · 01.02.2016, 16:27 **[ТС]**

Я прошу прощения за редкие ответы на возникающие вопросы.
Каждый отсчет соответствует 10 нс.
Я постараюсь предоставить файл в 2 раза больше (~128к точек), но это предел.
Чтобы не быть голословным, выкладываю результат о котором уже упоминал, в отличие от алгоритма, описанного в первом сообщении данной темы, в этом варианте использовался Quinn's estimator (the second estimator) вместо параболы.

@A_Santik · 01.02.2016, 19:07

128к файл не нужен, долго очень придётся считать

Вот не понятно, а что по оси Х отложено на левой картинке?
Давайте, еще раз разберёмся...
Модулирована (ЛЧМ) оказалась частота 1 МГц, как я понял из картинок.
А это не высокочастотка между импульсами! Между импульсами - 15 МГц!
А 1 Мгц соответствуют периоду между импульсами (1-2 и 2-3) где 1-самый "высокий" импульс.
Я "вручную" прикинул - у меня получилось 1030927Гц в конце и 1086956Гц в начале .... Т.е. наклон совсем в другую сторону...

@A_Santik · 02.02.2016, 02:01

...и ещё по выбору целевой частоты. Почему именно 1 МГц?
То что это максимум в спектре - это не совсем логичный критерий. Почему бы не взять локальный максимум (сигнал с частотой 15 МГц)?
Поясню свою мысль - эту частоту можно эффективно отфильтровать полосовым фильтром.
Чем выше частота, тем точнее её можно вычислить за фиксированное время.
Потом выделенную частоту 15 МГц можно коррелировать с ЛЧМ с параметром р0... Правда это долго...
Можно, как в вашем методе попытаться - оконное Фурье с уточнением частоты.

@A_Santik · 02.02.2016, 09:23

Частоту уточнил 18.753 МГц
Вот понаблюдал за ней - у меня она не ЛЧМ!

@A_Santik · 02.02.2016, 10:09

Вот картинка в увеличенном масштабе:

@A_Santik · 02.02.2016, 11:15

Убрал логарифмирование амплитуды

Нам же "артефакты" не надо разглядывать...

@A_Santik · 02.02.2016, 11:30

Параметры расчёта, пока не забыл

:
Данные из файла "data.txt" 64k
Центральная частота окна:

Fortran
1
2
3
4
    do 100  k=18700000,18800000,500 
        Ft=k
        dff=50000 
......

Тип окна: Flat-top

@A_Santik · 02.02.2016, 11:54

Наблюдаем за 2 гармоникой...

@A_Santik · 02.02.2016, 15:19

Вернёмся к частоте 1 МГц.
Интересно, а эта частота линейно модулирована...
Надо посмотреть 14-ю гармонику. (Там вроде пик на FFT)
О, как интересно! Линейная часть - это явно 14 гармоника, а остальной Х - это что???

@Farmazon · 02.02.2016, 16:28 **[ТС]**

Ухх! Не спешите, пож-ста. Я не успеваю за Вами.
Хотелось бы получить тоже самое, но для начала с ЛЧМ сигналом мне нужно разобраться.

ЛЧМ 100 - 800 Гц, 1 секунда, частота дискретизации 5000Гц, сглаживание 100 мс окном Ханна (в начале и в конце).
Ширина частотного окна 2*dff=100 Гц (dff=50), шаг окна 10 Гц, центральная частота окна 10-1000.

1. Вопрос относительно сглаживания окном Ханна в начале и в конце. Реализовал, как показано на Рис. 1. Аналогично и последние 100 мс. Это правильно?
2. У меня ЛЧМ 100 - 800 Гц, но FFT дает результат только до 600 Гц. Брал формулу отсюда https://en.wikipedia.org/wiki/Chirp - Рис. 2. Куда подевались 200 Гц?
3. Конечный результат вот такой, учел Ваши замечания: оси поменял местами, отображаю ln(magnitude) - Рис. 3.
Смущает размазывание в начале и в конце.

P.S. Линейная зависимость частоты от времени была взята в качестве предположения. И полученный результат не является конечным и возможно ошибочен!

@A_Santik · 02.02.2016, 17:14

1. Со сглаживанием что-то не так, судя по спектру....
Спектр с нормальным сглаживанием не имеет таких выбросов на вершине.
Окно на 100 мс уже должно выходить на 1 и оставаться 1 до времени Т-100мс, и далее плавно спадать!
2. чтобы 100-800 получить формула должна быть sin{2*Pi*(100t+700/(2T)*t^2)}, где Т- длительность
700 - это Конечная частота- Начальная
3. "Размазывание" - это однозначно из-за неправильного сглаживания.

Добавлено через 6 минут
ЛЧМ со сглаживанием
f0 -начальная частота (100),df -приращение частоты (700),T00 - длительность, Tcon время сглаживания (0.1 сек)

Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
Subroutine Sweep(f0,df,T00,Tcon)
use PeltonGlobals
real*8 f0,df,xr,t,T00,Tcon,ex,er,Pi
integer*4 i,Ncon
Pi=3.14159265358
do i=1,NN
t=(i-1)/Fd
xr=dsin(2*3.1415926*(f0*t+df*t*t/(2*T00)))
X(i)=cmplx(xr,0.0)
end do
do i=1,NN
t=(i-1)/Fd
if(t .LE. Tcon) then
    ex=Pi*t/Tcon        
    er=0.5*(1.-dcos(ex))    
X(i)=X(i)*er
end if
if(t .GE. (T00-Tcon)) then
    ex=Pi*(2.-(T00-t)/Tcon)     
    er=0.5*(1.-dcos(ex))    
X(i)=X(i)*er
end if
end do
return
end

Добавлено через 15 минут
P.S. NN,Fd, X(i) - глобальные (число точек, частота дискретизации, рабочий массив)

@A_Santik · 02.02.2016, 17:40

Вот такое окно должно быть в начале. Аналогично спадающее в конце.

@Farmazon · 03.02.2016, 10:03 **[ТС]**

Сглаживание на Рис.1, зеленая кривулина сглаживает красную кривулину до 0.1, сглаженная кривулина изображена синим, в конце аналогично, но зеркально.
Результирующая "сигара" на Рис.2.

Мой вопрос насчет частоты разрешился благодаря Вам очень просто - я элементарно забыл про начальную фазу - голова садовая.

@A_Santik · 03.02.2016, 10:20

Ну, по крайней мере теперь видно, что это 100-800 Гц

Всплеск почему-то виден на 0.1 и 0.9....
Амплитуда в логарифмическом масштабе? Попробуй без логарифма - основной сигнал будет лучше видно, а гармоник на сигнале нет пока...
Кстати, получить гармоники просто - ограничь свой сигнал сверху по уровню 0.9 - появятся и чётные и нечётные гармоники.
Потом с параметром dff можно "поиграть" для получения максимального разрешения.

@A_Santik · 03.02.2016, 10:23

P.S. И цветовой масштаб можно изменить...

@Farmazon · 03.02.2016, 11:52 **[ТС]**

Кстати, получить гармоники просто - ограничь свой сигнал сверху по уровню 0.9 - появятся и чётные и нечётные гармоники.

Ограничить сигнал - это значит, что все значения, которые превышают порог 0.9 должны быть исключены из рассмотрения (т.е. n < 5000 станет) или их приравнять просто к 0.9?

@A_Santik · 03.02.2016, 12:04

Те выборки, которые > 0.9 приравнять 0.9.

@Farmazon · 04.02.2016, 10:29 **[ТС]**

Попытки вменяемую шкалу градиента цвета сделать пока не увенчались успехом.
Сигнал ограничил до 0.9 - Рис.1, результат в этом случае Рис.2.

@A_Santik · 04.02.2016, 11:01

Вот я вижу, что логарифмирования уже нет. А зря! Т.к. гармоники имеют очень низкую амплитуду, и их без логарифмирования просто не разглядеть! Обычно берут 20*log10(cdabs(FX(i)), чтобы в дБ получать значения.
cdabs - модуль комплексного числа. log10 - десятичный логарифм.
При таком несимметричном ограничении конечно же появится постоянная составляющая (частота 0 Гц)
Поэтому лучше картинку смотреть с 50 Гц.
Можно конечно сделать симметричное ограничение (-0.9), но тогда в сигнале будут присутствовать только нечётные гармоники (3,5,7...)

Добавлено через 6 минут
P.S. какой-то паразитный выброс разглядел на Т=0.1 - что-то ещё не так со сглаживанием!

Новые блоги и статьи Все статьи Все блоги /
Лето уходит kumehtar 17.08.2026	Мысли в слух kumehtar 17.08.2026 Забавно, насколько сейчас стала доступна информация. Например о магии, духовном развитии, медитациях, и других подобных направлениях, ранее зачастую тайных, передаваемых от учителя к ученику. Хотя. . .	Перемещение строк из ТЧ в другой документ с учетом текущего пробега Maks 17.08.2026 Реализация из решения ниже выполнена на примере нетипового документа "Автозапчасти", с ТЧ "Шины". За основу взят алгоритм отсюда: https:/ / www. cyberforum. ru/ blogs/ 359708/ 10838. html Задача: . . .	Саморегулирующийся социальный контракт для сервера cross-section. Hrethgir 14.08.2026 С кодом конечно таких глубоких размышлений пока не было, впрочем я уже привык к алгоритмизации. Суть предмета записи: снова в диалоге с нейросетью (я взял пока себе ник для учётки админа - Rector). . . .
Часы электронные Uhbif79 12.08.2026 Выкладываю программу часов. Программа позволяет: 1. Использовать системное время и дату, 2. Есть возможность вводить время и дату вручную. 3. Реализованы 2 будильника: начало и конец рабочего дня. . . .	Часы с будильником на основе класса QLCDNumber Uhbif79 12.08.2026 Всем добрый день, выкладываю программу часов с будильником на основе класса QLCDNumber. Здесь я пробовал самостоятельно создавал классы, впервые столкнулся с видимостью переменной одного класса из. . .	Установка MinGW GCC 16.2 и CMake 8Observer8 10.08.2026 VK Видео: https:/ / vkvideo. ru/ video-240781534_456239017 YouTube: eY5-5PyI9NM Текстовая версия	Неделя из жизни имитационной модели склада: мои кривые руки растут, откуда надо anaschu 10.08.2026 Неделя из жизни имитационной модели склада: как я почти написал неправильную логику и что с этим делать Работаю сейчас над учебно-рабочим проектом: строю в AnyLogic имитационную модель процессов. . .

@Farmazon 9 / 1 / 1 Регистрация: 27.11.2015 Сообщений: 51

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

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	01.02.2016, 13:01
	Ещё один вариант: рассматривать не сам сигнал, а его гармоники. Если изменение частоты присутствует в исходном сигнале, то будет присутствовать и его гармониках, но df/dt будет гораздо больше. К сожалению на 64k выборок - у меня проблема с выводом результатов. Я Surfer8 использую... 10-15 мин прорисовка картинки занимает 0

@Farmazon 9 / 1 / 1 Регистрация: 27.11.2015 Сообщений: 51
	01.02.2016, 16:27 [ТС]
	Я прошу прощения за редкие ответы на возникающие вопросы. Каждый отсчет соответствует 10 нс. Я постараюсь предоставить файл в 2 раза больше (~128к точек), но это предел. Чтобы не быть голословным, выкладываю результат о котором уже упоминал, в отличие от алгоритма, описанного в первом сообщении данной темы, в этом варианте использовался Quinn's estimator (the second estimator) вместо параболы. Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	01.02.2016, 19:07
	128к файл не нужен, долго очень придётся считать Вот не понятно, а что по оси Х отложено на левой картинке? Давайте, еще раз разберёмся... Модулирована (ЛЧМ) оказалась частота 1 МГц, как я понял из картинок. А это не высокочастотка между импульсами! Между импульсами - 15 МГц! А 1 Мгц соответствуют периоду между импульсами (1-2 и 2-3) где 1-самый "высокий" импульс. Я "вручную" прикинул - у меня получилось 1030927Гц в конце и 1086956Гц в начале .... Т.е. наклон совсем в другую сторону... 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	02.02.2016, 02:01
	...и ещё по выбору целевой частоты. Почему именно 1 МГц? То что это максимум в спектре - это не совсем логичный критерий. Почему бы не взять локальный максимум (сигнал с частотой 15 МГц)? Поясню свою мысль - эту частоту можно эффективно отфильтровать полосовым фильтром. Чем выше частота, тем точнее её можно вычислить за фиксированное время. Потом выделенную частоту 15 МГц можно коррелировать с ЛЧМ с параметром р0... Правда это долго... Можно, как в вашем методе попытаться - оконное Фурье с уточнением частоты. Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	02.02.2016, 09:23
	Частоту уточнил 18.753 МГц Вот понаблюдал за ней - у меня она не ЛЧМ! Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	02.02.2016, 10:09
	Вот картинка в увеличенном масштабе: Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	02.02.2016, 11:15
	Убрал логарифмирование амплитуды Нам же "артефакты" не надо разглядывать... Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	02.02.2016, 11:54
	Наблюдаем за 2 гармоникой... Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	02.02.2016, 15:19
	Вернёмся к частоте 1 МГц. Интересно, а эта частота линейно модулирована... Надо посмотреть 14-ю гармонику. (Там вроде пик на FFT) О, как интересно! Линейная часть - это явно 14 гармоника, а остальной Х - это что??? Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	02.02.2016, 17:40
	Вот такое окно должно быть в начале. Аналогично спадающее в конце. Миниатюры 1

@Farmazon 9 / 1 / 1 Регистрация: 27.11.2015 Сообщений: 51
	03.02.2016, 10:03 [ТС]
	Сглаживание на Рис.1, зеленая кривулина сглаживает красную кривулину до 0.1, сглаженная кривулина изображена синим, в конце аналогично, но зеркально. Результирующая "сигара" на Рис.2. Мой вопрос насчет частоты разрешился благодаря Вам очень просто - я элементарно забыл про начальную фазу - голова садовая. Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	03.02.2016, 10:20
	Ну, по крайней мере теперь видно, что это 100-800 Гц Всплеск почему-то виден на 0.1 и 0.9.... Амплитуда в логарифмическом масштабе? Попробуй без логарифма - основной сигнал будет лучше видно, а гармоник на сигнале нет пока... Кстати, получить гармоники просто - ограничь свой сигнал сверху по уровню 0.9 - появятся и чётные и нечётные гармоники. Потом с параметром dff можно "поиграть" для получения максимального разрешения. Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	03.02.2016, 10:23
	P.S. И цветовой масштаб можно изменить... 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	03.02.2016, 12:04
	Те выборки, которые > 0.9 приравнять 0.9. 0

@Farmazon 9 / 1 / 1 Регистрация: 27.11.2015 Сообщений: 51
	04.02.2016, 10:29 [ТС]
	Попытки вменяемую шкалу градиента цвета сделать пока не увенчались успехом. Сигнал ограничил до 0.9 - Рис.1, результат в этом случае Рис.2. Миниатюры 0

@A_Santik 148 / 129 / 18 Регистрация: 29.04.2015 Сообщений: 625
	04.02.2016, 11:01
	Вот я вижу, что логарифмирования уже нет. А зря! Т.к. гармоники имеют очень низкую амплитуду, и их без логарифмирования просто не разглядеть! Обычно берут 20log10(cdabs(FX(i)), чтобы в дБ получать значения. cdabs - модуль комплексного числа. log10 - десятичный логарифм. При таком несимметричном* ограничении конечно же появится постоянная составляющая (частота 0 Гц) Поэтому лучше картинку смотреть с 50 Гц. Можно конечно сделать симметричное ограничение (-0.9), но тогда в сигнале будут присутствовать только нечётные гармоники (3,5,7...) Добавлено через 6 минут P.S. какой-то паразитный выброс разглядел на Т=0.1 - что-то ещё не так со сглаживанием! 1