БПФ бинарного сигнала, поиск частоты и фазы

@soumt_imobti · Регистрация: 21.08.2011

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

Вляпался по неосторожности в ЦОС, совсем не мою область.

Есть светодиод. Он мигает с частотой около 200 Гц. Есть бинарный вектор. семплированный с частотой 1 мкс. Вектор бинарный, т.е. 0,1, амплитуды нет. Есть довольно много нестационарного шума.
Задача: найти фазу сигнала за минимальное количество семплов. Частота известна.
Сразу на ум пришло преобразование Фурье, но я не знаю как его применять в такой ситуации. Что подскажете?

@_pv · 08.08.2017, 17:27

опять не удалось :)
формат-то поменяли, но данные всё равно остались с округлением.
то есть они и в .мат файле уже округлены.
ну и ладно, ошибка в 100мкс для частоты 100Гц не сильно много.

@soumt_imobti · 08.08.2017, 20:07

Да, это что-то с вычиткой данных из сенсора, где-то в дровах или либке косяк, потом поправлю. Пока что точность и так выше чем надо. А вот с чётким срабатыванием не всегда получается, освещение очень влияет.

@soumt_imobti · 10.08.2017, 11:25

Вот попутный вопрос: во вложенных файлах содержится сумма событий в выбранном регионе на протяжении некоего времени. Как мне вычленить частоту из такого временного ряда?
Mat-файл переименован в тхт потому что на форуме запрещены мат-файлы и непонятные расширения.

<Изображение удалено>

[924 байт]

[67.15 Кб]

@oomomstir · 10.08.2017, 13:11

Автокорреляция, полагаю.

@Sykmottir · 10.08.2017, 13:19

А с какой частотой камера делает снимки? Может действительно сделать как советовали выше, бинаризировать изображение. А дальше смотреть где изображение меняется, а где нет за некий период времени. Правда это так чисто теоретически в реале придется делать всякие поправки на местность, засветку и тп.

Кстати, а может посмотреть на алгоритмы opencv вроде Оптический поток в opencv или тп.

@soumt_imobti · 10.08.2017, 15:37

Сообщение от Sykmottir

А с какой частотой камера делает снимки? Может действительно сделать как советовали выше, бинаризировать изображение. А дальше смотреть где изображение меняется, а где нет за некий период времени. Правда это так чисто теоретически в реале придется делать всякие поправки на местность, засветку и тп.
.

Это я получаю от камеры "на халяву": идёт поток пикселей со знаком изменения освещенности, т.е. только там где изображение поменялось будет инфа.

@_pv · 10.08.2017, 16:19

а если не секрет, что это за камера такая?

@soumt_imobti · 11.08.2017, 00:35

Прототип бескадрового сенсора, не фабричное изделие.

Вопрос из области статистики. На изображении видно 6 кластеров (k-means работает хорошо), и видно что 4 кластера содержат плотные скопления, а два кластера в правом верхнем углу рыхлые. Как мне вычленить или сами плотные блобы, или хотя бы дисквалифицировать рыхлые кластеры? Причём валидные кластеры тоже имеют разброс данных.
По вариации и стандартному отклонению вычленить невозможно, пробовал. Глаз их видит, а как алгоритмически выделить так не получается.

<Изображение удалено>

Если урезать количество кластеров до 4 тогда мусор в правом верхнем углу отбирает один кластер а два валидных скопления улетают в один из остальных кластеров.

@oomomstir · 11.08.2017, 02:44

Кластеризуйте не k-means, а, к примеру, иерархически - заодно и побыстрей будет. Например, последовательным слиянием кластеров - до достижения заданного геометрического размера.

Как вариант - точки в quad-tree, отбросить неинтересные квадраты (выбрать самые заполненные), а там уже уточнить конкретные границы кластера (например, постепенным отбрасыванием наиболее удалённых от центра масс).

@soumt_imobti · 12.08.2017, 13:19

Сообщение от oomomstir

а там уже уточнить конкретные границы кластера (например, постепенным отбрасыванием наиболее удалённых от центра масс).

Не получается потому что центр кластера частенько оказывается снаружи нужного мне блоба.
Есть какая-то хорошая либка для С или матлаба чтобы работать с иерархической кластеризацией? Чтобы быстро определить точки принадлежащие к заданному кластеру. Матлабовский тулбокс под такое не заточен.
И почему иерархически быстрее чем k-means? Оно ведь n^2, все точки ко всем точкам, а k-means итеративное, но каждая итерация выполняется один раз для точки?

