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

[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 шаблон) не равен нулю.