介绍
在处理波形数据时,常常需要对数据进行预处理,例如去均值,滤波等。本文利用matlab,通过实例来介绍常见的几种预处理方法:去均值、去线性趋势和波形尖灭以及带通滤波。
去均值:去除波形数据的平均值。
去线性趋势:将数据拟合成一条直线,然后从数据中减去该直线所表征的线性趋势。
波形尖灭:将波形数据的首尾两端由其原始值不断光滑地减小到0。
带通滤波:只保留特定频段的波形,同时屏蔽其他频段的波形。
实例
首先,我们给出一个原始波形:
dt = 0.01;
t = [0:dt:10-dt]';
data = 10*sin(2*pi*t)+8*cos(8*pi*t)+1*t+10;
该波形由一个振幅为10,周期为1秒的正弦函数、一个振幅为8,周期为0.25秒的余弦函数以及一个一次函数组合而成。
此时波形的平均值为:
>> mean(data)
ans =
10.4995
去均值
这里我们直接调用函数datapreproc
来实现,函数源码见下文。
通过指定rmean
为true来去除均值。
preproc_data = datapreproc(data,dt,'rmean',true,'rtrend',false,'tap_width',0);
figure(1)
subplot(211)
plot(t,data)
ylim([-40 40])
title('raw')
subplot(212)
plot(t,preproc_data)
ylim([-40 40])
title('rmean')
此时波形的平均值为:
>> mean(preproc_data)
ans =
-4.2011e-15
近似为0,但仍然可以看到线型趋势(一次函数的作用)。
去线性趋势
通过指定rtrend
为true来去除线性趋势。
preproc_data = datapreproc(data,dt,'rmean',true,'rtrend',true,'tap_width',0);
figure(1)
subplot(211)
plot(t,data)
ylim([-40 40])
title('raw')
subplot(212)
plot(t,preproc_data)
ylim([-40 40])
title('rmean, rtrend')
可以看出,此时波形从左往右缓缓增加的线性趋势
1
∗
t
1*t
1∗t也被去除了。
波形尖灭
通过指定tap_width
的值来指定波形尖灭的宽度,当tap_width
指定为0时,不进行波形尖灭。
preproc_data = datapreproc(data,dt,'rmean',true,'rtrend',true,'tap_width',0.1);
preproc_data2 = datapreproc(data,dt,'rmean',true,'rtrend',true,'tap_width',0.3);
figure(1)
subplot(311)
plot(t,data)
ylim([-40 40])
title('raw')
subplot(312)
plot(t,preproc_data)
ylim([-40 40])
title('rmean, rtrend, tap_width=0.1','Interpreter','none')
subplot(313)
plot(t,preproc_data2)
ylim([-40 40])
title('rmean, rtrend, tap_width=0.3','Interpreter','none')
滤波
滤波可能不能算作波形“预”处理的一部分,但也是非常常见的波形处理方法。
通过指定filtband
的值[freqmin freqmax]来对波形进行滤波。
filtband
的值指定为一个长度为2的一维数组,第一个数表示通带滤波的低转角频率(low corner frequency),第二个数表示通带滤波的高转角频率(high corner frequency)。
preproc_data = datapreproc(data,dt,'rmean',true,'rtrend',true,'filtband',[0.9,1.1],'tap_width',0.1);
figure(1)
subplot(211)
plot(t,data)
ylim([-40 40])
title('raw')
subplot(212)
plot(t,preproc_data)
ylim([-40 40])
title('rmean, rtrend, filtered: 0.9-1.1 Hz')
此时滤除了频率为4 Hz(周期为0.25秒)的余弦函数,保留了频率为1 Hz的正弦函数的波形信息。
函数源码
function preproc_data = datapreproc(raw_data,dt,varargin)
% Data preprocessing
% usage: dataout = datapreproc(datain,sampling interval,'filtband',[0.02 0.05],'tap_width',0)
% usage: dataout = datapreproc(datain,sampling interval,'filtband',[0.01 0.1],'rmean',false,'rtrend',false)
% perform rmean, rtrend and taper by default
% tap_width can be 0, then the taper operation will not be executed.
%
% Yuechu Wu
% 12131066@mail.sustech.edu.cn
% 2022-09-29
%
% added default parameters
% 2023-03-02, Yuechu Wu
% set up input parser
p = inputParser;
addParameter(p,'filtband',0);
addParameter(p,'tap_width',0.01);
addParameter(p,'rmean',true);
addParameter(p,'rtrend',true);
parse(p,varargin{:});
filtband = p.Results.filtband;
tap_width = p.Results.tap_width;
rmean = p.Results.rmean;
rtrend = p.Results.rtrend;
raw_data(isnan(raw_data)) = 0; % change NaN to 0
if rmean
rmean_data = raw_data - mean(raw_data);
else
rmean_data = raw_data;
end
if rtrend
rtrend_data = detrend(rmean_data);
else
rtrend_data = rmean_data;
end
if filtband
filt_data = bpfilt(rtrend_data,dt,filtband(1),filtband(2));
else
filt_data = rtrend_data;
end
if tap_width
taper_data = filt_data.*tukeywin(length(filt_data),tap_width);
preproc_data = taper_data;
else
preproc_data = filt_data;
end
return
function filt_data = bpfilt(data,dt,lowfreq,highfreq)
% Bandpass filtering of time series.
% usage: filt_data = bpfilt(data, sampling interval [s], low frequency [Hz], high frequency [Hz])
%
% Yuechu Wu
% 12131066@mail.sustech.edu.cn
% 2022-10-08
fn = 1/2/dt;
[b,a] = butter(2,[lowfreq/fn,highfreq/fn]);
filt_data = filtfilt(b,a,data);
return