主要内容

信号平滑

这个例子展示了如何使用移动平均滤波器和重采样来隔离时间的周期性成分对每小时温度读数的影响,以及从开环电压测量中去除不必要的线噪声。该示例还展示了如何使用中值滤波器平滑时钟信号的电平,同时保留边缘。该示例还展示了如何使用Hampel过滤器去除较大的离群值。

动机

平滑是我们在数据中发现重要模式的方法,同时忽略不重要的东西(如噪音)。我们使用滤波来实现平滑。平滑的目标是产生缓慢的价值变化,以便更容易看到数据中的趋势。

有时当你检查输入数据时,你可能希望平滑数据以看到信号的趋势。在我们的例子中,我们有一组温度读数,以摄氏度为单位,每小时在波士顿的洛根机场,2011年1月的整个月。

负载bostemp天=一句子* 24)/ 24;情节(天,tempC)轴ylabel (“临时(\ circC)”)包含(时间自2011年1月1日起(天))标题(“洛根机场干球温度(来源:美国国家海洋和大气管理局)”

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含一个类型为线的对象。

注意,我们可以直观地看到一天中的时间对温度读数的影响。如果你只对一个月里每天的温度变化感兴趣,那么每小时的波动只会造成噪音,使每天的变化难以识别。为了去除一天中时间的影响,我们现在想要使用移动平均滤波器平滑我们的数据。

移动平均滤波器

在其最简单的形式中,长度为N的移动平均滤波器取波形每个连续N个样本的平均值。

为了对每个数据点应用移动平均滤波器,我们构造滤波器的系数,使每个点的权重相等,并对总平均值贡献1/24。这就得到了每24小时内的平均温度。

hoursPerDay = 24;coeff24hMA = ones(1,小时sperday)/小时sperday;avg24hTempC = filter(coeff24hMA, 1, tempC);情节(天,[tempC avg24hTempC])传说(每小时的临时的《平均24小时》(延迟播出)“位置”“最佳”) ylabel (“临时(\ circC)”)包含(时间自2011年1月1日起(天))标题(“洛根机场干球温度(来源:美国国家海洋和大气管理局)”

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含2个类型为线的对象。这些对象表示小时温度,24小时平均(延迟)。

过滤器延迟

注意,过滤后的输出延迟了大约12个小时。这是由于移动平均滤波器有延迟。

任何长度为N的对称滤波器都有(N-1)/2个样本的延迟。我们可以手工解释这次延误。

fDelay =(长度(coeff24hMA) 1) / 2;情节(天,tempC,...avg24hTempC days-fDelay / 24日)轴传奇(每小时的临时的24小时平均的“位置”“最佳”) ylabel (“临时(\ circC)”)包含(时间自2011年1月1日起(天))标题(“洛根机场干球温度(来源:美国国家海洋和大气管理局)”

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含2个类型为线的对象。这些对象表示小时温度,24小时平均。

提取的平均差异

另外,我们还可以使用移动平均滤波器来更好地估计一天中的时间如何影响整体温度。要做到这一点,首先,从每小时的温度测量中减去平滑数据。然后,将不同的数据划分为天数,并取当月所有31天的平均值。

图deltaTempC = tempC - avg24hTempC;deltaTempC =重塑(deltaTempC, 24, 31).';情节(桥,意味着(deltaTempC))轴标题(“平均气温与24小时平均气温之差”)包含(“一小时(从午夜开始)”) ylabel (“温差(\ circC)”

图中包含一个轴对象。标题为“平均温度差与24小时平均值”的轴对象包含一个类型为“线”的对象。

提取包络峰值

有时我们也想对我们的温度信号每天的高低变化有一个平稳变化的估计。为此,我们可以使用信封功能连接在24小时期间的子集检测到的极端高和低。在本例中,我们确保每个极端高和极端低之间至少有16个小时。我们还可以通过取两个极端之间的平均值来了解高点和低点的趋势。

[envHigh, envLow] = envelope(tempC,16, 16)“高峰”);envMean = (envHigh + envLow) / 2;情节(天,tempC,...天,envHigh,...天,envMean,...天,envLow)轴传奇(每小时的临时的“高”“的意思是”“低”“位置”“最佳”) ylabel (“临时(\ circC)”)包含(时间自2011年1月1日起(天))标题(“洛根机场干球温度(来源:美国国家海洋和大气管理局)”

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含4个类型为线的对象。这些对象代表每小时温度,高,平均,低。

加权移动平均滤波器

其他类型的移动平均过滤器并不对每个样本平均加权。

另一个常见的过滤器遵循的二项展开 1 / 2 1 / 2 n .这种类型的滤波器在n值较大时近似于一条正态曲线。它对滤除高频噪声很有用,当n值较小时 1 / 2 1 / 2 然后将输出进行迭代卷积 1 / 2 1 / 2 规定的次数在这个例子中,总共使用五个迭代。

H = [1/2 /2];binomialCoeff = conv (h, h);n = 1:4 binomialeff = conv(binomialeff,h);结束figure fDelay = (length(binomialeff)-1)/2;binomialMA = filter(binomialeff, 1, tempC); / /输出情节(天,tempC,...binomialMA days-fDelay / 24日)轴传奇(每小时的临时的“二项加权平均”“位置”“最佳”) ylabel (“临时(\ circC)”)包含(时间自2011年1月1日起(天))标题(“洛根机场干球温度(来源:美国国家海洋和大气管理局)”

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含2个类型为线的对象。这些对象表示每小时温度,二项式加权平均。

另一个类似于高斯膨胀滤波器的滤波器是指数移动平均滤波器。这种类型的加权移动平均滤波器易于构造,不需要大的窗口大小。

你用0到1之间的参数来调整指数加权移动平均滤波器。alpha值越高,平滑度越低。

α= 0.45;exponentialMA = filter(alpha, [1 alpha-1], tempC);情节(天,tempC,...binomialMA days-fDelay / 24日,...days-1/24 exponentialMA)轴传奇(每小时的临时的...“二项加权平均”...指数加权平均的“位置”“最佳”) ylabel (“临时(\ circC)”)包含(时间自2011年1月1日起(天))标题(“洛根机场干球温度(来源:美国国家海洋和大气管理局)”

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含3个类型的对象线。这些对象代表每小时温度,二项式加权平均,指数加权平均。

放大一天的读数。

轴([3 4 -5 2])

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含3个类型的对象线。这些对象代表每小时温度,二项式加权平均,指数加权平均。

Savitzky-Golay过滤器

您会注意到,通过平滑数据,极端值在某种程度上被剪掉了。

为了更接近地跟踪信号,您可以使用加权移动平均滤波器,它试图以最小二乘的方式,在特定数量的样本上拟合一个特定顺序的多项式。

为了方便,您可以使用该函数sgolayfilt实现一个Savitzky-Golay平滑滤波器。使用sgolayfilt,指定数据的奇数长度段和严格小于段长度的多项式阶。的sgolayfilt函数在内部计算平滑多项式系数,执行延迟对齐,并在数据记录的开始和结束处理瞬态效应。

cubicMA = sgolayfilt(tempC, 3,7);quarticMA = sgolayfilt(tempC, 4,7); / /调整时间quticma = sgolayfilt(tempC, 5, 9);plot(天,[tempC cubicMA quarticMA quinticMA]每小时的临时的“马Cubic-Weighted”“马Quartic-Weighted”...“马Quintic-Weighted”“位置”“东南”) ylabel (“临时(\ circC)”)包含(时间自2011年1月1日起(天))标题(“洛根机场干球温度(来源:美国国家海洋和大气管理局)”轴([3 5 -5 2])

图中包含一个轴对象。轴对象的标题为Logan机场干球温度(来源:NOAA)包含4个类型为线的对象。这些对象代表每小时温度,三次加权MA,四次加权MA,五次加权MA。

重采样

有时为了正确地应用移动平均线而对信号进行重新采样是有益的。

在我们的下一个例子中,我们在60hz交流电源线噪声干扰的情况下采样了模拟仪器输入端的开环电压。我们以1khz的采样率对电压进行采样。

负载openloop60hertzfs = 1000;t = (0:numel(openLoopVoltage)-1) / fs;情节(t, openLoopVoltage) ylabel (“电压(V)”)包含(“时间(s)”)标题(开环电压测量的

图中包含一个轴对象。标题为“开环电压测量”的轴对象包含一个类型为line的对象。

让我们尝试使用移动平均滤波器来去除线噪声的影响。

如果你构造一个均匀加权的移动平均滤波器,它将移除任何与滤波器持续时间相关的周期性成分。

在1000hz采样时,一个60hz的完整周期中大约有1000 / 60 = 16.667个样本。让我们尝试“集中”并使用17点过滤器。这将给我们在1000hz / 17 = 58.82 Hz的基本频率的最大滤波。

情节(t, sgolayfilt (openLoopVoltage 1 17) ylabel (“电压(V)”)包含(“时间(s)”)标题(开环电压测量的)传说(“58.82 Hz的移动平均滤波器”...“位置”“东南”

图中包含一个轴对象。标题为“开环电压测量”的轴对象包含一个类型为line的对象。这个对象表示58.82 Hz的移动平均滤波器。

注意,虽然电压被显著地平滑,它仍然包含一个小的60赫兹纹波。

我们可以显著地减少纹波,如果我们对信号进行重新采样,这样我们就可以通过移动平均滤波器捕获60 Hz信号的一个完整周期。

如果我们在17 * 60hz = 1020hz重新采样信号,我们可以使用17点移动平均滤波器来去除60hz的线噪声。

fsResamp = 1020;vResamp = reample (openLoopVoltage, fsResamp, fs); / /重新采样tResamp = (0:numel(vResamp)-1) / fsResamp; / /重删vAvgResamp = sgolayfilt (vResamp 1 17);情节(tResamp vAvgResamp) ylabel (“电压(V)”)包含(“时间(s)”)标题(开环电压测量的)传说(“60赫兹移动平均滤波器”...“位置”“东南”

图中包含一个轴对象。标题为“开环电压测量”的轴对象包含一个类型为line的对象。这个对象表示运行在60hz的移动平均滤波器。

中值滤波器

移动平均,加权移动平均,和Savitzky-Golay过滤器平滑所有他们过滤的数据。然而,这可能并不总是我们想要的。例如,如果我们的数据是从时钟信号中获取的,并且有我们不希望平滑的锐边,该怎么办?到目前为止讨论的过滤器并不是很好:

负载clockexyMovingAverage = conv (x)的(5 - 1)/ 5,“相同”);ySavitzkyGolay = sgolayfilt (x, 3、5);情节(t x,...t yMovingAverage...t, ySavitzkyGolay)传说(原始信号的“移动平均”“Savitzky-Golay”

图中包含一个轴对象。轴对象包含3个类型为line的对象。这些物体代表原始信号,移动平均,萨维茨基-戈莱。

移动平均和Savitzky-Golay分别滤除时钟信号边缘附近的欠校正和过校正。

保持边缘的简单方法是使用中值过滤器:

yMedFilt = medfilt1 (x 5“截断”);情节(t x,...t, yMedFilt)传说(原始信号的“中值滤波”

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表原始信号,中值滤波。

通过Hampel过滤器去除离群值

许多过滤器对异常值很敏感。与中值滤波器密切相关的一个滤波器是汉普尔滤波器。这个过滤器有助于在不过度平滑数据的情况下去除信号中的异常值。

为了看到这一点,加载一段火车哨声的录音,并添加一些人工噪声峰值:

负载火车y(1:400:结束)= 2.1;情节(y)

图中包含一个轴对象。axis对象包含一个类型为line的对象。

因为我们引入的每个尖峰都只有一个样本的持续时间,所以我们可以使用一个只有三个元素的中值过滤器来删除尖峰。

持有情节(medfilt1 (y, 3))传奇(原始信号的“过滤信号”

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表原始信号,滤波后的信号。

滤波器去除了尖波,但同时也去除了原始信号的大量数据点。Hampel过滤器的工作原理类似于中值过滤器,但是它只替换那些与本地中值相差几个标准差的值。

传说hampel (y, 13) (“位置”“最佳”

图中包含一个轴对象。轴对象包含3个类型为line的对象。这些对象代表原始信号,滤波信号,离群值。

只有异常值从原始信号中去除。

进一步的阅读

有关滤波和重采样的更多信息,请参阅信号处理工具箱。

参考译文:肯德尔、莫里斯·G.、艾伦·斯图尔特和J.基斯·奥德。高级统计理论,第3卷:设计和分析,和时间序列.第四版。伦敦:麦克米伦出版社,1983。

另请参阅

||||