|
|
  |
Офлафн оценка скорости по дискретным отсчетам, Какие есть методы? |
|
|
|
Aug 17 2013, 10:43
|
Знающий
   
Группа: Свой
Сообщений: 518
Регистрация: 12-04-07
Из: Санкт-Петербург
Пользователь №: 26 997

|
1) Если в самом сигнале шумов нет и все проблемы связаны только с квантованием, то прореживание, как мне кажется, даст эффект, близкий к сглаживанию фнч, но без смещения. А потом по сглаженному сигналу оценить производную - такая была исходная идея.
2) Да, действительно, на моих модельных данных filtfilt БИХ ФНЧ + первая разность дают приличный результат. Но, все же, чуть хуже и менее гладкий, чем сглаживающий сплайн с последующим дифференцированием. Возможно, дело в том, что хорошо подобрать один сглаживающий параметр проще, чем ФНЧ. Могу выложить исходные данные.
2) Если бы речь шла про онлайн, то есть эффективные методы, практически не уступающие по фильтрации фнч+дифференциатор, но почти не дающие задержки - тот же slliding differentiator. Задача в том, что оффлайн обработка должна дать лучше результат, чем то, что можно получить онлайн.
|
|
|
|
|
Aug 18 2013, 19:52
|
Знающий
   
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730

|
Вот скрипт иллюстрирующий оценку производной при помощи fir-дифферециатора. Код clear all;
NB=12;%Num bits M=20;%Num sins N=4000;%Seq length L=300;%Differentiators length Fd=1000;%Sam[ling freq Fmax=50;%Max freq
%Random phase p=randn(1,M);
f=(1:M)*Fmax/M; x=(0:N-1).'*2*pi*f/Fd; x=x+repmat(p,length(x),1);
s=sum(sin(x).')/M; %Sequence %s=s+randn(1,length(s))*0.001*std(s);
s=round(s*(2^NB))/(2^NB);% Quantization
ds=(1/M:1/M:1)*cos(x).';%Reference Derivative Sequence
h=firpm(L,[0 Fmax 1.6*Fmax Fd/2]/(Fd/2),[0 M 0 0],'differentiator');
plot(0:Fd/20000:Fd/2-Fd/20000,20*log10(abs(freqz(h,1,10000)))); grid on;
y=filter(h,1,[ s(end-L/2+1:end) s s(1:L/2+1)]);% Calculation Derivative from sequence y=y(L+1:end-1); disp(['Max error = ',num2str(max(abs(ds-y))/std(ds))]);%Compare with ref Точность оценки может быть 10^-6 и круче...
|
|
|
|
|
Aug 19 2013, 06:17
|
Знающий
   
Группа: Свой
Сообщений: 518
Регистрация: 12-04-07
Из: Санкт-Петербург
Пользователь №: 26 997

|
Спасибо за код. Смотрите, убираем вашу строчку Код s=round(s*(2^NB))/(2^NB);% Quantization а после вычисления ds добавляем: Код t=(0:N-1)/Fd; a=[2*ones(1e3,1); zeros(1e3,1); -1.95*ones(1e3,1); 2*zeros(1e3,1)]; vel=lsim(tf(1,[1 0]),a,t); pos=lsim(tf(1,[1 0]),vel,t);
s=pos'+s; ds=vel'+ds;
s=quant(s,range(s)/3e3); И результаты заметно ухудшаются. Как повысить точность оценки? Я, кстати, возможно погорячился на счет 50 Гц. Можно, наверное, 15 Гц ограничить.
|
|
|
|
|
Aug 19 2013, 12:09
|
Знающий
   
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730

|
Это была иллюстрация. Ваши поправки требуют изменений: Код clear all;
NB=12;%Num bits M=20;%Num sins N=4000;%Seq length L=200;%Differentiators length Fd=1000;%Sam[ling freq Fmax=50;%Max freq
%Random phase p=randn(1,M);
f=(1:M)*Fmax/M; x=(0:N-1).'*2*pi*f/Fd; x=x+repmat(p,length(x),1);
s=sum(sin(x).')/M; %Sequence %s=s+randn(1,length(s))*0.001*std(s);
%s=round(s*(2^NB))/(2^NB);% Quantization
ds=(2*pi*f/M)*cos(x).';%Reference Derivative Sequence
t=(0:N-1)/Fd; a=[2*ones(1e3,1); zeros(1e3,1); -1.95*ones(1e3,1); 2*zeros(1e3,1)]; vel=lsim(tf(1,[1 0]),a,t); pos=lsim(tf(1,[1 0]),vel,t);
s=pos'+0.2*s; ds=vel'+0.2*ds;
s=quant(s,range(s)/3e3);
h=firpm(L,[0 Fmax 1.5*Fmax Fd/2]/(Fd/2),[0 2*pi*Fmax 0 0],'differentiator');
plot(0:Fd/20000:Fd/2-Fd/20000,20*log10(abs(freqz(h,1,10000)))); grid on;
y=filter(h,1,[ s(end-L+1:end) s s(1:L)]);% Calculation Derivative from sequence y=y(L+1+L/2:end-L/2); disp(['Mean error = ',num2str(std((ds(L/2:end-L/2+1)-y(L/2:end-L/2+1)))/std((ds)))]);%Compare with ref
|
|
|
|
|
Aug 19 2013, 13:14
|
Гуру
     
Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883

|
Цитата(RHnd @ Aug 16 2013, 02:09)  Дополню: Допустим, у нас некоторая небольшая постоянная скорость. Исходный сигнал представляет из себя, соответственно. прямую. Однако из-за квантования он выглядит как пять точек на одном уровне, ступенька, пять точек на следующему уровне, ступенька, и так далее. Если все точки использовать для сглаживающего сплайна, то он будет стремиться минимизировать расстояние до всех точек и аппроксимирующая кривая будет так же похожа на ступеньки, только сглаженные. В результате оценка производной по такой кривой весьма паршивая. Если же прорядить исходный сигнал и построить сглаживающий сплайн по каждому k-ому отсчету, то это затруднение исчезает при хорошем выборе k. Но Вы же не можете знать, как на самом деле ведет себя система. Может, скорость не постоянная, и там, где Вы видите ступеньку, а думаете, что там прямая, там ступенька и есть на самом деле.
|
|
|
|
|
Aug 19 2013, 14:34
|
Знающий
   
Группа: Свой
Сообщений: 518
Регистрация: 12-04-07
Из: Санкт-Петербург
Пользователь №: 26 997

|
Цитата(thermit @ Aug 19 2013, 16:09)  Это была иллюстрация. Ваши поправки требуют изменений: Спасибо! У меня пара вопросов осталось, если не трудно: 1) Из каких соображений L=300 заменили на L=200? И отсечение с 1.6*Fmax на 1.5*Fmax? 2) Почему Код y=filter(h1,1,[ s(end-L+1:end) s s(1:L)]); а не, например, Код y2=filter(h1,1,[ s(1:L) s s(end-L+1:end)]); Мне упраженения показывают, что второй лучше дает результаты, нет колебаний в начале и конце отрезка, но я от ЦОС далек и не знаю. UPD: Есть, правда, один недостаток. Если теперь в примере убрать синусоидальные компоненты и оставить только ускоренное движение (s=pos'+0*s; ds=vel'+0*ds), то результрующая оценка скорости будет сильно колебаться вокруг истинного знаечния. Предполагаю, что это как раз из-за квантования. Как быть в этом случае? Цитата(Tanya @ Aug 19 2013, 17:14)  Но Вы же не можете знать, как на самом деле ведет себя система. Может, скорость не постоянная, и там, где Вы видите ступеньку, а думаете, что там прямая, там ступенька и есть на самом деле. В общем случае это совершенно справедливое замечание. В тех конкретных примерах, которые я держу в голове и моделирую, речь идет именно о невыскоих скоростях. Это можно косвенно оценить, например, по физическому смыслу эксперимента и по другим наблюдаемым состояниям. При оценке скорости я хотел задать это как ограничение на ускорение, что, в принципе, срабатывает при использовании сглаживающего сплайна. При использовании фильтров такое ограничение, видимо, должно регулироваться через настройку полосы пропускания.
|
|
|
|
|
Aug 19 2013, 15:59
|
Гуру
     
Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883

|
Цитата(RHnd @ Aug 19 2013, 18:34)  В общем случае это совершенно справедливое замечание. В тех конкретных примерах, которые я держу в голове и моделирую, речь идет именно о невыскоих скоростях. Это можно косвенно оценить, например, по физическому смыслу эксперимента и по другим наблюдаемым состояниям. При оценке скорости я хотел задать это как ограничение на ускорение, что, в принципе, срабатывает при использовании сглаживающего сплайна. При использовании фильтров такое ограничение, видимо, должно регулироваться через настройку полосы пропускания. Если сделать ограничения на ускорение и скорость, то такой алгоритм будет нелинейным... А все Ваши фильтры... И мелкие ступеньки могут скрывать все, что угодно. Их никаким образом невозможно "правильно приукрасить"
|
|
|
|
|
Aug 20 2013, 12:55
|
Знающий
   
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730

|
Цитата RHnd: 1) Из каких соображений L=300 заменили на L=200? И отсечение с 1.6*Fmax на 1.5*Fmax? Можно чуть повысить точность меняя х-ки дифференциатора. Но вообще-то - ловля блох... Цитата 2) Почему Код y=filter(h1,1,[ s(end-L+1:end) s s(1:L)]);
а не, например, Код y2=filter(h1,1,[ s(1:L) s s(end-L+1:end)]); Непринципиально. В любом случае в начале и в конце отсчеты числом равным половине длины дифференциатора будут кривыми. Цитата Есть, правда, один недостаток. Если теперь в примере убрать синусоидальные компоненты и оставить только ускоренное движение (s=pos'+0*s; ds=vel'+0*ds), то результрующая оценка скорости будет сильно колебаться вокруг истинного знаечния. Предполагаю, что это как раз из-за квантования. Как быть в этом случае? ФНЧ, ограничивающий полосу полезного сигнала например Код [b,a]=butter(6,0.25*Fmax/Fd); y=filtfilt(b,a,y);
|
|
|
|
|
Aug 20 2013, 14:06
|
Знающий
   
Группа: Свой
Сообщений: 518
Регистрация: 12-04-07
Из: Санкт-Петербург
Пользователь №: 26 997

|
Но если взять такой фнч, то при наличии синусоидальных компонент, он их съест. На данный момент я пришел к выводу, что при налчии среднечастотных колебаний стоит использовать один метод (дифференциатор со сдвигом), а если таких колебаний нет и появляются ступеньки из-за квантования, то другой метод (предварительный фнч или сгл. сплайны). К сожалению, пока у меня не получилось сформирвать метод, который хорошо работал бы и в том случае, и в другом. Жаль. Большое спасибо за помощь! Цитата(thermit @ Aug 20 2013, 16:55)  Непринципиально. В любом случае в начале и в конце отсчеты числом равным половине длины дифференциатора будут кривыми. Да, но степень кривизны получается разная. Возможно, еще лучше поставить x(1)*ones(1,L) - тогда хотя бы разрыва фильтруемого сигнала не будет. А вообще, конечно, к недостаткам метода надо отнести, что L точек теряется. На больших данных это фигня, конечно, но все же.
|
|
|
|
|
Aug 20 2013, 14:59
|
Гуру
     
Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883

|
Цитата(RHnd @ Aug 19 2013, 21:17)  Tanya, простите, а в чем смысл вашего коментария? Вы считаете, что методы (сглаживающий сплайн, фильтр) не позволяют ограничить ускорение и/или высокочастотные компоненты аппроксимирующей кривой? Или считаете, что неправильно их ограничивать основываясь на физическом смысле процесса? Это я к тому, что вот, допустим, имеется длинная ступенька в один квант расстояния. Никакими сглаживаниями нельзя восстановить "правильную скорость". Ограничение на скорость даст только небольшое уменьшение пространства. У Вас ведь реальная система подчиняется "правилам", поэтому все ступеньки уже должны быть "правильными". Вы же, насколько я понимаю Вас, не считаете измерения зашумленными, когда всякие там методы фильтрации были бы уместны... Через Ваши ступеньки можно провести сколько угодно сколь угодно гладких кривых, удовлетворяющих "правилам" - с ограниченной первой и второй производными. Почти всегда.
|
|
|
|
|
Aug 20 2013, 16:20
|
Знающий
   
Группа: Свой
Сообщений: 518
Регистрация: 12-04-07
Из: Санкт-Петербург
Пользователь №: 26 997

|
Я не пойму, к чему вы ведете. Вы считаете, что рассмотренные методы по каким-то причинам нельзя использовать? Или хотите предложить более удачные решения?
Да, существует множество кривых с ограниченными производными, которые после квантования дадут одинаковый набор измерений, это очевидно. Если я вижу набор ступенек, то буду искать движение с малым интегралом квадрата ускорения (ну или с отстутсвующими среднечастотными составляющими). Логика простая - проще предположить отсутствие колебаний, чем такие специальные колебания переменной частоты, что их не видно в измерениях. В зависимости от штрафа за ускорение или полосы пропускания я, естественно, буду получать разные кривые, но они будут близки друг к другу. Если же я увижу осцилляции в измерениях, то буду искать движение, в котором есть среднечастотные компоненты.
Придумалась такая возможная идея. Параметризуем алгоритм аппроксимации некоторым параметром Q . Это может быть или штраф за ускорение, или полоса пропускания фнч. Чем Q меньше, тем более сглаженный и низкочастотный сигнал получается на выходе. Далее этот выходной сигнал квантуем с известным шагом квантования. И считаем отклонение исходного сигнала от проквантованного сглаженного. Далее два варианта: либо найти наименший Q, при котором имеем нулевой отклонение, либо минимизировать некоторый критерий, сочетающий норму отклонения и штраф за величину Q. Т.е. найти самую низкочастотную кривую, которая после квантования будет близка к исходному сигналу. UPD: Действительно, я сейчас промоделировал - работает, по первым прикидкам, очень неплохо. И для сигналов с среднечастотными колебательными компонентами, и без таких - один и тот же код выдает очень хорошую оценку скорости. Уф.
|
|
|
|
Guest_TSerg_*
|
Aug 21 2013, 06:57
|
Guests

|
Это не задача ЦОС, это задача идентификации динамической системы и построения оптимального наблюдателя за вектором состояния системы.
|
|
|
|
|
  |
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|