信号平滑
这个例子展示了如何使用移动平均滤波器和重采样来隔离一天中时间的周期性分量对每小时温度读数的影响,以及从开环电压测量中去除不需要的线路噪声。该示例还展示了如何通过使用中值滤波器在保持边缘的同时平滑时钟信号的电平。该示例还展示了如何使用Hampel过滤器来删除较大的异常值。
动机
平滑是我们如何发现数据中的重要模式,同时忽略不重要的东西(即噪声)。我们使用过滤来执行平滑。平滑的目标是产生数值的缓慢变化,以便更容易看到数据中的趋势。
有时,当您检查输入数据时,您可能希望平滑数据,以便看到信号的趋势。在我们的例子中,我们有一组2011年1月在波士顿Logan机场每小时采集的温度读数,以摄氏度为单位。
负载bostemp天数= (1:31*24)/24;plot(days, tempC)轴紧ylabel (“临时(\ circC)”)包含(自2011年1月1日起已经过时间(天))标题(洛根机场干球温度(来源:NOAA))
请注意,我们可以直观地看到一天中的时间对温度读数的影响。如果你只对一个月的每日温度变化感兴趣,那么每小时的波动只会产生噪音,这可能会使每日的变化难以辨别。为了消除一天中时间的影响,我们现在想使用移动平均过滤器来平滑我们的数据。
移动平均滤波器
在其最简单的形式中,长度为N的移动平均滤波器取每N个连续波形样本的平均值。
为了对每个数据点应用移动平均过滤器,我们构造过滤器的系数,使每个点的权重相等,并对总平均值贡献1/24。这就得到了每24小时内的平均温度。
hoursPerDay = 24;coeff24hMA = ones(1, hoursPerDay)/hoursPerDay;avg24hTempC = filter(coeff24hMA, 1, tempC);plot(days,[tempC avg24hTempC])每小时的临时的,“24小时平均(延迟)”,“位置”,“最佳”) ylabel (“临时(\ circC)”)包含(自2011年1月1日起已经过时间(天))标题(洛根机场干球温度(来源:NOAA))
过滤器延迟
注意,过滤后的输出延迟了大约12个小时。这是因为我们的移动平均滤波器有一个延迟。
任何长度为N的对称滤波器都将具有(N-1)/2个样本的延迟。我们可以手动解释这个延迟。
fDelay = (length(coeff24hMA)-1)/2;情节(天,tempC,...avg24hTempC days-fDelay / 24日)轴紧传奇(每小时的临时的,“24小时平均”,“位置”,“最佳”) ylabel (“临时(\ circC)”)包含(自2011年1月1日起已经过时间(天))标题(洛根机场干球温度(来源:NOAA))
提取平均差异
或者,我们也可以使用移动平均滤波器来更好地估计一天中的时间如何影响整体温度。要做到这一点,首先,从每小时的温度测量中减去平滑数据。然后,将不同的数据分割成天,并对当月所有31天取平均值。
图deltaTempC = tempC - avg24hTempC;deltaTempC =重塑(deltaTempC, 24,31).';plot(1:24, mean(deltaTempC))轴紧标题(“24小时平均温度差”)包含(“一天中的几个小时(从午夜开始)”) ylabel (“温差(\circC)”)
提取峰值包络
有时,我们也希望对每天的温度信号的高低变化有一个平稳的估计。为此,我们可以使用信封
函数连接24小时内子集检测到的极端高点和低点。在本例中,我们确保每个极端高点和极端低点之间至少有16个小时。我们还可以通过取两个极端之间的平均值来了解高点和低点的趋势。
[envHigh, envLow] =信封(tempC,16,“高峰”);envMean = (envHigh+envLow)/2;情节(天,tempC,...天,envHigh,...天,envMean,...天,envLow)轴紧传奇(每小时的临时的,“高”,“的意思是”,“低”,“位置”,“最佳”) ylabel (“临时(\ circC)”)包含(自2011年1月1日起已经过时间(天))标题(洛根机场干球温度(来源:NOAA))
加权移动平均滤波器
其他类型的移动平均滤波器对每个样本的权重不相等。
的二项式展开是另一种常见的过滤器 .对于较大的n值,这种类型的滤波器近似于正态曲线。对于较小的n值,它对于滤除高频噪声很有用。要找到二项式滤波器的系数,请进行卷积 然后对输出进行迭代卷积 次数:规定次数在这个例子中,总共使用五次迭代。
H = [1/2 /2];binomialCoeff = conv(h,h);为n = 1:4 binomialCoeff = conv(binomialCoeff,h);结束figure fDelay = (length(binomialCoeff)-1)/2;binomialMA = filter(binomialCoeff, 1, tempC);情节(天,tempC,...binomialMA days-fDelay / 24日)轴紧传奇(每小时的临时的,“二项式加权平均”,“位置”,“最佳”) ylabel (“临时(\ circC)”)包含(自2011年1月1日起已经过时间(天))标题(洛根机场干球温度(来源:NOAA))
另一个类似于高斯展开滤波器的滤波器是指数移动平均滤波器。这种类型的加权移动平均滤波器易于构造,不需要大的窗口大小。
通过在0到1之间的alpha参数来调整指数加权移动平均滤波器。alpha值越高,平滑程度越低。
Alpha = 0.45;exponentialMA = filter(alpha, [1 alpha-1], tempC);情节(天,tempC,...binomialMA days-fDelay / 24日,...days-1/24 exponentialMA)轴紧传奇(每小时的临时的,...“二项式加权平均”,...“指数加权平均”,“位置”,“最佳”) ylabel (“临时(\ circC)”)包含(自2011年1月1日起已经过时间(天))标题(洛根机场干球温度(来源:NOAA))
放大一天的阅读量。
轴([3 4 -5 2])
Savitzky-Golay过滤器
您将注意到,通过平滑数据,极端值在某种程度上被剪切了。
为了更紧密地跟踪信号,您可以使用加权移动平均滤波器,它尝试在最小二乘意义上拟合指定数量样本上指定顺序的多项式。
为了方便起见,您可以使用该函数sgolayfilt
实现一个Savitzky-Golay平滑滤波器。使用sgolayfilt
,指定数据的奇数长度段和严格小于段长度的多项式顺序。的sgolayfilt
函数在内部计算平滑多项式系数,执行延迟对齐,并照顾数据记录开始和结束时的瞬态效应。
cubicMA = sgolayfilt(tempC, 3,7);quarticMA = sgolayfilt(tempC, 4,7);quinticMA = sgolayfilt(tempC, 5, 9);plot的用法和样例:每小时的临时的,“马Cubic-Weighted”,“马Quartic-Weighted”,...“马Quintic-Weighted”,“位置”,“东南”) ylabel (“临时(\ circC)”)包含(自2011年1月1日起已经过时间(天))标题(洛根机场干球温度(来源:NOAA))轴([3 5 -5 2])
重采样
有时为了正确地应用移动平均线,对信号重新采样是有益的。
在我们的下一个例子中,我们采样了一个模拟仪器输入端的开环电压,该电压存在来自60 Hz交流电源线噪声的干扰。我们以1khz的采样率对电压进行采样。
负载openloop60hertzFs = 1000;t =(0:数字(openLoopVoltage)-1) / fs;情节(t, openLoopVoltage) ylabel (“电压(V)”)包含(“时间(s)”)标题(“开环电压测量”)
让我们尝试通过使用移动平均滤波器来去除线噪声的影响。
如果你构造一个均匀加权的移动平均滤波器,它将删除任何与滤波器持续时间有关的周期性成分。
当以1000 Hz采样时,在60 Hz的完整周期中大约有1000 / 60 = 16.667个样本。让我们尝试“四舍五入”并使用17点过滤器。这将使我们在1000hz / 17 = 58.82 Hz的基频下获得最大滤波。
情节(t, sgolayfilt (openLoopVoltage 1 17) ylabel (“电压(V)”)包含(“时间(s)”)标题(“开环电压测量”)传说(“工作在58.82 Hz的移动平均滤波器”,...“位置”,“东南”)
注意,虽然电压明显平滑,但它仍然包含一个小的60 Hz纹波。
如果我们重新采样信号,通过移动平均滤波器捕获60 Hz信号的完整周期,我们可以显著减少纹波。
如果我们以17 * 60 Hz = 1020 Hz重新采样信号,我们可以使用17点移动平均滤波器去除60 Hz线噪声。
fsResamp = 1020;vResamp = resample(openLoopVoltage, fsResamp, fs);tResamp = (0: number (vResamp)-1) / fsResamp;vAvgResamp = sgolayfilt(vResamp,1,17);情节(tResamp vAvgResamp) ylabel (“电压(V)”)包含(“时间(s)”)标题(“开环电压测量”)传说(“移动平均滤波器工作在60赫兹”,...“位置”,“东南”)
中值滤波器
移动平均、加权移动平均和Savitzky-Golay滤波器平滑了他们过滤的所有数据。然而,这可能并不总是我们想要的。例如,如果我们的数据来自一个时钟信号,并且具有我们不希望平滑的尖锐边缘,该怎么办?到目前为止讨论的过滤器工作得不太好:
负载clockexyMovingAverage = conv(x,ones(5,1)/5,“相同”);ySavitzkyGolay = sgolayfilt(x,3,5);情节(t x,...t yMovingAverage...t, ySavitzkyGolay)传说(原始信号的,“移动平均”,“Savitzky-Golay”)
移动平均和Savitzky-Golay滤波器分别在时钟信号的边缘处进行了欠校正和过校正。
一个简单的方法来保持边缘,但仍然平滑的水平是使用中值过滤器:
yMedFilt = medfilt1(x,5,“截断”);情节(t x,...t, yMedFilt)传说(原始信号的,“中值滤波”)
通过Hampel滤波器去除异常值
许多过滤器对异常值很敏感。一个与中值滤波器密切相关的滤波器是Hampel滤波器。该滤波器有助于从信号中去除异常值,而不会过度平滑数据。
为了看到这一点,加载火车汽笛的音频记录,并添加一些人工噪音峰值:
负载火车Y (1:400:end) = 2.1;情节(y)
由于我们引入的每个尖峰的持续时间只有一个样本,所以我们可以使用只有三个元素的中值过滤器来去除尖峰。
持有在情节(medfilt1 (y, 3))从传奇(原始信号的,“过滤信号”)
过滤器去掉了尖峰,但也去掉了原始信号的大量数据点。Hampel过滤器的工作原理类似于中值过滤器,但是它只替换与局部中值相差几个标准差的值。
传说hampel (y, 13) (“位置”,“最佳”)
只有异常值从原始信号中去除。
进一步的阅读
有关滤波和重采样的更多信息,请参阅信号处理工具箱。
参考资料:肯德尔,莫里斯G.,艾伦斯图尔特,和J.基思奥德。高级统计理论,第3卷:设计和分析,和时间序列.第四版。伦敦:麦克米伦出版社,1983年。