信号处理中简单实用的方法——消除趋势项

最小二乘法拟合消除趋势项

趋势项又分为线性趋势项和多项式趋势项，在许多实际信号获取后都有一些基线的漂移，这可能是采集系统引起的，也可能是信号本身引起的，但在信号处理之前要消除这种漂移，称为消除趋势项。
在实际处理中，信号中的趋势项往往比较复杂。本次介绍最小二乘法拟合消除多项式的趋势项。

消除趋势项函数

在MATLAB的工具箱中已有消除线性趋势项的detrend函数；本次介绍以最小二乘法拟合消除趋势项的polydetrend 函数。
函数：detrend
功能：消除线性趋势项
调用格式：y=detrend（x)
说明：输入参数x是带有线性趋势项的信号序列，输出参数y是消除趋势项的序列。
函数：polydetrend
功能：最小二乘法拟合消除多项式的趋势项
调用格式：[y,xtrend]=polydetrend(x,fs,m)
说明：输入参数x是带有趋势项的信号，fs是采样频率，m是调用本函数时设置的多项
式阶次。输出参数y是消除趋势项后的信号序列，xtrend是叠加在信号上的趋势项序列。
函数polydetrend的程序清单如下：

function [y,xtrend]=polydetrend(x, fs, m)
x=x(:);                 % 把语音信号x转换为列数据
N=length(x);            % 求出x的长度
t= (0: N-1)'/fs;        % 按x的长度和采样频率设置时间序列
a=polyfit(t, x, m);     % 用最小二乘法拟合语音信号x的多项式系数a
xtrend=polyval(a, t);   % 用系数a和时间序列t构成趋势项
y=x-xtrend;             % 从语音信号x中清除趋势项

基线漂移的修正案例：读人已知的心电图数据，用最小二乘法拟合消除基线的漂移。程序清单如下：

% 
clear all; clc; close all;

load ecgdata2.mat                   % 读入心电图数据
N=length(y);                        % 数据长度
time=(0:N-1)/fs;                    % 计算出时间刻度
[x,xtrend]=polydetrend(y, fs, 3);   % 用多项式拟合法求出趋势项及消除后的序列
% 作图
subplot 311; plot(time,y,'k')
title('输入心电信号'); ylabel('幅值');
axis([0 max(time) -2000 6000]); grid;
subplot 312; plot(time,xtrend,'k','linewidth',1.5);
title('趋势项信号'); ylabel('幅值');
axis([0 max(time) -2000 6000]); grid;
subplot 313; plot(time,x,'k'); 
title('消除趋势项心电信号'); ylabel('幅值');
xlabel('时间/s');
axis([0 max(time) -2000 6000]); grid;
set(gcf,'color','w');

运行结果如下图所示：

心电图数据链接如下：

https://download.csdn.net/download/qq_42233059/86405498

除了以最小二乘法拟合消除趋势项以外，还有其他方法可以消除趋势项，只是最小二乘法拟合的方法用得较多。其他的方法主要是一些平滑或滤波的方法，只取信号中的低频信号。例如可以用sgolay滤波器求取趋势项。sgolay滤波器是由Savitzky A和Golay M在1964年提出的一种基于多项式拟合的最佳形式的低通滤波器。下面给出用sgolay 滤波器对心电图数据消除趋势项的方法。

基线漂移的修正案例：读人已知的心电图数据，用sgolay 滤波器消除基线的漂移。程序清单下：

clear all; clc; close all;

load ecgdata2.mat                   % 读入心电图数据
N=length(y);                        % 数据长度
time=(0:N-1)/fs;                    % 计算出时间刻度
y1=sgolayfilt(y,3,1001);            % 用sgolay滤波器求出趋势项
x=y-y1;                             % 计算消除趋势项后的序列
% 作图
subplot 311; plot(time,y,'k')
title('输入心电信号'); ylabel('幅值');
axis([0 max(time) -2000 6000]); grid;
subplot 312; plot(time,y1,'k','linewidth',1.5);
title('趋势项信号'); ylabel('幅值');
axis([0 max(time) -2000 6000]); grid;
subplot 313; plot(time,x,'k'); 
title('消除趋势项心电信号'); ylabel('幅值');
xlabel('时间/s');
axis([0 max(time) -2000 6000]); grid;
set(gcf,'color','w');

运行结果如下：

线性函数拟合消除线性趋势项

基线漂移的修正案例：读人已知的一组实验数据，用线性函数拟合消除基线的漂移。程序清单如下：

clear all; clc; close all;

load qldata.mat                  % 读入数据
N=length(y);                     % 数据长度
time=(0:N-1)/fs;                 % 时间刻度
% 第一部分
Y=fft(y);                        % FFT
n2=1:N/2+1;                      % 取正频率索引序列
freq=(n2-1)*fs/N;                % 频率刻度
% 作图
subplot 211; plot(time,y,'k'); ylim([0 15]); grid;
title('有趋势项的数据')
xlabel('时间/s'); ylabel('幅值');
subplot 212; plot(freq,abs(Y(n2)),'k')
title('有趋势项的数据频谱')
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');
pause
% 第二部分
x=detrend(y);                    % 消除趋势项
X=fft(x);                        % FFT
% 作图
figure
subplot 211; plot(time,x,'k'); ylim([-5 5]); grid;
title('消除趋势项后的数据')
xlabel('时间/s'); ylabel('幅值');
subplot 212; plot(freq,abs(X(n2)),'k');
title('消除趋势项后的数据频谱')
xlabel('频率/Hz'); ylabel('幅值');
set(gcf,'color','w');

运行结果如下：

实验数据qldata.mat下载链接如下：

https://download.csdn.net/download/qq_42233059/86405503

参考文献：MATLAB数字信号处理85个实用案例精讲——入门到进阶；宋知用（编著）