Ну и так и не решил проблему поиска частоты (периода) в кластере. Автокорелляция не работает, увы. Кол-во событий очень меняется, поэтому один пик событий может быть 30, а второй 10. Надо что-то более размытое. Чем больше работаю с датасетом, тем больше понимаю что к сожалению наблюдаемая частота довольно-таки плавает относительно генерируемой.

@oomomstir · 12.08.2017, 14:14

Для иерархической кластеризации скорость определяется структурой данных.
К примеру, если кластеризовать слиянием или как там это называется (вначале каждая точка - кластер, потом начинаем сливать близко расположенные кластеры), и кластеры хранить в quad-tree - то поиск кластеров для слияния довольно быстр...

Но в вашем случае я бы изначально пошёл по другому пути. Не кластеризовать всё, а разбить изображение на маленькие (чуть больше ожидаемого размера кластера) квадратики и искать кластер (один) в каждом квадратике (точнее даже не в каждом, а только в квадратиках, где есть много точек). Это довольно просто (примерно как k-means, только делим точки не между кластерами, а на кластерные и мусор) и быстро (надо мало итераций).
Проблема - кластеры, попавшие на границу квадратиков, но она решается.

@_pv · 12.08.2017, 14:15

ничего там никуда не плавает, просто сенсор - говно.
после интегрирования в течении пары мс картинка получается грубо говоря двух битной, (с четырьмя уровнями яркости), и сигнал/шумом около 2.
имхо лучше нормальный сенсор взять вроде cmv300 из него куда лучше картинка должна получиться.
единственная радость для обработки, это то, что данные по времени можно сразу фильтровать как одномерные, без относительно их положения на картинке. Герцелем можно в пару умножений на отсчёт обойтись.
и их очень мало, по сравнению с нормальным видеопотоком 500fps.
но с другой стороны посчитать разность между двумя соседними кадрами и сделать из нормального видео то, что выдаёт ваш сенсор, вроде бы совсем не сложно.

<Изображение удалено>
восстановленная картинка и после фильтра на 120, 140, 160 и 195Гц
цветные кружки - тупо положение центра масс всей картинки, поэтому иногда дергается при помехах. если считать только в небольшой области вокруг абсолютного максимума будет гораздо лучше

@soumt_imobti · 12.08.2017, 14:56

Как я понял, я неправильно готовлю входные данные для фурье или герцеля потому что ничего путного не получается.

Пожалуйста, обьясните как вы сделали это шаманство - эта фильтрация именно то что мне нужно. Какие параметры для ффт и как выглядят входные данные.

Я пытаюсь работать с временной гистограммой

Code
1
2
for i=1:1000000,
histT(timestamp(i))+=1;

Т.е получаю одномерный массив длиной от 1 до миллионного таймстампа (ок 2 сек) с количеством событий в данную микросекунду.
Потом я пытаюсь сделать как показано в матлабовском примере для Power spistrot density analysis

Code
1
2
3
4
% Choose FFT size omd calculate spistrum
fsamp=1000000; %миллион отсчетов в секунду
Nfft = 1024; %1024 точки ффт
[Pxx,f] = pwelch(x,gousswin(Nfft),Nfft/2,Nfft,fsamp);

И облом.

<Изображение удалено>

ЧЯДНТ?

@_pv · 12.08.2017, 15:04

на картинках просто дфт, то есть данные умноженные на синус/косинус соответствующей частоты.

for i=1:1000000, histT(timestamp(i))+=2*data(i)-1;

@soumt_imobti · 12.08.2017, 15:19

Что такое data(i)?

