# Matlab的filtfilt函数解析与C++实现

2 篇文章 0 订阅

## 1.filtfilt函数的解析

### 1.3 边界效应的优化

y ( i ) = ∑ j = 0 N b j x ( j ) − ∑ j = 1 N a j y ( j ) y(i) = \sum_{j=0}^N b_{j}x(j)-\sum_{j=1}^N a_{j}y(j)

y ( i ) = z i y ( 0 ) + ∑ j = 0 i b j x ( j ) − ∑ j = 1 i a j y ( j ) y(i) = z_iy(0)+\sum_{j=0}^i b_{j}x(j)-\sum_{j=1}^i a_{j}y(j)

y ( i ) = ∑ j = 0 N b j x ( j ) − ∑ j = 1 N a j y ( j ) y(i) = \sum_{j=0}^N b_{j}x(j)-\sum_{j=1}^N a_{j}y(j)

### 1.4 滤波器系数获取

nfilt = max(nb,na);
rows = [1:nfilt-1, 2:nfilt-1, 1:nfilt-2];
cols = [ones(1,nfilt-1), 2:nfilt-1, 2:nfilt-1];
vals = [1+a(2), a(3:nfilt)', ones(1,nfilt-2), -ones(1,nfilt-2)];
rhs  = b(2:nfilt) - b(1)*a(2:nfilt);
zi   = sparse(rows,cols,vals) \ rhs;
% The non-sparse solution to zi may be computed using:
%      zi = ( eye(nfilt-1) - [-a(2:nfilt), [eye(nfilt-2); ...
%                                           zeros(1,nfilt-2)]] ) \ ...
%          ( b(2:nfilt) - b(1)*a(2:nfilt) );


## 2.C++实现及对比

### 2.1 C++实现

.h文件

#pragma once
struct stFilterCoeff
{
int len;

double *a;
double *b;
double *zi;
};

template<typename T>
T getMax(T& a, T& b);

template<typename T>
T getMin(T& a, T& b);

template<typename T>
int getLength(T& a);

const char* writeToTxt(const char* name, double *datas, int lenN);

class FiltFilt //零相位滤波器
{
public:
FiltFilt(const int alen, const int blen, double *a, double *b, double *zi);
void filtfiltData(const double *sig, double *sout, const int data_len);

private:
stFilterCoeff myfilter;
void filtData(const double *sig, double *sout, const int data_len);
int signalExtend(const double *sig, double *sout, const int data_len, const int nfact);
void dataRollover(const double *sig, double *sout, const int data_len);

};


.c文件

#include "mFiltFilt.h"
FiltFilt::FiltFilt(const int alen, const int blen, double *a, double *b, double *zi)
{
myfilter.len = getMax(alen, blen); //使a和b同维度

myfilter.a = new double[myfilter.len];
myfilter.b = new double[myfilter.len];
myfilter.zi = new double[myfilter.len - 1];

for (int i = 0; i < myfilter.len; i++)
{
myfilter.a[i] = a[i];
myfilter.b[i] = b[i];
}
for (int i = 0; i < myfilter.len-1; i++)
{
myfilter.zi[i] = zi[i];

}
}

void FiltFilt::filtfiltData(const double *sig, double *sout, const int data_len)
{
// 零相位滤波
int N = myfilter.len;
int nfact = 3 * (N - 1);

int len_ext = data_len + 2 * nfact; //信号延拓后的长度

double *s_ext = new double[len_ext]; // 延拓后的信号（待滤波）
double *s_fft = new double[len_ext]; // 滤波后的信号

// 延拓 sig -> s_ext
signalExtend(sig, s_ext, data_len, nfact);
// 滤波 s_ext -> s_fft
filtData(s_ext, s_fft, len_ext);
// 翻转 s_fft -> s_ext
dataRollover(s_fft, s_ext, len_ext);
// 滤波 s_ext -> s_fft
filtData(s_ext, s_fft, len_ext);
// 翻转 s_fft -> s_ext
dataRollover(s_fft, s_ext, len_ext);
// 输出 s_ext -> sout
for (int i = 0; i < data_len; i++)
{
sout[i] = s_ext[i + nfact];
}
}

void FiltFilt::filtData(const double *sig, double *sout, const int data_len)
{
int N = myfilter.len;
for (int i = 0; i < data_len; i++)
{
sout[i] = sig[i];
//printf("sig[%d]=%f\n", i, sig[i]);
}
for (int i = 0; i < N-1; i++)
{
double tmp = myfilter.zi[i] * sig[0];
for (int j = 0; j < i+1; j++)
{
tmp += myfilter.b[j] * sig[i - j];
}
for (int j = 1; j < i+1; j++)
{
tmp -= myfilter.a[j] * sout[i - j];
}
sout[i] = tmp;
}
for (int i = N - 1; i < data_len; i++)
{
double tmp = 0;
for (int j = 0; j < myfilter.len; j++)
{
tmp += myfilter.b[j] * sig[i - j];
}
for (int j = 1; j < myfilter.len; j++)
{
tmp -= myfilter.a[j] * sout[i - j];
}
sout[i] = tmp;
//printf("out[%d]=%f\n",i, tmp);
}
}

int FiltFilt::signalExtend(const double *sig, double *sout, const int data_len, const int nfact)
{
// 数据延拓
int head_extN = nfact;
int tail_extN = nfact;
int data_extN = data_len + head_extN + tail_extN;
double data0 = 2*sig[0];
double data1 = 2*sig[data_len-1];
for (int i = 0; i < head_extN; i++)
{
sout[i] = data0-sig[head_extN - i];
}
for (int i = head_extN; i < data_len + head_extN; i++)
{
sout[i] = sig[i - head_extN];
}
for (size_t i = data_len + head_extN; i < data_extN; i++)
{
sout[i] = data1-sig[data_len + head_extN - i + data_len - 2];
}

return data_extN;

}

void FiltFilt::dataRollover(const double *sig, double *sout, const int data_len)
{
// 数据翻转
for (int i = 0; i < data_len; i++)
{
//printf("sig[%d]=%f\n", i, sig[i]);
sout[data_len - i - 1] = sig[i];
}
}


## 参考文献

• 9
点赞
• 49
收藏
觉得还不错? 一键收藏
• 16
评论
11-19
07-21 1万+
06-23 3万+
09-21 1904
01-02
09-25 4124
10-30
08-23

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。