Нахождение самого похожего участка в сигнале, БПФ, БПХ и другие методы Опции

[getch]

Приветствую всех!

Я от ЦОС далек, но задача, по-моему, найдет здесь понимание.

Постановка задачи:
1. Есть ВР (временной ряд (вещественный), LEN точек), который с определенной частотой дополняется новым членом.
2. Есть Шаблон - N точек.
3. Надо найти в ВР участок, который наиболее сильно "похож" на Шаблон. Соответственно, в этом участке тоже N точек.

Критерием похожести я выбрал КК (коэффициент корреляции). Т.е. задача сводится к нахождению максимума абсолютного значения КК.

Нашел два оптимизированных алгоритма расчета быстрой корреляции: через БПФ (вещественное) и БПХ (Хартли).

Честно говоря, математическую часть алгоритмов слабо понимаю. И у меня возникло несколько вопросов:
1. Имеет ли смысл для моей задачи (с точки зрения скорости выполнения) использовать алгоритм БПФ для длин векторов, не являющихся степенями двойки? Это работы Бейли, Просекова и другие.
2. БПФ Винограда?
3. Возможно, есть другой эффективный критерий похожести, отличный от КК.
4. Задача, очевидно, стандартная и максимально эффективный велосипед уже давно изобретен. Подскажите, как официально она называется?
5. Какие могут быть оптимизации, чтобы не проделывать полностью вычисления БПФ? Поясню ниже на примере двух случаев:

1. Шаблон - это крайние N точек исходного ВР из Len точек. Пришло новое значение ВР, т.е. его длина стала (Len + 1) точек. Шаблон, соответсвенно "сместился" вправо (изменился) тоже на одно значение. БПФ заново пересчитывать?
2. Шаблон - это какой-то участок из N точек исходного ВР. Если я сдвигаю шаблон на одну точку в какую-то сторону, БПФ надо будет полностью пересчитывать?

И еще один вопрос, если надо найти похожий участок не только на шаблон, но и на "перевернутый шаблон" - точки в обратном направлениии (Было {1, 2, 5, 7}, стало {7, 5, 2, 1}), то можно ли из "нормального" БПФ перейти эффективно к "перевернутому" БПФ?

Уверен, моя косноязычная формулировка задачи вызовет улыбку у форумчан, особенно у гуру. Но это все же хорошие эмоции. Готов учиться и вникать.