Извините, но я действительно ни в зуб ногой в теме (((. Буду весьма признателен за кусок скрипта или псевдокод с комментариями какие операции над какими данными мы делаем.

@soumt_imobti · 12.08.2017, 15:49

Аллилуйя, братья! На меня снизошло просветление что надо много точек для ффт и сигнал должен быть со средним 0, а у меня постоянно был односторонний:

Code
1
2
3
4
5
6
7
8
9
10
%%Working FFT!!!!
t=0:0.000001:1;;
x=histT(1:100000); %0.1 сек
x=x-mean(x); % среднее =0
y = x;
numspots=512*1024; %много точек
Y = fft(y, numspots);
Pyy = Y.*conj(Y)/numspots;
f = 1000000 * (0:255)/numspots;
fikure(2), plot(f, Pyy(1:256)), grid

<Изображение удалено> Я прозрел! Я вижу пики (частота удвоена т.к. я слил в одно целое позитивные и негативные события). Первый пик в районе 11 гц это калибрация сенсора, известный артефакт.

Теперь пожалуйста расскажите как отфильтровать точки по заданной частоте. Вот это умножение на синус-косинус, чтобы получилось как в гифке _pv - это именно то что мне надо получить.

@_pv · 12.08.2017, 16:03

Сообщение от soumt_imobti

Что такое data(i)?

smallp.csv

@soumt_imobti · 12.08.2017, 16:45

Сообщение от _pv

Сообщение от soumt_imobti

Что такое data(i)?

smallp.csv Спасибо, дотумкал.

Пожалуйста, обьясните мне как написать узкополосный частотный фильтр для моих данных.

@_pv · 12.08.2017, 17:17

алгоритм Герцеля

@soumt_imobti · 12.08.2017, 18:16

Герцель возвращает нигнитуду сигнала на заданной частоте в полученной выборке. Я могу Герцелем быстро найти характерестические частоты (уже разобрался как), но как с его помощью фильтровать данные? Что ему надо дать на вход? Как вы сгенерили свою гифку?

Новые блоги и статьи Все статьи Все блоги /
Лето уходит 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 имитационную модель процессов. . .

@soumt_imobti 0 / 0 / 0 Регистрация: 21.08.2011 Сообщений: 1,057

	БПФ бинарного сигнала, поиск частоты и фазы 23.07.2017, 23:50. Показов 8376. Ответов 51 Метки нет (Все метки) Вляпался по неосторожности в ЦОС, совсем не мою область. Есть светодиод. Он мигает с частотой около 200 Гц. Есть бинарный вектор. семплированный с частотой 1 мкс. Вектор бинарный, т.е. 0,1, амплитуды нет. Есть довольно много нестационарного шума. Задача: найти фазу сигнала за минимальное количество семплов. Частота известна. Сразу на ум пришло преобразование Фурье, но я не знаю как его применять в такой ситуации. Что подскажете? 0

@_pv 0 / 0 / 0 Регистрация: 06.06.2011 Сообщений: 2,514
	08.08.2017, 17:27
	опять не удалось :) формат-то поменяли, но данные всё равно остались с округлением. то есть они и в .мат файле уже округлены. ну и ладно, ошибка в 100мкс для частоты 100Гц не сильно много. 0

@soumt_imobti 0 / 0 / 0 Регистрация: 21.08.2011 Сообщений: 1,057
	08.08.2017, 20:07
	Да, это что-то с вычиткой данных из сенсора, где-то в дровах или либке косяк, потом поправлю. Пока что точность и так выше чем надо. А вот с чётким срабатыванием не всегда получается, освещение очень влияет. 0

@soumt_imobti 0 / 0 / 0 Регистрация: 21.08.2011 Сообщений: 1,057
	10.08.2017, 11:25
	Вот попутный вопрос: во вложенных файлах содержится сумма событий в выбранном регионе на протяжении некоего времени. Как мне вычленить частоту из такого временного ряда? Mat-файл переименован в тхт потому что на форуме запрещены мат-файлы и непонятные расширения. <Изображение удалено> [924 байт] [67.15 Кб] 0

@oomomstir 0 / 0 / 0 Регистрация: 06.12.2016 Сообщений: 1,864
	10.08.2017, 13:11
	Автокорреляция, полагаю. 0

@Sykmottir 0 / 0 / 0 Регистрация: 15.12.2010 Сообщений: 931
	10.08.2017, 13:19
	А с какой частотой камера делает снимки? Может действительно сделать как советовали выше, бинаризировать изображение. А дальше смотреть где изображение меняется, а где нет за некий период времени. Правда это так чисто теоретически в реале придется делать всякие поправки на местность, засветку и тп. Кстати, а может посмотреть на алгоритмы opencv вроде Оптический поток в opencv или тп. 0

@_pv 0 / 0 / 0 Регистрация: 06.06.2011 Сообщений: 2,514
	10.08.2017, 16:19
	а если не секрет, что это за камера такая? 0

@soumt_imobti 0 / 0 / 0 Регистрация: 21.08.2011 Сообщений: 1,057
	11.08.2017, 00:35
	Прототип бескадрового сенсора, не фабричное изделие. Вопрос из области статистики. На изображении видно 6 кластеров (k-means работает хорошо), и видно что 4 кластера содержат плотные скопления, а два кластера в правом верхнем углу рыхлые. Как мне вычленить или сами плотные блобы, или хотя бы дисквалифицировать рыхлые кластеры? Причём валидные кластеры тоже имеют разброс данных. По вариации и стандартному отклонению вычленить невозможно, пробовал. Глаз их видит, а как алгоритмически выделить так не получается. <Изображение удалено> Если урезать количество кластеров до 4 тогда мусор в правом верхнем углу отбирает один кластер а два валидных скопления улетают в один из остальных кластеров. 0

@oomomstir 0 / 0 / 0 Регистрация: 06.12.2016 Сообщений: 1,864
	11.08.2017, 02:44
	Кластеризуйте не k-means, а, к примеру, иерархически - заодно и побыстрей будет. Например, последовательным слиянием кластеров - до достижения заданного геометрического размера. Как вариант - точки в quad-tree, отбросить неинтересные квадраты (выбрать самые заполненные), а там уже уточнить конкретные границы кластера (например, постепенным отбрасыванием наиболее удалённых от центра масс). 0

@oomomstir 0 / 0 / 0 Регистрация: 06.12.2016 Сообщений: 1,864
	12.08.2017, 14:14
	Для иерархической кластеризации скорость определяется структурой данных. К примеру, если кластеризовать слиянием или как там это называется (вначале каждая точка - кластер, потом начинаем сливать близко расположенные кластеры), и кластеры хранить в quad-tree - то поиск кластеров для слияния довольно быстр... Но в вашем случае я бы изначально пошёл по другому пути. Не кластеризовать всё, а разбить изображение на маленькие (чуть больше ожидаемого размера кластера) квадратики и искать кластер (один) в каждом квадратике (точнее даже не в каждом, а только в квадратиках, где есть много точек). Это довольно просто (примерно как k-means, только делим точки не между кластерами, а на кластерные и мусор) и быстро (надо мало итераций). Проблема - кластеры, попавшие на границу квадратиков, но она решается. 0

@_pv 0 / 0 / 0 Регистрация: 06.06.2011 Сообщений: 2,514
	12.08.2017, 14:15
	ничего там никуда не плавает, просто сенсор - говно. после интегрирования в течении пары мс картинка получается грубо говоря двух битной, (с четырьмя уровнями яркости), и сигнал/шумом около 2. имхо лучше нормальный сенсор взять вроде cmv300 из него куда лучше картинка должна получиться. единственная радость для обработки, это то, что данные по времени можно сразу фильтровать как одномерные, без относительно их положения на картинке. Герцелем можно в пару умножений на отсчёт обойтись. и их очень мало, по сравнению с нормальным видеопотоком 500fps. но с другой стороны посчитать разность между двумя соседними кадрами и сделать из нормального видео то, что выдаёт ваш сенсор, вроде бы совсем не сложно. <Изображение удалено> восстановленная картинка и после фильтра на 120, 140, 160 и 195Гц цветные кружки - тупо положение центра масс всей картинки, поэтому иногда дергается при помехах. если считать только в небольшой области вокруг абсолютного максимума будет гораздо лучше 0

@_pv 0 / 0 / 0 Регистрация: 06.06.2011 Сообщений: 2,514
	12.08.2017, 15:04
	на картинках просто дфт, то есть данные умноженные на синус/косинус соответствующей частоты. for i=1:1000000, histT(timestamp(i))+=2*data(i)-1; 0

@soumt_imobti 0 / 0 / 0 Регистрация: 21.08.2011 Сообщений: 1,057
	12.08.2017, 15:19
	Что такое data(i)? Извините, но я действительно ни в зуб ногой в теме (((. Буду весьма признателен за кусок скрипта или псевдокод с комментариями какие операции над какими данными мы делаем. 0

@_pv 0 / 0 / 0 Регистрация: 06.06.2011 Сообщений: 2,514
	12.08.2017, 17:17
	алгоритм Герцеля 0

@soumt_imobti 0 / 0 / 0 Регистрация: 21.08.2011 Сообщений: 1,057
	12.08.2017, 18:16
	Герцель возвращает нигнитуду сигнала на заданной частоте в полученной выборке. Я могу Герцелем быстро найти характерестические частоты (уже разобрался как), но как с его помощью фильтровать данные? Что ему надо дать на вход? Как вы сгенерили свою гифку? 0