getch 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 · Жалоба Приветствую всех! Я от ЦОС далек, но задача, по-моему, найдет здесь понимание. Постановка задачи: 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 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 (изменено) · Жалоба подобное делал: Стабилизация сигнала, т.е. новый кадр осцилограмы должен быть как можно более похожый на старый, чтоб синусойда или голос не бегал как попало, т.е. очень хорошая стабилизация и очень интеллектуальный триггер. Видео как оно работает: http://www.youtube.com/watch?v=Q5qWt__qzqY Сделано так: Количество выборок Count нового сигнала Новый сигнал New - он в 2 раза больше нового, чтоб найти похожесть в окне размером в 2*Count Старый сигнал Old это тот сигнал который был выведен на экран Массив похожестей Sum - размера нового сигнала; for j:=0 to Count-1 do for i:=0 to Count-1 do sum[j]=sum[j]+Old[i+j]*New; далее нужно найти макс значение в массиве sum[j]; - начиная с него и выводим на экран Т.е. максимум суммы произведений старого и нового сигнала. Вроде так. стабилизация делается в функции find_sync 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(B); 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, B) 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,B); 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. Изменено 2 сентября, 2011 пользователем Porty Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
getch 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 · Жалоба подобное делал: Стабилизация сигнала, т.е. новый кадр осцилограмы должен быть как можно более похожый на старый, чтоб синусойда или голос не бегал как попало, т.е. очень хорошая стабилизация и очень интеллектуальный триггер. Видео как оно работает: 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 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 (изменено) · Жалоба Могли бы вы более подробно описать решенную Вами задачу и алгоритм-идею? За исходники спасибо. Еще не разбирался в них. Нужно было сделать так чтоб на моём осциллографе сигнал рисовался красиво а не дёргался как попало и в то же время отображался в режиме реального времении с макс. количеством кадров (до 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 будет максимально похожим. Его и рисуем как отстабилизированное изображение. И принимаем за максимально похожий Добавлено - если нужен максимально похожий но перевернутый относительно оси Х то нужно просто найти в массиве подобия не максимально большой а максимальной малый коэффициент. Изменено 2 сентября, 2011 пользователем Porty Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
getch 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 · Жалоба 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 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 · Жалоба Похоже, вы вычисляете для каждой точки свой КК в лобовую (~O(N^2)). Хотя это не совсем КК. Ведь КК был бы равен сумме попарных произведений, если у каждого варианта средняя была бы равна нулю и СКО - единице. Спорное утверждение. Ведь КК(шаблон; перевернутый по X шаблон) не равен нулю. Не спорю что спорное. Только что на практике проверил, нормально, но естественно не работает если сигнал имеет постоянное смещение. Тут нужны фильтры ВЧ. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Alexey Lukin 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 · Жалоба Критерием похожести я выбрал КК (коэффициент корреляции). Т.е. задача сводится к нахождению максимума абсолютного значения КК. 1. Имеет ли смысл для моей задачи (с точки зрения скорости выполнения) использовать алгоритм БПФ для длин векторов, не являющихся степенями двойки? Это работы Бейли, Просекова и другие. Вам нужны только БПФ со степенями двойки, независимо от длин сигнала. Сигналы дополняются нулями до степеней двойки. 2. БПФ Винограда? Не нужно. Используйте те стандартные программы, которые у вас имеются. Наиболее быстрые БПФ работают по стандратным алгоритмам Кули-Тьюки либо со смешанным основанием. 3. Возможно, есть другой эффективный критерий похожести, отличный от КК. Есть, но КК вычисляется наиболее быстро и в целом неплоха. 4. Задача, очевидно, стандартная и максимально эффективный велосипед уже давно изобретен. Подскажите, как официально она называется? КК. 5. Какие могут быть оптимизации, чтобы не проделывать полностью вычисления БПФ? Пара БПФ вычисляет значения КК не для одной позиции шаблона, а для целого ряда позиций. Для вычисления одной позиции вам никакого БПФ не нужно, т.к. КК вырождается в одно скалярное произведение. И еще один вопрос, если надо найти похожий участок не только на шаблон, но и на "перевернутый шаблон" - точки в обратном направлениии (Было {1, 2, 5, 7}, стало {7, 5, 2, 1}), то можно ли из "нормального" БПФ перейти эффективно к "перевернутому" БПФ? Достаточно взять комплексное сопряжение спектра, это соответствует развороту во времени. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
getch 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 · Жалоба Пара БПФ вычисляет значения КК не для одной позиции шаблона, а для целого ряда позиций. Для вычисления одной позиции вам никакого БПФ не нужно, т.к. КК вырождается в одно скалярное произведение. Так мне же и надо иметь (Len - N) значений КК, чтобы среди них выбрать максимальное абсолютное значение. Т.е. если делать в лоб, то надо вычислить (Len - N) скалярных произведений векторов длины N. Поэтому и интересуюсь быстрой корреляцией через БПФ, БПХ и другие методы. Достаточно взять комплексное сопряжение спектра, это соответствует развороту во времени. Можете пояснить более подробно для слабо-разбирающегося в ЦОС? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
ivan219 0 2 сентября, 2011 Опубликовано 2 сентября, 2011 (изменено) · Жалоба Так мне же и надо иметь (Len - N) значений КК, чтобы среди них выбрать максимальное абсолютное значение. Т.е. если делать в лоб, то надо вычислить (Len - N) скалярных произведений векторов длины N. Поэтому и интересуюсь быстрой корреляцией через БПФ, БПХ и другие методы. Тогда вам потребуется взять БПФ размеров 2 * (Len - N) округлёныый в большую сторону до степени двойки. Нужно взять БПФ сигнала и БПФ образца после чего сделать комплексно сопряженный БПФ обрасца. Для этого нужно поменять знак его мнимой части. И перемножить эти спектры. После чего взять обратный БПФ в итоге получите 2 * (Len - N) коэфицентов корреляции. Код могу привести на Delphi Можете пояснить более подробно для слабо-разбирающегося в ЦОС? Делается БПФ после чего у всех мнимых величин меняется знак на противоположный. Изменено 2 сентября, 2011 пользователем ivan219 Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
getch 0 8 сентября, 2011 Опубликовано 8 сентября, 2011 · Жалоба Написал алгоритм быстрой корреляции через БПФ. Как понял, быстрая корреляция, это то, что и ранее предлагалось - сумма попарных произведений. Но такой критерий похожести явно не работает в общем случае (когда у выборок МО != 0 и СКО != 1), т.к. максимум такой корреляции показывает совсем не там, где надо. Подскажите, как от быстрой корреляции перейти к Пирсону с той же эффективностью алгоритма O(N * log(N))? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
getch 0 9 сентября, 2011 Опубликовано 9 сентября, 2011 · Жалоба Врубился, как сделать. Надо лишь перед БПФ привести шаблон к нулевому среднему и СКО единица. А после вычисления через БПФ быстрой корреляции делить полученные элементы массива значений КК на СКО соответствующих выборок в сигнале. При этом вычисление этих СКО дается практически бесплатно - итерационно. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
getch 0 9 сентября, 2011 Опубликовано 9 сентября, 2011 · Жалоба Подскажите, есть ли возможность из RealFFT(Signal, SignalLen) получить RealFFT(Signal, SignalLen + 1)? Т.е. при приходе нового значения сигнала нужно ли проделывать полностью БПФ или же можно из предыдущего результата БПФ получить текущий? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Alexey Lukin 0 9 сентября, 2011 Опубликовано 9 сентября, 2011 · Жалоба Обычно это делается не по одному отсчёту, а по блокам. Если надо по одному отсчёту — то см. алгоритм Герцеля, он делает это со сложностью N. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
getch 0 9 сентября, 2011 Опубликовано 9 сентября, 2011 · Жалоба Герцеля смотрел ранее, но перед тем, как его применить, провел следующий эксперимент: 1. Сделал RealFFT(Signal, SignalLen), в массиве на выходе посмотрел значение элемента с номером SignalLen. 2. Сделал RealFFT(Signal, SignalLen + 1), в массиве на выходе посмотрел значение элемента также с номером SignalLen. Соответствующие элементы оказались не равны. Т.е. в общем случае получается, что нет никаких совпадений между RealFFT(Signal, SignalLen) и RealFFT(Signal, SignalLen + 1). Это так и должно быть (ведь тогда Герцель неприменим)? Или же это особенность вещественной реализации БПФ и надо переходить на комплексную, чтобы получить совпадение? Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Alexey Lukin 0 9 сентября, 2011 Опубликовано 9 сентября, 2011 · Жалоба Если сигнал поступает по 1 отсчёту, то ни Герцель, ни БПФ вам вообще не нужны. Считайте корреляцию в лоб. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться