выделение основных частот сигнала
Alex Mizrahi
подозреваю, состоит из нескольких основных синусоид и какого-то шума.
Я хочу выделить эти самые основные синусоиды.
Формально я это понимаю так -- есть последовательность {x_i}, i = -N..N.
Hужно найти такие параметры a, f, p, что вектор {y_i}
y_i = x_i - a * sin(f * (Pi*i/N) + p)
имеет найменьшую норму, то есть, минимизируем сумму квадратов отклонений от
синусоиды с амплитудой a, частотой f и фазой p. y_i как раз и будет "шум",
из которого потом при желании можно выделить другие синусоиды.
Hасколько я понимаю, тема спектрального анализа вообще очень обширная и
сложная. Hо меня пока что интересует самое простое решение, а потом
посмотрим... Так что идти в библиотеку и читать книги не предлагайте :)
Пока что я сделал вот что -- аналитически через решение уравнений нашёл
параметры синусоиды для идеального случая (если она там одна и нет шума). Hо
этот метод принимает шум за сигнал и пытается его интерпретировать, так что
получается неточно, хотя как первое приближение сойдёт.
Вот разложение сигнала через дискретное преобразование Фурье:
http://i.imgur.com/V3eoC.png (http://imgur.com/V3eoC.png если не работает)
Здесь "Ряд 1" -- исходные данные, а "Ряд 4" -- моя аналитическая
апроксимация.
Интуитивно у меня появилась идея "выпрямить" неровности на
Фурье-преобразованном фрагменте, потом преобразовать его обратно и искать
опять же аналитический ответ. Hо не совсем понятно как именно их выпрямлять,
есть какие-то идеи, но наверняка ж не я первый это делаю и есть уже готовые
решения...
Nickita A Startcev
11 Jun 10 , 18:38 Alex Mizrahi ����� � All:
AM> ���������� � ���� ��������� ���� "���������" ���������� ��
AM> �����-��������������� ���������, ����� ������������� ��� ������� �
AM> ������ ����� �� ������������� �����. H� �� ������ ������� ��� ������
AM> �� ����������, ���� �����-�� ����, �� ��������� � �� � ������ ���
AM> ����� � ���� ��� ������� �������...
����� ������� ���:
������� ���, ��������� ������ �������� (����� ��������� ��������) ��������
����� ������ ���.
����������� �������� ������ � ������� � ��������� ���� �������, � ����� �����
�������� �� (�������), ������������ �� ������������ ���������.
� ��������, a*sin(wt)+b*cos(wt) �������� � ������ ������ (��������) � ������
����� � ����������, �� ��� �� �������� - ��� � ����� ������� ���������.
������� ���������� ��������� �� ��������� �������, � ����ݣ���� �������� ����
���������.
ps: ��� � �������, ���� ��� ������ ������� - ��� ������������ �� ������.
. � ���������, H�����.
icq:240059686, lj-user:nicka_startcev
... � � ������ ������ � ������ WDT.
Alex Mizrahi
NAS> О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫, О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ (О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫)
NAS> О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫.
NAS> О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫, О©╫
NAS> О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫ (О©╫О©╫О©╫О©╫О©╫О©╫О©╫),
NAS> О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫.
О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫?
О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫вёО©╫О©╫О©╫О©╫
A := int(f(x)*sin(f*x), x=-Pi..Pi);
B := int(f(x)*cos(f*x), x=-Pi..Pi);
О©╫ О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ A^2 + B^2.
О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫ О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫.
О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫, О©╫О©╫О©╫О©╫ О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫, О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫
О©╫О©╫О©╫О©╫!
О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫ чёО©╫ -- О©╫.О©╫. О©╫ О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫
О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫, О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ (f < 1) О©╫
О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ A0 -- О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ (О©╫О©╫ О©╫О©╫ DC).
О©╫О©╫О©╫ О©╫О©╫О©╫О©╫ A0 О©╫ О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫ О©╫О©╫О©╫О©╫.
О©╫ О©╫О©╫ О©╫О©╫О©╫О©╫, О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫, О©╫О©╫
О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫ О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫ О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫
О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫.
NAS> ps: О©╫О©╫О©╫ О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫, О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ - О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫
О©╫О©╫О©╫О©╫О©╫О©╫.
О©╫О©╫О©╫-О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫ О©╫О©╫О©╫ О©╫О©╫ О©╫О©╫О©╫О©╫О©╫О©╫О©╫О©╫ :(
О©╫ О©╫О©╫О©╫О©╫...
Nickita A Startcev
14 Jun 10 , 13:32 Alex Mizrahi писал к Nickita A Startcev:
NAS>> Можно сделать так:
NAS>> сделать БПФ, посчитать модули амплитуд (сумму квадратов
NAS>> проекций) выделить самый высоки пик. Сворачивать исходный сигнал
NAS>> с синусом и косинусом этой частоты, а потом чуток шевелить ее
NAS>> (частоту),
NAS>> ориентируясь на максимизацию амплитуды.
AM> Амплитуду?
Интенсивность суммы квадратов свёрток (то что у тебя ниже обозначено как
A^2+B^2).
AM> Я считал интеграл свёрток
AM> A := int(f(x)*sin(f*x), x=-Pi..Pi);
AM> B := int(f(x)*cos(f*x), x=-Pi..Pi);
AM> А потом максимизировал A^2 + B^2.
не понял. f(x) и f*x - это разные f?
AM> Максимум оказался совсем не на нужной частоте.
AM> Более того, если у вектора съехавшая фаза, то максимум вообще
AM> находится в нуле!
AM> Я подозреваю проблема вот в чём -- т.к. у меня период функции больше
AM> чем исследуемый интервал,
Ага. Проблема именно в этом. Hа исследуемом интервале должно помещаться
несколько периодов.
NAS>> ps: сам я плотник, если где сильно прогнал - жду комментариев от
NAS>>
AM> спецов.
AM> Что-то мне кажется тут практически никого уже не осталось :(
AM> А жаль...
Посмотри в ЖЖ. Там Санитар Женя часто бывает.
. С уважением, Hикита.
icq:240059686, lj-user:nicka_startcev
... Добросовестно и усердно веселиться?!
Alex Mizrahi
AM>> A := int(f(x)*sin(f*x), x=-Pi..Pi);
AM>> B := int(f(x)*cos(f*x), x=-Pi..Pi);
AM>> А потом максимизировал A^2 + B^2.
NAS> не понял. f(x) и f*x - это разные f?
Ага :)
f(x) -- исследуемый фрагмент сигнала, f -- частота.
AM>> Максимум оказался совсем не на нужной частоте.
AM>> Более того, если у вектора съехавшая фаза, то максимум вообще
AM>> находится в нуле!
AM>> Я подозреваю проблема вот в чём -- т.к. у меня период функции больше
AM>> чем исследуемый интервал,
NAS> Ага. Проблема именно в этом. Hа исследуемом интервале должно
NAS> помещаться несколько периодов.
Hу да, я так понимаю, чем их больше (то есть чем выше частота сигнала) тем
больше точность.
Hо DC будет ненулевой для любой нецелой частоты, так что мне не совсем
понятно какова будет точность у этого метода.
То есть если частота не целая, то период не кратен длине отрезка и на
отрезке есть неуравновешенный "хвост" колебания.
Хвост не больше периода, так что на частотах 10<f<11 хвост будет не больше
1/20.
Hо на частоте f=1.5 хвост равен половине колебания! А эта частота уже выше
чем первая частота Фурье (f=1.0), так что не ясно где граница, после которой
уже хорошо.
С другой стороны, мой аналитический метод даёт погрешность не более одного
процента даже на низкой частоте f = 0.54.
Метод, кстати, тоже использует коэффициенты разложения Фурье -- я вычислил
интегралы свёрток:
A1 = int(a*sin(f*x+p) * sin(x), x = -Pi..Pi);
B1 = int(a*sin(f*x+p) * cos(x), x = -Pi..Pi);
И среднее A0:
A0 = int(a*sin(f*x+p),x = -Pi..Pi);
Т.о. имеем три неизвестных (f, a, p) и три уравнения, так что параметры
находятся однозначно.
Hо этот метод не учитывает шум вообще...
Может быть стоит использовать больше уравнений (коэффициенты A_i, B_i) и
решать их через М.H.К.
По-идее будет более точно, хотя физический смысл этого не вполне очевиден.
Hо там ещё нужно понять как это решить, с тригонометрией работать не так-то
просто...
Хотя, в принципе, оно мне не на столько уж и надо, попробую зайти с другого
бока...
AM>> Что-то мне кажется тут практически никого уже не осталось :(
AM>> А жаль...
NAS> Посмотри в ЖЖ.
NAS> Там Санитар Женя часто бывает.
Да вот я как раз подумал что это задача по его профилю :)