信号平滑
这个例子展示了如何使用移动平均滤波器和重采样来隔离每小时温度读数中时间的周期性成分的影响,以及从开环电压测量中去除不必要的线路噪声。该示例还展示了如何通过使用中值滤波器在保持边缘的同时平滑时钟信号的电平。该示例还展示了如何使用Hampel过滤器去除较大的异常值。
动机
平滑是我们在数据中发现重要模式的方法,同时去掉不重要的东西(如噪声)。我们使用滤波来执行平滑。平滑的目标是产生值的缓慢变化,以便更容易看到数据中的趋势。
有时,当您检查输入数据时,您可能希望平滑数据,以便看到信号的趋势。在我们的例子中,我们有一组2011年1月在波士顿Logan机场每小时采集的温度读数(以摄氏度为单位)。
负载bostemp天=一句子* 24)/ 24;情节(天,tempC)轴紧ylabel (“临时(\ circC)”)包含(“时间从2011年1月1日(天)开始”)标题(“洛根机场干球温度(来源:NOAA)”)
注意,我们可以直观地看到一天中的时间对温度读数的影响。如果您只对一个月的日温度变化感兴趣,每小时的波动只会产生噪声,这可能会使每日的变化难以辨别。为了消除时间的影响,我们现在想通过使用移动平均过滤器平滑我们的数据。
移动平均滤波器
在其最简单的形式中,长度为N的移动平均滤波器取波形中每N个连续样本的平均值。
为了对每个数据点应用移动平均过滤器,我们构造过滤器的系数,使每个点的权重相等,占总平均值的1/24。这就得到了每24小时内的平均温度。
hoursPerDay = 24;coeff24hMA = ones(1, hoursPerDay)/hoursPerDay;avg24hTempC = filter(coeff24hMA, 1, tempC);情节(天,[tempC avg24hTempC])传说(每小时的临时的,“24小时平均(延迟)”,“位置”,“最佳”) ylabel (“临时(\ circC)”)包含(“时间从2011年1月1日(天)开始”)标题(“洛根机场干球温度(来源:NOAA)”)
过滤器延迟
注意,过滤后的输出延迟了大约12个小时。这是由于我们的移动平均滤波器有一个延迟。
任何长度为N的对称滤波器的时延都是(N-1)/2个样本。我们可以手工解释这个延迟。
fDelay =(长度(coeff24hMA) 1) / 2;情节(天,tempC,...avg24hTempC days-fDelay / 24日)轴紧传奇(每小时的临时的,24小时平均的,“位置”,“最佳”) ylabel (“临时(\ circC)”)包含(“时间从2011年1月1日(天)开始”)标题(“洛根机场干球温度(来源:NOAA)”)
提取的平均差异
或者,我们也可以使用移动平均滤波器来更好地估计一天中的时间是如何影响整体温度的。要做到这一点,首先,从每小时的温度测量中减去平滑数据。然后,将不同的数据划分为天,并取当月所有31天的平均值。
figure deltaTempC = tempC - avg24hTempC;deltaTempC =重塑(deltaTempC, 24, 31).';情节(桥,意味着(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值越高平滑效果越差。
α= 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);情节(天,天)传说(每小时的临时的,“马Cubic-Weighted”,“马Quartic-Weighted”,...“马Quintic-Weighted”,“位置”,“东南”) ylabel (“临时(\ circC)”)包含(“时间从2011年1月1日(天)开始”)标题(“洛根机场干球温度(来源:NOAA)”轴([3 5 -5 2])
重采样
有时,为了恰当地应用移动平均线,重采样信号是有益的。
在我们的下一个例子中,我们在存在60hz交流电源线噪声干扰的情况下,对模拟仪器输入端的开环电压进行采样。我们以1khz的采样率对电压进行采样。
负载openloop60hertzfs = 1000;t = (0:numel(openLoopVoltage)-1) / fs;情节(t, openLoopVoltage) ylabel (“电压(V)”)包含(“时间(s)”)标题(开环电压测量的)
让我们尝试通过使用移动平均滤波器来去除线噪声的影响。
如果你构造一个统一加权移动平均滤波器,它将去除任何与滤波器持续时间有关的周期性分量。
当以1000赫兹采样时,在60赫兹的完整周期中大约有1000 / 60 = 16.667个样本。让我们尝试“四舍五入”并使用17点过滤器。这将使我们在1000hz / 17 = 58.82 Hz的基频下得到最大的滤波。
情节(t, sgolayfilt (openLoopVoltage 1 17) ylabel (“电压(V)”)包含(“时间(s)”)标题(开环电压测量的)传说(“以58.82赫兹工作的移动平均滤波器”,...“位置”,“东南”)
注意,虽然电压是显著的平滑,它仍然包含一个小的60赫兹纹波。
如果我们重新采样信号,通过移动平均滤波器捕获60hz信号的完整完整周期,我们可以显著减少纹波。
如果我们以17 * 60hz = 1020 Hz重采样信号,我们可以使用我们的17点移动平均滤波器去除60hz的线噪声。
fsResamp = 1020;vResamp = resample(openLoopVoltage, fsResamp, fs);tResamp = (0:numel(vResamp)-1) / fsResamp;vAvgResamp = sgolayfilt (vResamp 1 17);情节(tResamp vAvgResamp) ylabel (“电压(V)”)包含(“时间(s)”)标题(开环电压测量的)传说(“移动平均滤波器工作在60赫兹”,...“位置”,“东南”)
中值滤波器
移动平均,加权移动平均,和Savitzky-Golay过滤器平滑所有他们过滤的数据。然而,这可能并不总是我们想要的。例如,如果我们的数据来自一个时钟信号,并且有我们不希望平滑的尖锐边缘,该怎么办?到目前为止讨论的过滤器工作得不是很好:
负载clockexyMovingAverage = conv (x)的(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)传说(原始信号的,“中值滤波”)
用汉普尔滤波器去除离群点
许多过滤器对异常值很敏感。与中值滤波器密切相关的一个滤波器是汉佩尔滤波器。该滤波器有助于从信号中去除异常值,而无需过度平滑数据。
为了看到这个,加载一个火车汽笛的音频记录,并添加一些人工噪音峰值:
负载火车y(1:400:结束)= 2.1;情节(y)
由于我们引入的每个尖峰的持续时间只有一个样本,我们可以使用一个只有三个元素的中值滤波器来去除尖峰。
持有在情节(medfilt1 (y, 3))从传奇(原始信号的,“过滤信号”)
滤波器去除了峰值,但它也去除了原始信号的大量数据点。Hampel滤波器的工作原理类似于中值滤波器,但是它只替换与局部中值相差几个标准差的值。
传说hampel (y, 13) (“位置”,“最佳”)
只有异常值被从原始信号中去除。
进一步的阅读
有关滤波和重采样的更多信息,请参见信号处理工具箱。
参考:肯德尔,莫里斯·G.,艾伦·斯图尔特和j·基思·奥德。高级统计理论,Vol. 3:设计和分析,时间序列.伦敦:麦克米伦出版社,1983年出版。