信号处理中简单实用的方法

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

消除趋势项函数

在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中清除趋势项

基线漂移的修正

1.概述

在实际采集到的信号中经常会发生基线漂移，这可能是由多种原因造成的。本案例将介绍在心电图检测中，若由于病人身体的轻微活动造成了心电图信号的基线漂移，则在后期处理中就可以用最小二乘法拟合消除消除基线的漂移。

2.实例

读入已知的心电图ecgdata2.mat数据，用最小二乘法拟合消除基线的漂移。程序清单如下：

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');

3.案例延伸

除了以最小二乘法拟合消除趋势项以外，还有其他方法可以消除趋势项，只是最小二乘法拟合的方法用得较多。其他的方法主要是一些平滑或滤波的方法，只取信号中的低频信号。例如可以用sgolay滤波器求取趋势项。sgolay滤波器是由SavitzkyA和GolayM在1964年提出的一种基于多项式拟合的最佳形式的低通滤波器[23,将在4.5.3小节中详细地说明。下面给出用sgolay滤波器对中的心电图数据消除趋势项的方法。读人已知的心电图ecgdata2.mat数据,用sgolay 滤波器消除基线的漂移。程序清单如下:

clear all; clc; close all;
load ecgdata2.mat                   % 读入心电图数据
N=length(y);                        % 数据长度
time=(0:N-1)/fs;                    % 计算出时间刻度
y1=sgolayfilt(y,3,2001);            % 用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('消除趋势项心电信号'); ylabe