Офлафн оценка скорости по дискретным отсчетам, Какие есть методы? Опции

[RHnd]

Исходные данные - есть дискретные квантованные отсчеты некоторого сигнала. Мы знаем, что этот сигнал - выход некоторой механической системы, т.е. как минимум дважды дифференцируем. Более того, если надо, то мы можем априорно предположить некоторые ограничения на первую и вторую производные. Уровень квантования сигнала достаточно маленький по сравнению с диапазоном изменений сигнала. Частота дискретизации достаточно велика по сравнению со скоростью самого процесса. Скажем, дискретизация 1кГц при том, что в составе сигнала нет существенных компонент с частотами выше 50Гц.
Записанный процесс достаточно длинный - несколько сотен тысяч отсчетов, сигнал имеет совершенно произвольную форму (с учетом дифференцируемости), но преобладают низкие частоты.

Задача - как по имеющимся отсчетам оценить скорость? Какие есть оффлайновые методы? Если кто-то сможет дать ссылку на какой-то перечень или обзор, то было бы здорово. Рабочая среда - матлаб.

Я бы решал задачу так - прорядил бы отсчеты, взяв, например, каждый 20-й. Полученные отсчеты описал бы каким-то сплайном, а затем по найденному сплайну вычисли бы значения производной во всех исходных точках. Плюс - в матлабе это делается в две строчки, и если исходный сигнал достаточно гладкий, то получается приемлемая точность. Но при таком методе есть недостаток - выпадают промежуточные, прореженные точки, т.е. часть информации пренебрегается. Кто как решает, на какие методы стоит посмотреть?

[thermit]

Это была иллюстрация. Ваши поправки требуют изменений:

Код

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

[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)]);

Мне упраженения показывают, что второй лучше дает результаты, нет колебаний в начале и конце отрезка, но я от ЦОС далек и не знаю.

UPD: Есть, правда, один недостаток. Если теперь в примере убрать синусоидальные компоненты и оставить только ускоренное движение (s=pos'+0*s;
ds=vel'+0*ds), то результрующая оценка скорости будет сильно колебаться вокруг истинного знаечния. Предполагаю, что это как раз из-за квантования. Как быть в этом случае?

Но Вы же не можете знать, как на самом деле ведет себя система. Может, скорость не постоянная, и там, где Вы видите ступеньку, а думаете, что там прямая, там ступенька и есть на самом деле.

В общем случае это совершенно справедливое замечание.
В тех конкретных примерах, которые я держу в голове и моделирую, речь идет именно о невыскоих скоростях. Это можно косвенно оценить, например, по физическому смыслу эксперимента и по другим наблюдаемым состояниям. При оценке скорости я хотел задать это как ограничение на ускорение, что, в принципе, срабатывает при использовании сглаживающего сплайна. При использовании фильтров такое ограничение, видимо, должно регулироваться через настройку полосы пропускания.