Анализ зашумленных данных

Вот мое решение с использованием фильтра Калмана. Вам нужно будет настроить параметры (даже +- порядков), если вы хотите сгладить более или менее.

#!/usr/bin/env octave

% Kalman filter to smooth measures of altitude and estimate
% speed and acceleration. The continuous time model is more or less as follows:
% derivative of altitude := speed
% derivative of speed := acceleration
% acceleration is a Wiener process

%------------------------------------------------------------
% Discretization of the continuous-time linear system
% 
%   d  |x|   | 0 1 0 | |x|
%  --- |v| = | 0 0 1 | |v|   + "noise"
%   dt |a|   | 0 0 0 | |a|
%
%   y = [1 0 0] |x|     + "measurement noise"
%               |v|
%               |a|
%
st = 0.05;    % Sampling time
A = [1  st st^2/2;
     0  1  st    ;
     0  0  1];
C = [1 0 0];

%------------------------------------------------------------
% Fine-tune these parameters! (in particular qa and R)
% The acceleration follows a "random walk". The greater is the variance qa,
% the more "reactive" the system is expected to be, i.e.
% the more the acceleration is expected to vary
% The greater is R, the more noisy is your measurement instrument
% (less "accuracy" of the barometric altimeter);
% if you increase R, you will smooth the estimate more
qx = 1.0;                      % Variance of model noise for position
qv = 1.0;                      % Variance of model noise for speed
qa = 50.0;                     % Variance of model noise for acceleration
Q  = diag([qx, qv, qa]);
R  = 100.0;                    % Variance of measurement noise
                               % (10^2, if 10ft is the standard deviation)

load data.txt  % Put your measures in this file

est_position     = zeros(length(data), 1);
est_speed        = zeros(length(data), 1);
est_acceleration = zeros(length(data), 1);

%------------------------------------------------------------
% Kalman filter
xhat = [0;0;0];     % Initial estimate
P    = zeros(3,3);  % Initial error variance
for i=1:length(data),
   y = data(i);
   xpred = A*xhat;                                    % Prediction
   Ppred = A*P*A' + Q;                                % Prediction error variance
   Lambdainv = 1/(C*Ppred*C' + R);
   xhat  = xpred + Ppred*C'*Lambdainv*(y - C*xpred);  % Update estimation
   P = Ppred - Ppred*C'*Lambdainv*C*Ppred;            % Update estimation error variance
   est_position(i)     = xhat(1);
   est_speed(i)        = xhat(2);
   est_acceleration(i) = xhat(3);
end

%------------------------------------------------------------
% Plot
figure(1);
hold on;
plot(data, 'k');               % Black: real data
plot(est_position, 'b');       % Blue:  estimated position
plot(est_speed, 'g');          % Green: estimated speed
plot(est_acceleration, 'r');   % Red:   estimated acceleration
pause

Вы можете попробовать пропустить данные через фильтр нижних частот. Это сгладит высокочастотный шум. Может быть, простой FIR.

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

Вы пытались выполнить среднее значение окна прокрутки ваших значений? В основном вы выполняете окно, скажем, 10 значений (от 0 до 9) и вычисляете его среднее значение. затем вы прокручиваете окно на одну точку (от 1 до 10) и пересчитываете. Это сгладит значения, сохраняя количество точек относительно неизменным. Окна большего размера дают более плавные данные за счет потери более высокочастотной информации.

Вы можете использовать медиану вместо среднего, если ваши данные содержат выбросы.

Вы также можете попробовать использовать автокорреляцию.

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

Что касается создания некоторого чувства скорости и ускорения во время полета, это должно быть простым усреднением скорости из нескольких разных результатов. Что-то вроде: EsitimatedV = Vmeasured*(1/n) + (1 - 1/n)*EstimatedV. Установите n в зависимости от того, насколько быстро вы хотите изменить скорость.

Я ничего не знаю о ракетах. Я нарисовал ваши точки, и они выглядят прекрасно.

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

Предположение:

1. Контролируйте максимальную высоту на протяжении всего полета.
2. Непрерывно следите за апогеем, (скажем, просто) сравнивая самые последние несколько точек с текущим максимумом.
3. Пока вы не достигнете максимума, с фиксированным (0,0) и некоторым произвольным набором узлов вычислите набор естественных сплайнов до текущей высоты. Используйте остатки относительно сплайнов, чтобы решить, какие данные отбросить. Пересчитайте сплайны.
4. По максимуму сохраняйте самые последние рассчитанные сплайны. Начните вычислять новый набор сплайнов для кривой за апогеем.

Модель ARIMA и поиск автокорреляции в остатке — стандартная процедура. Модель волатильности другая.