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

28

Я довольно новичок в 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 за добавление моих изображений в мой пост :)

Майкл М
источник
Добро пожаловать в DSP! Является ли красная кривая на втором изображении «сглаженной» версией исходных (зеленых) данных?
Фонон
Красная кривая - это (надеюсь!) Кривая скорости, сгенерированная из fft с последующей фильтрацией, за которой следует омега-арифметика (деленная на 2 * pi f j), за которой следует inv. FFT
Майкл М
1
Возможно, если вы включите более точное математическое выражение или псевдокод, то, что вы сделали, немного прояснит ситуацию.
Фонон
Добавил некоторые теперь, вот общее чувство кода ..
Майкл М
1
Мой вопрос будет: что вы ожидаете увидеть в данных? Вы не узнаете, есть ли у вас хороший подход, если вы не знаете что-то о базовом сигнале, который вы ожидаете увидеть после фильтрации. Кроме того, код, который вы показали, сбивает с толку. Хотя вы не показываете инициализацию fzмассива, похоже, что вы применяете фильтр верхних частот.
Джейсон Р

Ответы:

13

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

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

Если является сигналом с помехами, а является сигналом, который должен быть оценен, минимизируемая функция имеет вид , где - оператор конечных разностей.yxμxy2+Dx1D

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

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

введите описание изображения здесь введите описание изображения здесь введите описание изображения здесь

Даниэль Р. Пипа
источник
Очень нравится этот ответ, собираюсь попробовать. Извините, что я так долго отвечал!
Майкл М,
Отличный ответ. Спасибо за подробности. Я ищу версию C этого кода. Кто-нибудь здесь перенес этот код Matlab на C, которым они хотели бы поделиться? Спасибо.
Рене Лимбергер
Что значит кусочная константа?
Тилапримера
6

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

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

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

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

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

Даниэль Р. Пипа
источник
1
Я бы не согласился с вашим первым заявлением. Вы предполагаете, что интересующий сигнал покрывает всю полосу пропускания входной последовательности, что маловероятно. В этом случае все еще возможно получить улучшенное отношение сигнал / шум, используя линейную фильтрацию, исключая внеполосный шум. Если сигнал сильно передискретизирован, то при таком простом подходе вы можете получить значительное улучшение.
Джейсон Р
Это правда, и это достигается с помощью фильтра Винера, когда вы знаете статистику вашего сигнала и вашего шума.
Даниэль Р. Пипа
Хотя теория шумоподавления вейвлетов сложна, реализация так же проста, как и описанный вами подход. Он включает только банки фильтров и пороговые значения.
Даниэль Р. Пипа
1
Я исследую это сейчас, опубликую мой прогресс выше, спасибо вам и Phonon за всю вашу помощь!
Майкл М
@DanielPipa У меня нет доступа к рассматриваемым пакетам matlab. Можете ли вы предоставить документ или другую ссылку, которая описывает метод, который соответствует вашему коду Matlab.
Джон Робертсон
0

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

Здесь я изменил его код 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)
ryanjdillon
источник