利用matlab对波形进行去均值、去线性趋势和波形尖灭以及带通滤波

波形预处理

• • 介绍
  • 实例
  • • 去均值
    • 去线性趋势
    • 波形尖灭
    • 滤波
  • 函数源码
  • 参考资料

介绍

在处理波形数据时，常常需要对数据进行预处理，例如去均值，滤波等。本文利用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

参考资料

1. SAC Docs
2. ObsPy Tutorial