Анализ шумных данных
Я недавно запустил ракету с барометрическим высотомером, который является с точностью до примерно 10 футов (вычислен через данные, полученные во время полета). Зарегистрированные данные находятся в инкрементах времени 0,05 секунд на образец, и график высоты по сравнению со временем в значительной степени походит на него, должен при уменьшении по всему полету.
Проблема состоит в том, когда я пытаюсь вычислить другие значения, такие как скорость или ускорение от данных, точность измерений делает расчетные значения в значительной степени бесполезными. Какие методы я могу использовать для сглаживания данных так, чтобы я мог вычислить (или приблизительный) рыночную стоимость для скорости и ускорения? Важно, чтобы крупные события остались на месте вовремя, прежде всего 0 для для первой записи и самой высокой точки во время рейса (2707).
Высотные данные следуют и измеряются в ft над уровнем земли. Первый раз был бы 0.00, и каждый образец спустя 0.05 секунды после предыдущего образца. Скачок в начале полета происходит из-за технической проблемы, которая произошла во время старта, и удаление скачка оптимально.
Я первоначально пытался использовать линейную интерполяцию, насчитывая соседние точки данных, но потребовалось много повторений для сглаживания данных достаточно для интеграции, и выравнивание кривой удалило важный апогей и события уровня земли.
Вся справка значительно ценится. Обратите внимание, что это не полный набор данных, и я ищу предложения на лучших способах проанализировать данные, не, чтобы кто-то ответил с преобразованным набором данных. Было бы хорошо использовать алгоритм на борту будущих ракет, которые могут предсказать текущую высоту/скорость/ускорение, не зная полные полетные данные, хотя это не требуется.
00000
00000
00000
00076
00229
00095
00057
00038
00048
00057
00057
00076
00086
00095
00105
00114
00124
00133
00152
00152
00171
00190
00200
00219
00229
00248
00267
00277
00286
00305
00334
00343
00363
00363
00382
00382
00401
00420
00440
00459
00469
00488
00517
00527
00546
00565
00585
00613
00633
00652
00671
00691
00710
00729
00759
00778
00798
00817
00837
00856
00885
00904
00924
00944
00963
00983
01002
01022
01041
01061
01080
01100
01120
01139
01149
01169
01179
01198
01218
01238
01257
01277
01297
01317
01327
01346
01356
01376
01396
01415
01425
01445
01465
01475
01495
01515
01525
01545
01554
01574
01594
01614
01614
01634
01654
01664
01674
01694
01714
01724
01734
01754
01764
01774
01794
01804
01814
01834
01844
01854
01874
01884
01894
01914
01924
01934
01954
01954
01975
01995
01995
02015
02015
02035
02045
02055
02075
02075
02096
02096
02116
02126
02136
02146
02156
02167
02177
02187
02197
02207
02217
02227
02237
02237
02258
02268
02278
02278
02298
02298
02319
02319
02319
02339
02349
02359
02359
02370
02380
02380
02400
02400
01914
02319
02420
02482
02523
02461
02502
02543
02564
02595
02625
02666
02707
02646
02605
02605
02584
02574
02543
02543
02543
02543
02543
02543
02554
02543
02554
02554
02554
02554
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02543
02533
02543
02543
02543
02543
02543
02543
02543
02543
02533
02523
02523
02523
02523
02523
02523
02523
02523
02543
02523
02523
02523
02523
02523
02523
02523
02523
02513
02513
02502
02502
02492
02482
02482
02482
02482
02482
02482
02482
02482
02482
02482
02482
02482
02482
02482
02482
02472
02472
02472
02461
02461
02461
02461
02451
02451
02451
02461
02461
02451
02451
02451
02451
02451
02451
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02441
02431
02441
02431
02441
02431
02420
02431
02420
02420
02420
02420
02420
02420
02420
02420
02420
02420
02420
02420
02410
02420
02410
02410
02410
02410
02400
02400
02410
02400
02400
02400
02400
02400
02400
02400
02400
02400
02400
02400
02400
02390
02390
02390
02380
02380
02380
02380
02380
02380
02380
02380
02380
02380
02380
02380
02380
02370
02370
02380
02370
02359
02359
02359
02359
02359
02359
02359
02359
02359
02359
02359
02359
02359
02359
02349
02349
02349
02349
02349
02339
02339
02339
02339
02339
02339
02339
02339
02339
02339
02339
02339
02339
6 ответов
Вот мое решение, использующее фильтр Калмана . Вам нужно будет настроить параметры (даже +- порядки величины), если вы хотите сгладить более или менее.
#!/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
пробовали ли вы выполнить среднее значение окна прокрутки ваших значений? Обычно вы создаете окно, например, из 10 значений (от 0 до 9) и вычисляете его среднее значение. затем прокручиваете окно на одну точку (от 1 до 10) и пересчитываете. Это сгладит значения, сохраняя при этом количество точек относительно неизменным. Большие окна дают более гладкие данные за счет потери более высокочастотной информации.
Вы можете использовать медианное значение вместо среднего, если ваши данные показывают выбросы выбросов.
Вы также можете попробовать с Автокорреляция .
Вы можете попробовать пропустить данные через фильтр нижних частот. Это сгладит высокочастотный шум. Может быть, простой FIR.
Кроме того, вы можете извлечь основные события из необработанных данных, но использовать полиномиальную подгонку для данных скорости и ускорения.
Один из способов, которым вы можете подойти к анализу ваших данных, - это попытаться сопоставить ее с какой-либо моделью, сгенерировать функцию, а затем протестировать ее соответствие вашему набору данных .... Это может быть довольно сложно и, вероятно, не нужно ... но дело в том, что вместо генерации данных об ускорении / скорости непосредственно из ваших данных вы можете сопоставить их с вашей моделью (довольно просто для ракеты - некоторое ускорение вверх с последующим медленным спуском с постоянной скоростью.) По крайней мере, как я бы сделал это в физическом эксперименте.
Что касается создания некоторого ощущения скорости и ускорения во время полета, это должно быть таким же простым усреднением скорость из нескольких разных результатов. Что-то вроде: Оценочная V = V измеренная * (1 / n) + (1 - 1 / n) * Оценка V. Установите n в зависимости от того, насколько быстро вы хотите, чтобы ваша скорость регулировалась.
Я ничего не знаю о ракетах. Я нанес ваши точки, и они выглядят прекрасно.
На основании того, что я вижу на этом графике, позвольте мне предположить, что обычно существует единственный апогей и что функция, которая дала начало вашим точкам, не имеет производной по времени в этом апогее.
Предложение:
- Следите за максимальной высотой на протяжении всего полета.
- Постоянно следите за апогеем, (скажем, просто) сравнивая несколько последних точек с текущим максимумом.
- Пока вы не достигнете максимума, с помощью ( 0,0) фиксированный и произвольный набор узлов вычисляет набор естественных сплайнов до текущей высоты. Используйте остатки по сплайнам, чтобы решить, какие данные отбросить. Пересчитайте шлицы.
- На максимуме сохраните самые последние рассчитанные шлицы. Начните вычислять новый набор сплайнов для кривой за апогеем.
Модель ARIMA и поиск автокорреляции в остатке являются стандартной процедурой. Другая модель волатильности.