matlab 函数调试实验,C语言实验matlab中filter和filtfilt函数

本学期的创造工房，课题为使用睡眠观测仪器所采集的混合了心拍数、呼吸数以及体动等杂音的原始数据，利用滤波器将心拍和呼吸数据抽出。本次工房我的分工为设计算法，由于时间仓促此次工房仅使用高低通滤波器进行简单滤波而未采用更为准确的小波算法，其实，利用matlab中butter、filter、filtfilt函数便可简单实现，但考虑到服务器部署matlab的成本，遂决定将此部分交由C程序来实现。

之前还考虑将butter函数也实现了，后来发现数学太差无从下手，不过因为特定滤波器中差分方程系数固定故此部分仍可由matlab实现。

滤波器设计如下

%0.1-0.4 Hz bandpass filter for respieration waveform

extraction

[RBh, RAh] = butter(2, 0.1/50, 'high');

[RBl, RAl] = butter(2, 0.4/50, 'low');

%1-3 Hz bandpass filter for pulse waveform extraction

[HBh, HAh] = butter(2, 1/50, 'high');

[HBl, HAl] = butter(2, 3/50, 'low');

%respiration waveform extraction

resp = filtfilt(RBh, RAh, data);

resp = filtfilt(RBl, RAl, resp);

%pulse waveform extraction

pulse = filtfilt(HBh, HAh, data);

pulse = filtfilt(HBl, HAl, pulse);

参考资料

filtfilt

http://www.mathworks.co.jp/help/toolbox/signal/ref/filtfilt.html

filterアルゴリズム

http://www.mathworks.co.jp/help/ja_JP/techdoc/ref/filter.html

C言语实现如下

#include "stdafx.h"

#include

#include

#include

#define EPS 0.000001

//filter函数

void filter(const double* x, double* y, int xlen, double* a,

double* b, int nfilt, double* zi)//nfilt为系数数组长度{

double tmp;

int i,j;

//normalization

if(

(*a-1.0>EPS) || (*a-1.0

){

tmp=*a;

for(i=0;i

b[i]/=tmp;

a[i]/=tmp;