Выбор правильного фильтра для данных акселерометра

Я довольно новичок в DSP, и провел некоторые исследования возможных фильтров для сглаживания данных акселерометра в Python. Пример типа данных, которые я получу, можно увидеть на следующем рисунке:

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

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

БПФ дал некоторые интересные результаты. Одна из моих попыток состояла в том, чтобы БПФ ускорить сигнал, а затем визуализировать низкие частоты, чтобы получить абсолютное значение БПФ 0. Затем я использовал омега-арифметику и обратное БПФ, чтобы получить график скорости. Результаты были следующими:

Это хороший способ идти о вещах? Я пытаюсь устранить общую шумовую природу сигнала, но необходимо определить очевидные пики, например, около 80 секунд.

Я также устал использовать фильтр низких частот для исходных данных акселерометра, который отлично справился со сглаживанием, но я не совсем уверен, куда идти дальше. Любое руководство о том, куда идти отсюда, будет очень полезно!

РЕДАКТИРОВАТЬ: немного кода:

for i in range(len(fz)): 
    testing = (abs(Sz[i]))/Nz

    if fz[i] < 0.05:
        Sz[i]=0

Velfreq = []
Velfreqa = array(Velfreq)
Velfreqa = Sz/(2*pi*fz*1j)
Veltimed = ifft(Velfreqa)
real = Veltimed.real

Итак, по сути, я выполнил БПФ на моих данных акселерометра, дав Sz, отфильтровывающий высокие частоты с помощью простого фильтра кирпичной стены (я знаю, что это не идеально). Затем я использую арифметику омега на БПФ данных. Также большое спасибо datageist за добавление моих изображений в мой пост :)

Как отметил @JohnRobertson в « Сумке хитростей для шумоподавления сигналов при поддержании резких переходов» , шумоподавление Total Variaton (TV) является еще одной хорошей альтернативой, если ваш сигнал является кусочно-постоянным. Это может иметь место для данных акселерометра, если ваш сигнал продолжает изменяться между различными плато.

Ниже приведен код Matlab, который выполняет шумоподавление ТВ в таком сигнале. Код основан на статье «Метод расширенного лагранжиана для полного восстановления видео с вариациями» . Параметры и должны быть отрегулированы в соответствии с уровнем шума и характеристиками сигнала.$\mu$$\rho$

Если является сигналом с помехами, а является сигналом, который должен быть оценен, минимизируемая функция имеет вид , где - оператор конечных разностей.$y$$x$$\mu\|{x-y}\|^2+\|{Dx}\|_1$$D$

function denoise()

f = [-1*ones(1000,1);3*ones(100,1);1*ones(500,1);-2*ones(800,1);0*ones(900,1)];
plot(f);
axis([1 length(f) -4 4]);
title('Original');
g = f + .25*randn(length(f),1);
figure;
plot(g,'r');
title('Noisy');
axis([1 length(f) -4 4]);
fc = denoisetv(g,.5);
figure;
plot(fc,'g');
title('De-noised');
axis([1 length(f) -4 4]);

function f = denoisetv(g,mu)
I = length(g);
u = zeros(I,1);
y = zeros(I,1);
rho = 10;

eigD = abs(fftn([-1;1],[I 1])).^2;
for k=1:100
    f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
    v = D(f)+(1/rho)*y;
    u = max(abs(v)-1/rho,0).*sign(v);
    y = y - rho*(u-D(f));
end

function y = D(x)
y = [diff(x);x(1)-x(end)];

function y = Dt(x)
y = [x(end)-x(1);-diff(x)];

Полученные результаты:

Проблема в том, что у вашего шума плоский спектр. Если вы предполагаете белый гауссов шум (который оказывается хорошим предположением), его плотность спектра мощности постоянна. Грубо говоря, это означает, что ваш шум содержит все частоты. Вот почему любой частотный подход, например, DFT или фильтры нижних частот, не является хорошим. Какими будут ваши частоты среза, поскольку ваш шум распространяется по всему спектру?

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

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

В MATLAB введите «wavemenu», а затем «SWT denoising 1-D». Затем «Файл», «Пример анализа», «Шумовые сигналы», «с Хааром на уровне 5, шумные блоки». В этом примере используется вейвлет Хаара, который должен нормально работать для вашей задачи.

Я не очень хорош в Python, но я считаю, что вы можете найти несколько пакетов NumPy, которые выполняют шумоподавление вейвлетом Хаара.

По предложению Даниэля Пипа я взглянул на шум вейвлета и нашел эту прекрасную статью Франсиско Бланко-Сильвы.

Здесь я изменил его код Python для обработки изображений, чтобы работать с 2D (акселерометр), а не с 3D (изображение) данными.

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

def wavelet_denoise(data, wavelet, noise_sigma):
    '''Filter accelerometer data using wavelet denoising

    Modification of F. Blanco-Silva's code at: https://goo.gl/gOQwy5
    '''
    import numpy
    import scipy
    import pywt

    wavelet = pywt.Wavelet(wavelet)
    levels  = min(15, (numpy.floor(numpy.log2(data.shape[0]))).astype(int))

    # Francisco's code used wavedec2 for image data
    wavelet_coeffs = pywt.wavedec(data, wavelet, level=levels)
    threshold = noise_sigma*numpy.sqrt(2*numpy.log2(data.size))

    new_wavelet_coeffs = map(lambda x: pywt.threshold(x, threshold, mode='soft'),
                             wavelet_coeffs)

    return pywt.waverec(list(new_wavelet_coeffs), wavelet)

Где:

• wavelet- имя строки используемой формы вейвлета (см. pywt.wavelist(), например 'haar')
• noise_sigma - стандартное отклонение шума от данных
• data - массив значений для фильтрации (например, данные оси x, y или z)