Python零相位滤波 零相位滤波器在线设计

零相位(双边)滤波器设计–C++/Matlab

Date	Version	Comments
2019/11/12	V0.1	Init
2021/03/22	V0.2	修改延拓后滤波起点

Matlab滤波器设计

借助Matlab进行单边滤波器的设计比较简单，通过fdatool命令打开滤波器设计分析工具，按照自己的要求设计滤波器，点击Design Filter可以看到Filter Specifications的幅值响应。

Fs是采样频率，Fc是截止频率，选择了低通滤波器。点击菜单栏的Analysis-> Filter Coefficients导出滤波器系数，可以看到有3个Section,因为这里选择的滤波器为5阶，并且是直接Ⅱ型的，通过若干个两阶滤波器级联得到的。由于$5\div 2=2......1$，所以有三个Section。分别是传递函数为
$H(Z)=\frac{Gain\sum_{k=0}^{M}b[k]z^{-k}}{1+\sum_{k=1}^{N}a[k]z^{-k}}$
的三个系统串联得到的。其中，$M=N=3$。另外，通过Edit->Convert to Single Section可以将系数转化为一个多阶的形式。系数的长度将会是Order+1，也就是5阶时候，得到的a,b分别是6个数。

同时，获得滤波器系数的方式还有[b,a]=butter(5,0.2,'low');等方式，可以获得巴特沃兹滤波器，切比雪夫I/II型，椭圆滤波器，Bessel滤波器等的系数。可以输入help butter查看。以butter滤波器为例，第一个参数为阶数，第二个是截止频率（圆频率），采样频率的一半为1，假如采样率100Hz,那么一半时候(50Hz)为1，10Hz就是10/50=0.2，最后一个参数为滤波方式。与fdatool的方式得到的滤波器系数是一样的。

得到系数之后，可以通过Y=filter(b,a,x)的方式得到单边滤波结果。但是由于filter函数有边界效应，在信号有直流分量时候，在信号的起始处会有边界效应。
如果安装了WFDB的话，可以使用心电信号查看一下，也可以输入自己的信号查看滤波效果：

[signal,Fs,tm]=rdsamp('mitdb/100',[1],1000);
signal=signal+100;
[b,a]=butter(5,40/(Fs/2),'low');% 设计一个40Hz的滤波器
filtered=filter(b,a,signal);
h=figure();
subplot(211);
plot(signal,'b');hold on;plot(filtered,'LineWidth',2);legend('源信号','40Hz滤波');
subplot(212);
plot(signal,'b');hold on;plot(filtered,'LineWidth',2);legend('源信号','40Hz滤波');

在第二行中使用signal=signal+100;在信号上加一个直流偏置。效果如下图：

上图看到起始点基本上从零开始的，把上图局部放大得到下图，看到起始点与源信号差距很大，并且存在一个相位变化，整体的信号波峰向后移动了。

为了抑制这个边界问题还有相位变化，可以使用matlab的零相位滤波器函数filtfilt，例如bfiltered=filtfilt(b,a,signal);，滤波器系数还是之前的方式获得。

[signal,Fs,tm]=rdsamp('mitdb/100',[1],1000);
signal=signal+100;
[b,a]=butter(5,40/(Fs/2),'low');% 设计一个40Hz的滤波器
filtered=filter(b,a,signal);
plot(signal,'b');hold on;plot(filtered);
bfiltered=filtfilt(b,a,signal);
plot(bfiltered,'k','LineWidth',2);
legend('源信号','单边滤波器','双边滤波器');

得到的图为：

这样没有边界效果，也没有了相位问题。为了探究filtfilt函数是如何实现的，open filtfilt可以看到其中核心处理部分的代码是:

for ii=1:L
    % Single channel, data explicitly concatenated into one vector
    y = [2*y(1)-y(nfact+1:-1:2); y; 2*y(end)-y(end-1:-1:end-nfact)]; %#ok<AGROW>
    
    % filter, reverse data, filter again, and reverse data again
    y = filter(b(:,ii),a(:,ii),y,zi(:,ii)*y(1));
    y = y(end:-1:1);
    y = filter(b(:,ii),a(:,ii),y,zi(:,ii)*y(1));
    
    % retain reversed central section of y
    y = y(end-nfact:-1:nfact+1);
end

首先对信号向两端进行了延拓，为了遏制这个边界问题，然后用filter滤波处理一次，信号前后翻转过来，将翻转的信号送入filter进行处理，最后将信号翻转回来。

滤波器设计通过数学公式实现

在matlab中实现的核心代码是filter函数，但是在C语言中，需要使用滤波器的差分方程的方式实现，用数学方式实现的程序可在任何电脑上运行，甚至嵌入式中运行。差分方程可以写为：
$y(n)=\sum_{i=0}^{M}b_ix(n-i)+\sum_{i=1}^{N}a_iy(n-i)$
其中，$a$和$b$分别就是滤波器系数，$x(n)$为输入的第$n$个点，$y(n)$为输出的第$n$个点。由于模型的$a_i$都是负反馈,所以实际编程时候，都要做减法，而不是加法。以$M=N=2$,$n=0$为例计算：
$y(0)=b_0x(0)+b_1x(-1)+b_2x(-2)-a_1y(-1)-a_2y(-2)$
这个公式中，出现了-1和-2的下标，查看网上部分资源时候，直接将为负数的下标去掉，在信号本身在0附近时候，没有问题；但是信号有一个较大的直流偏量时候，就出现了起始位置的误差大。去边界效应一个简单有效的方式就是进行边界延拓，filtfilt的延拓方式是：

y = [2*y(1)-y(nfact+1:-1:2); y; 2*y(end)-y(end-1:-1:end-nfact)];

其实，没太明白为什么要这样延拓，在我写代码过程中，直接对称延拓效果也不错，由于在计算过程中只用到了当前信号点前的数据，没有后面的数据，所以这里不用向后延拓（因为这里是用的单边滤波）：

y=[y(nfact+1:-1:2),y];

延拓的长度nfact=3*(length(b)-1),发现，使用到的长度实际上没有这么长，nfact=2*(N-1);足够了，但是按理说，延拓得越长，那么对边界效应对信号起始段的影响越小。
用matlab试验一下效果,写一个自己的滤波器函数：

function s = mfilter(sig,b,a)
% 根据差分公式（在Latex中显示）：
% $$y(n)=\sum_{i=0}^{M}b_ix(n-i)+\sum_{i=1}^{N}a_iy(n-i)$$
% 进行滤波，滤波器系数b,a可以通过fdatool获得，或者
% [b,a]=butter(5,0.2,'low');等命令获得相应类型的滤波器系数
% 2019/11/11    V0.1
    if numel(a)~=numel(b)
        % 检查滤波器系数变量
        error(['滤波器系数长度不等！len(a)=',num2str(numel(a)),',len(b)=',num2str(numel(b))]);
    end
    % 信号延拓，防止信号有直流偏置时候，有很强的边界效应
    N=numel(a);
    nfact=2*(N-1);
    % 延拓的方式：上面一种参考了filtfilt函数的延拓方式，下面一种对称延拓
    %sig=[2*sig(1)-sig(nfact+1:-1:2),sig,2*sig(end)-sig(end-1:-1:end-nfact)];
    sig=[sig(nfact+1:-1:2),sig];
    s=sig;
    % 开始滤波
    % 由于信号延拓过了，在计算时候，只对中间的有效段进行的计算
    for i = N+1:length(s)-nfact
        tmp=0;
        for j=1:N
            tmp=tmp+b(j)*sig(i-j+1);
        end
        for j=2:N
            tmp = tmp - a(j)*s(i-j+1);
        end
        s(i)=tmp;
    end
    s=s(nfact+1:length(s)-nfact);