P.S. Если возможно, разъясните на более понятном языке, что имелось в виду в этом посте (Сообщение #28)?

[Porty]

подобное делал:
Стабилизация сигнала, т.е. новый кадр осцилограмы должен быть как можно более похожый на старый, чтоб синусойда или голос не бегал как попало, т.е. очень хорошая стабилизация и очень интеллектуальный триггер.

Видео как оно работает:
http://www.youtube.com/watch?v=Q5qWt__qzqY

Сделано так:
Количество выборок Count нового сигнала
Новый сигнал New[i] - он в 2 раза больше нового, чтоб найти похожесть в окне размером в 2*Count
Старый сигнал Old[i] это тот сигнал который был выведен на экран
Массив похожестей Sum[i] - размера нового сигнала;

for j:=0 to Count-1 do
for i:=0 to Count-1 do
sum[j]=sum[j]+Old[i+j]*New[i];

далее нужно найти макс значение в массиве sum[j]; - начиная с него и выводим на экран

Т.е. максимум суммы произведений старого и нового сигнала.

Вроде так.



стабилизация делается в функции find_sync

CODE

unit oscilogramm;
interface
uses windows, sysutils, math, extctrls, graphics, controls, classes;

type
tOSC_mode=(auto, normal, accumulate);
tOSC=class(tobject)
private
bmp:tbitmap;
paintbox:tpaintbox;
data:array of smallint;
old_data:array of smallint;
cursor_x, cursor_y:integer;

procedure SlowFFTPaintBoxMouseDown(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer);
procedure SlowFFTPaintBoxMouseUp(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer);
procedure SlowFFTPaintBoxMouseMove(Sender: TObject; Shift: TShiftState; X, Y: Integer);

public
line_color:tcolor;
x_mult:integer;
data_window:integer;
mode:tOSC_mode;
stop:boolean;

constructor create(v_paintbox:tpaintbox; size:integer);
destructor destroy;override;
procedure add_data(v_data:psmallint; count:integer);
procedure draw;
end;

implementation
uses hsv2rgb;

constructor tosc.create;
begin
inherited create;
paintbox := v_paintbox;

bmp := tbitmap.create;
bmp.Width := paintbox.Width;
bmp.Height := paintbox.Height;

setlength(data, size);
fillchar(data[0], length(data)*sizeof(data[0]), 0);
setlength(old_data, size);

x_mult := 1;
line_color := rgb(0,255,0);
data_window := length(data);
mode := normal;

paintbox.OnMouseDown:=Self.SlowFFTPaintBoxMouseDown;
paintbox.OnMouseMove:=Self.SlowFFTPaintBoxMouseMove;
// paintbox.OnMouseUp:=Self.SlowFFTPaintBoxMouseUp;
end;

destructor tosc.destroy;
begin
bmp.Free; bmp:=nil;
setlength(data,0);

inherited destroy;
end;

procedure tosc.add_data(v_data:psmallint; count:integer);
begin
if self=nil then exit;
if length(data)=0 then exit;
if stop then exit;

count:=min(length(data), count);
if count<length(data) then
move(data[count], data[0], (length(data)-count)*sizeof(data[0]));
move(v_data^, data[length(data)-count], count*sizeof(data[0]));
end;

procedure _inf(var v);
begin
end;

//САМА ПРОЦЕДУРА СТАБИЛИЗАЦИИ!!! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
function find_sync(etalon, where_find:psmallint; count:integer):integer; //САМА ПРОЦЕДУРА СТАБИЛИЗАЦИИ!!! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

function array_mult(a,b:psmallint; count:integer):double;
begin
result:=0.0;
while count>0 do
begin
result:=result+integer(a^)*integer(b^);

inc(a);
inc(;
dec(count);
end;
end;

function array_max_pos(a:pdouble;count:integer):integer;
var
max_val:double;
pos:integer;
begin
pos:=0;
result:=0;
max_val:=MinDouble;
while count>0 do
begin
if a^>max_val then
begin
result:=pos;
max_val:=a^;
end;
inc(a);
inc(pos);
dec(count);
end;
end;

var
results:array of double;
k:integer;
begin
setlength(results, count);

for k:=0 to count-1 do
begin
results[k]:=array_mult(etalon, where_find, count);
dec(where_find);
end;

result:=array_max_pos(@(results[0]), count);

setlength(results, 0);
end;

procedure tosc.draw;
var
maxx,maxy:integer;
pos:integer;
getx:integer;
xadd:integer;
scale:double;
midle:integer;

function gety:integer;
begin
gety:=0;
getx:=getx+xadd;

if pos<0 then exit;
if pos>=length(data) then exit;

gety:=round(midle-scale*data[pos]);
inc(pos);
end;

var
k,r,g,b:integer;
s:string;
line_cnt:integer;
line_pos:integer;
color_cnt:integer;
color_pos:integer;
begin
if self=nil then exit;
if length(data)=0 then exit;
if paintbox=nil then exit;

maxx := paintbox.Width;
maxy := paintbox.Height;
pos := 0;
xadd := min(maxx, max(1, x_mult));
getx := -xadd;
midle := maxy div 2;
scale := (1/(1 shl (8*sizeof(data[0])-1)))*midle*((midle-2)/midle);

bmp.Canvas.Brush.Color := rgb(0,0,0);
bmp.Canvas.Brush.Style := bsSolid;
bmp.Canvas.Pen.Color := rgb(0,0,0);
bmp.Canvas.Pen.Style := psSolid;
bmp.Canvas.Rectangle(-1, -1, maxx+1, maxy+1);


bmp.Canvas.Pen.Color:=rgb(64,64,64);
bmp.Canvas.Pen.Style:=psSolid;
bmp.Canvas.MoveTo(0, midle);
bmp.Canvas.LineTo(maxx, midle);

bmp.Canvas.Pen.Color:=rgb(50, 50, 100);
bmp.Canvas.Pen.Style:=psSolid;
bmp.Canvas.Font.Name:='Courier New';
bmp.Canvas.Font.Size:=12;
bmp.Canvas.Font.color:=rgb(128,128,128);

bmp.Canvas.MoveTo(0, cursor_y);
bmp.Canvas.LineTo(maxx, cursor_y);
bmp.Canvas.MoveTo(cursor_x, 0);
bmp.Canvas.LineTo(cursor_x, maxy);
s:='';
if mode=Normal then s:=s+'Stab';
if mode=Auto then s:=s+'Free';
if mode=accumulate then s:=s+'Summ';
s:=s+' ';
bmp.Canvas.TextOut(maxx-bmp.Canvas.TextWidth(s),0, s);
s:=' '+inttostr(cursor_x)+':'+inttostr(round((midle-cursor_y)/scale));
if cursor_x<>-1 then
bmp.Canvas.TextOut(0, 0, s);

if (mode=normal) or (mode=Auto) then
begin
pos := max(0, length(data)-1-(maxx div xadd));
if mode=normal then
begin
pos := pos - find_sync(@old_data[0], @data[pos], (maxx div xadd)+1);
if pos<0 then pos:=0;
for k:=0 to maxx div xadd do
old_data[k]:=round(old_data[k]*0.9+0.1*data[pos+k]);

//move(data[pos], old_data[0], (maxx div xadd)+1);
end;

bmp.Canvas.MoveTo(0, midle);
bmp.Canvas.Pen.Color:=line_color;
bmp.Canvas.Pen.Style:=psSolid;

bmp.Canvas.MoveTo(getx, gety);
while getx<=maxx do
bmp.Canvas.LineTo(getx, gety);
end;

if mode=accumulate then
begin
line_cnt := maxx div xadd;
color_cnt := length(data) div line_cnt;

for color_pos:=0 to color_cnt-1 do
begin
{ if mode=spectr then
HSVToRGB(1-color_pos/color_cnt, 1, round(255*(color_pos/color_cnt)), r, g,
else}
begin
r:=round(256*(color_pos/color_cnt));
b:=r;
g:=r;
end;
bmp.Canvas.MoveTo(0, midle);
bmp.Canvas.Pen.Color:=rgb(r,g,;
bmp.Canvas.Pen.Style:=psSolid;
pos:=line_cnt*color_pos;

bmp.Canvas.MoveTo(0, gety);
for line_pos:=1 to line_cnt-1 do
bmp.Canvas.LineTo(line_pos*xadd, gety);
end;
end;

paintbox.Canvas.Draw(0,0, bmp);
end;

procedure tosc.SlowFFTPaintBoxMouseUp(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer);
begin
cursor_x:=-1;
cursor_y:=-1;
end;

procedure tosc.SlowFFTPaintBoxMouseDown(Sender: TObject; Button: TMouseButton; Shift: TShiftState; X, Y: Integer);
begin
if ssmiddle in Shift then
stop:=not stop;
if not (ssright in Shift) then exit;

if mode=normal then mode:=auto else
if mode=auto then mode:=accumulate else
if mode=accumulate then mode:=normal;
end;

procedure tosc.SlowFFTPaintBoxMouseMove(Sender: TObject; Shift: TShiftState; X, Y: Integer);
begin
if ssleft in Shift then
begin
cursor_x:=x;
cursor_y:=y;
end
else
begin
cursor_x:=-1;
cursor_y:=-1;
end;

end;

end.

Сообщение отредактировал Porty - Sep 2 2011, 11:38

[getch]

подобное делал:
Стабилизация сигнала, т.е. новый кадр осцилограмы должен быть как можно более похожый на старый, чтоб синусойда или голос не бегал как попало, т.е. очень хорошая стабилизация и очень интеллектуальный триггер.

Видео как оно работает:
http://www.youtube.com/watch?v=Q5qWt__qzqY

Спасибо за ответ! Вроде, задачи разные. Но судя по описанию к видео (ниже), у вас то, что нужно.

Алгоритм отлично ищет подобия в сигнале в режиме реального времени и их отображает:
1. При низком уровне сигнала.
2. При очень сложной структуре сигнала, буквы А Б голосом
3. Не загружает процессор (менее 0.7%) даже для сигнала в 200к выборок в сек

Могли бы вы более подробно описать решенную Вами задачу и алгоритм-идею?
За исходники спасибо. Еще не разбирался в них.

По своей задаче поясню на простом примере. Шаблон имеет 1000 точек. Исходный ВР имеет 100 000 точек. Это значит, что шаблон надо "сравнить" (как вариант - через КК) с 99 000 (100 000 - 1000) вариантами. Или же, если шаблон является частью ВР, то с 98 000 вариантами.

Надеюсь, понятно объяснил.

[Porty]

Могли бы вы более подробно описать решенную Вами задачу и алгоритм-идею?
За исходники спасибо. Еще не разбирался в них.

Нужно было сделать так чтоб на моём осциллографе сигнал рисовался красиво а не дёргался как попало и в то же время отображался в режиме реального времении с макс. количеством кадров (до 25фпс) т.е. простой способ когда отрисовывают раз в секунду а остальное время картинка статична не катит никак. Рисовать постоянно не стирая а затухая градиентом цвета переходя от белого до чёрного тоже не катит - выходит сказачная каша.

Оставался единственный способ - найти похожий сигнал чтоб был близко похож.

Сделал так :
1. храним старую картинку, например в 1000 точек,
2. и получаем новые выборки, например 1000 новых выборок.
3. Но так как нужно искать подобие в каком то диапазоне а выводить нужно те же 1000 точек, то берём за новый диапазон побольше точек, например в 2 раза. т.е. 2000 точек.
3. Далее берём нулевое смещение и перемножаем выборки старого и нового между собой а потом их результат складываем:
Старый0 * Новый0 + Старый1*Новый1 + Старый2 * Новый2 + ... + Старый999*Новый999
Это будет первый коэфициент подобия.

4. Вычисляем второй, просто сдвинув индекс нового на 1
Старый0 * Новый1 + Старый1*Новый2 + Старый2 * Новый3 + ... + Старый999*Новый1000
Это уже второй коэфициент подобия

5. Вычисляем остальные так же, просто сдвинув индекс нового на N
Старый0 * Новый(0+N) + Старый1*Новый(1+N) + Старый2 * Новый(2+N) + ... + Старый999*Новый(999+N)
Это уже N коэфициент подобия

6. Ищем в коэфициентах подобия максимальное значение. Их будет 1000 шт.

7. Максимум и будет самым похожим смещением, т.е. если например максимум 123, то максимально похожим будет подмассив в новом массиве из 2000 тысячь точек, под массив от 123 до 1123 будет максимально похожим. Его и рисуем как отстабилизированное изображение. И принимаем за максимально похожий




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

Сообщение отредактировал Porty - Sep 2 2011, 13:22

[getch]

5. Вычисляем остальные так же, просто сдвинув индекс нового на N
Старый0 * Новый(0+N) + Старый1*Новый(1+N) + Старый2 * Новый(2+N) + ... + Старый999*Новый(999+N)
Это уже N коэфициент подобия

6. Ищем в коэфициентах подобия максимальное значение. Их будет 1000 шт.

7. Максимум и будет самым похожим смещением, т.е. если например максимум 123, то максимально похожим будет подмассив в новом массиве из 2000 тысячь точек, под массив от 123 до 1123 будет максимально похожим. Его и рисуем как отстабилизированное изображение. И принимаем за максимально похожий

Похоже, вы вычисляете для каждой точки свой КК в лобовую (~O(N^2)). Хотя это не совсем КК. Ведь КК был бы равен сумме попарных произведений, если у каждого варианта средняя была бы равна нулю и СКО - единице.

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

Спорное утверждение. Ведь КК(шаблон; перевернутый по X шаблон) не равен нулю.

[Porty]

Похоже, вы вычисляете для каждой точки свой КК в лобовую (~O(N^2)). Хотя это не совсем КК. Ведь КК был бы равен сумме попарных произведений, если у каждого варианта средняя была бы равна нулю и СКО - единице.


Спорное утверждение. Ведь КК(шаблон; перевернутый по X шаблон) не равен нулю.

Не спорю что спорное. Только что на практике проверил, нормально, но естественно не работает если сигнал имеет постоянное смещение. Тут нужны фильтры ВЧ.

[Alexey Lukin]

Критерием похожести я выбрал КК (коэффициент корреляции). Т.е. задача сводится к нахождению максимума абсолютного значения КК.
1. Имеет ли смысл для моей задачи (с точки зрения скорости выполнения) использовать алгоритм БПФ для длин векторов, не являющихся степенями двойки? Это работы Бейли, Просекова и другие.

Вам нужны только БПФ со степенями двойки, независимо от длин сигнала. Сигналы дополняются нулями до степеней двойки.

2. БПФ Винограда?

Не нужно. Используйте те стандартные программы, которые у вас имеются. Наиболее быстрые БПФ работают по стандратным алгоритмам Кули-Тьюки либо со смешанным основанием.

3. Возможно, есть другой эффективный критерий похожести, отличный от КК.

Есть, но КК вычисляется наиболее быстро и в целом неплоха.

4. Задача, очевидно, стандартная и максимально эффективный велосипед уже давно изобретен. Подскажите, как официально она называется?

КК.

5. Какие могут быть оптимизации, чтобы не проделывать полностью вычисления БПФ?

Пара БПФ вычисляет значения КК не для одной позиции шаблона, а для целого ряда позиций. Для вычисления одной позиции вам никакого БПФ не нужно, т.к. КК вырождается в одно скалярное произведение.

И еще один вопрос, если надо найти похожий участок не только на шаблон, но и на "перевернутый шаблон" - точки в обратном направлениии (Было {1, 2, 5, 7}, стало {7, 5, 2, 1}), то можно ли из "нормального" БПФ перейти эффективно к "перевернутому" БПФ?

Достаточно взять комплексное сопряжение спектра, это соответствует развороту во времени.

[getch]

Пара БПФ вычисляет значения КК не для одной позиции шаблона, а для целого ряда позиций. Для вычисления одной позиции вам никакого БПФ не нужно, т.к. КК вырождается в одно скалярное произведение.

Так мне же и надо иметь (Len - N) значений КК, чтобы среди них выбрать максимальное абсолютное значение. Т.е. если делать в лоб, то надо вычислить (Len - N) скалярных произведений векторов длины N. Поэтому и интересуюсь быстрой корреляцией через БПФ, БПХ и другие методы.

Достаточно взять комплексное сопряжение спектра, это соответствует развороту во времени.

Можете пояснить более подробно для слабо-разбирающегося в ЦОС?

[ivan219]

Так мне же и надо иметь (Len - N) значений КК, чтобы среди них выбрать максимальное абсолютное значение. Т.е. если делать в лоб, то надо вычислить (Len - N) скалярных произведений векторов длины N. Поэтому и интересуюсь быстрой корреляцией через БПФ, БПХ и другие методы.

Тогда вам потребуется взять БПФ размеров 2 * (Len - N) округлёныый в большую сторону до степени двойки.
Нужно взять БПФ сигнала и БПФ образца после чего сделать комплексно сопряженный БПФ обрасца. Для этого нужно поменять знак его мнимой части. И перемножить эти спектры. После чего взять обратный БПФ в итоге получите 2 * (Len - N) коэфицентов корреляции.
Код могу привести на Delphi

Можете пояснить более подробно для слабо-разбирающегося в ЦОС?
Делается БПФ после чего у всех мнимых величин меняется знак на противоположный.

Сообщение отредактировал ivan219 - Sep 2 2011, 19:13

[getch]

Написал алгоритм быстрой корреляции через БПФ.
Как понял, быстрая корреляция, это то, что и ранее предлагалось - сумма попарных произведений. Но такой критерий похожести явно не работает в общем случае (когда у выборок МО != 0 и СКО != 1), т.к. максимум такой корреляции показывает совсем не там, где надо.
Подскажите, как от быстрой корреляции перейти к Пирсону с той же эффективностью алгоритма O(N * log(N))?

[getch]

Врубился, как сделать. Надо лишь перед БПФ привести шаблон к нулевому среднему и СКО единица.
А после вычисления через БПФ быстрой корреляции делить полученные элементы массива значений КК на СКО соответствующих выборок в сигнале.
При этом вычисление этих СКО дается практически бесплатно - итерационно.

[getch]

Подскажите, есть ли возможность из RealFFT(Signal, SignalLen) получить RealFFT(Signal, SignalLen + 1)?
Т.е. при приходе нового значения сигнала нужно ли проделывать полностью БПФ или же можно из предыдущего результата БПФ получить текущий?

[Alexey Lukin]

Обычно это делается не по одному отсчёту, а по блокам.
Если надо по одному отсчёту — то см. алгоритм Герцеля, он делает это со сложностью N.

[getch]

Герцеля смотрел ранее, но перед тем, как его применить, провел следующий эксперимент:
1. Сделал RealFFT(Signal, SignalLen), в массиве на выходе посмотрел значение элемента с номером SignalLen.
2. Сделал RealFFT(Signal, SignalLen + 1), в массиве на выходе посмотрел значение элемента также с номером SignalLen.
Соответствующие элементы оказались не равны. Т.е. в общем случае получается, что нет никаких совпадений между RealFFT(Signal, SignalLen) и RealFFT(Signal, SignalLen + 1).
Это так и должно быть (ведь тогда Герцель неприменим)? Или же это особенность вещественной реализации БПФ и надо переходить на комплексную, чтобы получить совпадение?

[Alexey Lukin]

Если сигнал поступает по 1 отсчёту, то ни Герцель, ни БПФ вам вообще не нужны. Считайте корреляцию в лоб.