为加速度计数据选择正确的滤波器

我对DSP相当陌生，并且已经对可能的过滤器进行了一些研究，以平滑python中的加速度计数据。下图显示了一个疾病类型的示例：

本质上，我正在寻求有关平滑此数据以最终转换的建议它分为速度和位移。我了解到手机的加速度计非常嘈杂。

我认为目前无法使用Kalman滤波器，因为我无法握住设备来参考数据产生的噪声（我读到将设备放平并找到噪声的关键。这些读数会产生多少噪声？）

FFT产生了一些有趣的结果。我的尝试之一是对加速度信号进行FFT，然后渲染低频使其FFT绝对值为0。然后，我使用了ω算术和逆FFT来获得速度图。结果如下：

这是处理事情的好方法吗？我正在尝试消除信号的整体噪声性质，但是需要识别明显的峰值，例如在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

因此，从本质上讲，IVE对我的加速度计数据执行了FFT，从而使用简单的砖墙滤波器将Sz滤波后的高频信号去除了（我知道它并不理想）。然后我在数据的FFT上使用欧米茄算法。
非常感谢datageist将我的图像添加到我的帖子中：）

欢迎使用DSP！第二张图片中的红色曲线是原始（绿色）数据的“平滑”版本吗？

红色曲线是（希望如此！）由fft生成的速度曲线，然后进行滤波，然后是ω算术（除以2 * pifj），然后是inv。 fft

也许，如果您为所执行的操作添加了更精确的数学表达式或伪代码，将会使您的工作变得清楚一些。

现在添加了一些，多数民众赞成在代码的一般感觉..

我的问题是：您希望在数据中看到什么？除非您对滤波后希望看到的基础信号有所了解，否则您将不知道自己是否有一个好的方法。另外，您显示的代码令人困惑。尽管您没有显示fz数组的初始化，但似乎您正在应用高通滤波器。

#1 楼

正如@JohnRobertson在“在保持急剧过渡的同时对信号进行降噪的技巧包”中指出的那样，如果您的信号是分段恒定的，则总可变噪声（TV）降噪是另一个不错的选择。如果您的信号在不同的平稳区域之间保持变化，则加速度计数据可能就是这种情况。该代码基于论文“用于整体变化视频恢复的增强拉格朗日方法”。必须根据噪声水平和信号特性来调整参数$ \ mu $和$ \ rho $。

如果$ y $是噪声信号而$ x $是要估计的信号，要最小化的函数是$ \ mu \ | {xy} \ | ^ 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)];

结果：

$ \ begingroup $
真的很喜欢这个答案，然后继续尝试。抱歉，我花了这么长时间回复！
$ \ endgroup $
– Michael M
2012年8月28日在16:51

$ \ begingroup $
很好的答案。感谢您的详细信息。我正在寻找此代码的C版本。这里有人将此matlab代码移植到他们想共享的C吗？谢谢。
$ \ endgroup $
– Rene Limberger
16年1月12日在18:24

$ \ begingroup $
分段常数是什么意思？
$ \ endgroup $
– tilaprimera
18年6月1日在20:25

#2 楼

问题在于您的噪声频谱平坦。如果假设白色高斯噪声（事实证明是一个很好的假设），则其功率谱密度是恒定的。粗略地说，这意味着您的噪声包含所有频率。这就是为什么任何频率方法，例如DFT或低通滤波器不是一个很好的选择。由于噪声遍布整个频谱，您的截止频率是多少？

这个问题的一个答案是维纳滤波器，该滤波器需要了解噪声和所需信号的统计信息。基本上，噪声信号（信号+噪声）会在预期比您的信号更平滑的频率上衰减，并在预期比您的信号更平滑的频率上放大。但是，我建议使用非线性处理的更现代方法，例如小波去噪。这些方法提供了极好的结果。基本上，噪声信号首先分解为小波，然后将小系数归零。由于小波具有多分辨率特性，因此该方法有效（但DFT无效）。也就是说，在小波变换定义的频带中分别处理信号。

在MATLAB中，键入“ wavemenu”，然后键入“ SWT去噪一维”。然后是“文件”，“示例分析”，“噪声信号”，“ Haar处于第5级噪声块”。此示例使用Haar小波，该方法应该可以很好地解决您的问题。

$ \ begingroup $
我不同意你的第一句话。您假设感兴趣的信号覆盖了输入序列的整个带宽，这不太可能。在这种情况下，仍可以使用线性滤波来获得改善的信噪比，从而消除带外噪声。如果信号高度过采样，那么使用这种简单的方法可能会获得很大的改进。
$ \ endgroup $
–Jason R
2012年8月7日在1:37

$ \ begingroup $
是的，当您知道信号和噪声的统计信息时，这是通过维纳滤波器实现的。
$ \ endgroup $
–丹尼尔·R·皮帕（Daniel R. Pipa）
2012年8月7日上午10:57

$ \ begingroup $
尽管小波去噪背后的理论很复杂，但实现方法却与您描述的方法一样简单。它仅涉及滤波器组和阈值。
$ \ endgroup $
–丹尼尔·R·皮帕（Daniel R. Pipa）
2012年8月7日在12:08

$ \ begingroup $
我现在正在对此进行研究，感谢您和Phonon迄今为止的所有帮助，我将在上面发布我的进展！
$ \ endgroup $
– Michael M
2012年8月7日14:11

$ \ begingroup $
@DanielPipa我无权访问有问题的matlab软件包。您能否提供描述与您的Matlab代码相对应的方法的论文或其他参考资料。
$ \ endgroup $
–约翰·罗伯逊
2012年8月8日在22:44

#3 楼

根据丹尼尔·皮帕（Daniel Pipa）的建议，我看了小波去噪，发现了弗朗西斯科·布兰科·席尔瓦（Francisco Blanco-Silva）的这篇出色文章。

在这里，我修改了他的Python代码用于图像处理，以处理2D（加速度计）而不是3D（图像）数据。

注意，在Francisco的示例中，软阈值的阈值是“弥补”的。考虑这一点并为您的应用程序进行修改。

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轴数据）

编程黑洞网

为加速度计数据选择正确的滤波器

评论

#1 楼

评论

#2 楼

评论

#3 楼