end

调用方式：

[signal,Fs,tm]=rdsamp('mitdb/100',[1],1000);
signal=signal+100;
[b,a]=butter(5,40/(Fs/2),'low');% 设计一个40Hz的滤波器
filtered=mfilter(signal',b,a);
plot(signal,'b');hold on;plot(filtered,'LineWidth',2);plot(filter(b,a,signal),'*');
legend('源信号','myFilter','系统filter');

那么可以看到，在信号起始段的误差问题已经没有了，并且，后面段的mfilter的图与filter函数滤波的图基本重叠了，没有任何区别。但是相位问题还是存在的，按照filtfilt的逻辑，滤波->翻转->滤波->翻转的方式来进行的话。写一个函数s=bmfilter(sig,b,a):

function s=bmfilter(sig,b,a)
% 零相位滤波,通过函数mfilter
% 滤波->翻转->滤波->翻转
% 2019/11/11     V0.1
    s=mfilter(sig,b,a);
    s=s(end:-1:1);
    s=mfilter(s,b,a);
    s=s(end:-1:1);
end

那么调用方式同样更新一下：

[signal,Fs,tm]=rdsamp('mitdb/100',[1],1000);
signal=signal+100;
[b,a]=butter(5,40/(Fs/2),'low');% 设计一个40Hz的滤波器
filtered=bmfilter(signal',b,a);
plot(signal,'b');hold on;plot(filtered,'LineWidth',2);plot(filter(b,a,signal));
legend('源信号','myFilter','系统filter');

C语言实现

原理没有问题，数学公式也非常清楚，那么写C语言代码就简单了，甚至可以通过MATLAB Coder自动转化为C语言代码（需要将变量检查部分注释）。

// mfilter.cpp : 定义控制台应用程序的入口点。
//

#include <iostream>
#include <fstream>
#include "butterFilter.h"

using namespace std;

int main()
{
	int len = 6;
	double a[6] = { 1, -2.75244113481752, 3.33871807933893, -2.13725305171939, 0.713552009831375, -0.0984430331073669 };
	double b[6] = { 0.00200415217268850, 0.0100207608634425, 0.0200415217268850, 0.0200415217268850, 0.0100207608634425, 0.00200415217268850 };
	ifstream infile("data.dat");
	ofstream outfile("filterdata.txt");
	ofstream outBi("filterdataBi.txt");
	double *data = new double[1000];
	double *fdata = new double[1000];
	for (int i = 0; i < 1000; i++)
	{
		infile >> data[i];
	}
	ButterFilter mybutterFilter(len,a,b);
	mybutterFilter.FiltData(data, fdata, 1000);
	printf("滤波完成!查看filterdata.txt文件!\n");
	for (int i = 0; i < 1000; i++)
	{
		outfile << fdata[i] << endl;
	}

	mybutterFilter.BiDirectionalFilter(data, fdata, 1000);
	printf("双向滤波完成!查看filterdataBi.txt文件!\n");
	for (int i = 0; i < 1000; i++)
	{
		outBi << fdata[i] << endl;
	}
	delete[]data;
	delete[]fdata;

	//getchar();
    return 0;
}

// ButterFilter.cpp
#include "butterFilter.h"

void ButterFilter::FiltData(const double * sig, double * sout,const int data_len)
{
	int N = myfilter.len;
	int nfact = 2 * (N - 1);
	//信号延拓
	double* s_ext = new double[data_len + nfact];	// 延拓后的输入信号
	double* tmps = new double[data_len + nfact];	// 作为处理过程中的输出信号
	int data_len_ext = signalExtend(sig, s_ext, data_len, nfact);
	for (int i = 0; i < data_len_ext; i++)
	{
		tmps[i] = s_ext[i];
	}
	for (int i = N; i < data_len_ext; i++)
	{
		double tmp = 0;
		for (int j = 0; j < N; j++)
		{
			tmp += myfilter.b[j] * s_ext[i - j];
		}
		for (int j = 1; j < N; j++)
		{
			tmp -= myfilter.a[j] * tmps[i - j];
		}
		tmps[i] = tmp;
		//printf("out[%d]=%f\n",i, tmp);
	}
	//输出
	for (int i = 0; i < data_len; i++)
	{
		sout[i] = tmps[i + nfact];
	}
	delete[]tmps;
	delete[]s_ext;
}
void ButterFilter::BiDirectionalFilter(const double * sig, double * sout, const int data_len)
{
	double* tmp1 = new double[data_len];
	double* tmp2 = new double[data_len];
	FiltData(sig, tmp1, data_len);
	dataRollingOver(tmp1, tmp2, data_len);
	FiltData(tmp2, tmp1, data_len);
	dataRollingOver(tmp1, sout, data_len);
	delete[]tmp1;
	delete[]tmp2;
}

int ButterFilter::signalExtend(const double *sig, double * sout, const int data_len, const int nfact)
{
	for (int i = 0; i < nfact; i++)
	{
		sout[i] = sig[nfact - i];
	}
	for (int i = nfact; i < data_len+ nfact; i++)
	{
		sout[i] = sig[i- nfact];
	}
	return data_len + nfact;
}


void ButterFilter::dataRollingOver(const double * sig,double * sout,const int data_len)
{
	for (int i = 0; i < data_len; i++)
	{
		sout[data_len - i-1] = sig[i];
	}
}

ButterFilter::ButterFilter()
{
}

ButterFilter::ButterFilter(const int len, double *a, double *b)
{
	myfilter.len = len;
	myfilter.a = new double(myfilter.len);
	myfilter.b = new double(myfilter.len);
	
	for (int i = 0; i < len; i++)
	{
		myfilter.a[i] = a[i];
		myfilter.b[i] = b[i];
	}
}

ButterFilter::~ButterFilter()
{
}

//butterFilter.h
#pragma once
struct FilterCoeff
{
	int len;
	double *a;
	double *b;
};

class ButterFilter
{
public:
	
	ButterFilter();
	ButterFilter(const int len, double * a, double * b);
	~ButterFilter();
	void FiltData(const double * sig, double * sout, const int data_len);
	void BiDirectionalFilter(const double * sig, double * sout, const int data_len);

private:
	FilterCoeff myfilter;
	int signalExtend(const double * sig, double * sout, const int data_len, const int nfact);
	void dataRollingOver(const double * sig, double * sout, const int data_len);

};

并且将matlab的数据文件保存下来，保存为data.dat。使用g++ mfilter.cpp butterFilter.cpp -g -o mfilter.exe命令编译为可执行文件后即可执行。

PS D:\me\Desktop\f> g++ mfilter.cpp butterFilter.cpp -g -o mfilter.exe
PS D:\me\Desktop\f> .\mfilter.exe
滤波完成!查看filterdata.txt文件!
双向滤波完成!查看filterdataBi.txt文件!
PS D:\me\Desktop\f>

PS

程序在Windows使用Visual Studio编译时候，发现有处有未经处理的异常: 堆已损坏的问题，网上说是数据越界操作，查找了好久都没有看到问题。如果有朋友能解决了，烦请告诉一声。

Python

Python中可以使用数学方式做，也可以用scipy包的函数。y = scipy.signal.filtfilt(b, a, data),其中b,a可以由b, a = butter(filter_order, [low, high], btype="band")获得。

signal_freq=250
	lowcut =0.2
	lowcut =25
    nyquist_freq = 0.5 * signal_freq
    low = lowcut / nyquist_freq
    high = lowcut / nyquist_freq
    b, a = butter(filter_order, [low, high], btype="band")
    y = scipy.signal.filtfilt(b, a, data)

参考：

1. 单片机中（C语言）IIR滤波器的实现，
2. IIR数字滤波器的实现（C语言），