本文详细介绍了基于C语言的巴特沃思IIR滤波器设计，包括滤波器类型、巴特沃思滤波器原型、C代码实现及其与MATLAB结果的比较，特别关注了零极点计算和双线性变换过程。
参考文档
主要参考资料来源数字信号处理 基于计算机的方法 4版（中文版）和Matlab函数butter(巴特沃斯IIR滤波器) 两部分。
IIR滤波器和分类
-
IIR滤波器
即无限冲激响应的滤波器
缺点 是 相位响应非线性（相位畸变）且在高阶时不稳定（过度或发散响应）
优点 是相较于FIR而言高效、实时、频率选择等。 -
分类
如今主要的IIR滤波器有以下几类：
（这篇文章主要介绍巴特沃斯IIR滤波器的C实现）
C代码
void butter(int n, double Wn[], int type, int analog, double *ab, double *bb)
{
double fs = 2;
double u[2] = {0.0, 0.0};
// step 1: get analog, pre-warped frequencies
if (!analog)
{
if (type == 1 || type == 2)
{
fs = 2;
u[0] = 2 * fs * tan(PI * Wn[0] / fs);
}
else
{
fs = 2;
u[0] = 2 * fs * tan(PI * Wn[0] / fs);
u[1] = 2 * fs * tan(PI * Wn[1] / fs);
}
}
else if (type == 3 || type == 4)
{
if (type == 1 || type == 2)
{
u[0] = Wn[0];
}
else
{
u[1] = Wn[1];
}
}
// step 2: convert to low-pass prototype estimate
double Bw = 0.0;
if (type == 1 || type == 2)
{
Wn = u;
}
else if (type == 3 || type == 4)
{
Bw = u[1] - u[0];
Wn[0] = sqrt(u[0] * u[1]); // center
Wn[1] = 0.0;
}
// step 3: Get N-th order Butterworth analog lowpass prototype
Complex *z = (Complex *)malloc(0 * sizeof(Complex));
Complex *p = (Complex *)malloc(n * sizeof(Complex));
double k = 0;
buttap(n, p, &k);
// Transform to state-space
int a_size = n;
double *a = (double *)malloc(sizeof(double) * n * n);
double *b = (double *)malloc(sizeof(double) * n);
double *c = (double *)malloc(sizeof(double) * n);
double d;
zp2ss(z, p, 0, n, k, a, b, c, &d);
switch (type)
{
default:
case 1: // Lowpass
lp2lp(n, a, b, Wn[0]);
break;
case 2: // highpass
lp2hp(n, a, b, c, &d, Wn[0]);
break;
case 3: // Bandpass
lp2bp(n, &a, &b, &c, &d, Wn[0], Bw);
a_size = 2 * n;
break;
case 4: // Bandstop
lp2bs(n, &a, &b, &c, &d, Wn[0], Bw);
a_size = 2 * n;
break;
}
if (!analog)
{
bilinear(a_size, a, b, c, &d, fs);//双线性变换
}
hhbg(a, a_size);
double *u_real = (double *)malloc(sizeof(double) * a_size);
double *v_imag = (double *)malloc(sizeof(double) * a_size);
int jt = 60;
hhqr(a, a_size, EPS, jt, u_real, v_imag);//EPS为2.24E-18;
Complex *p1 = (Complex *)malloc(sizeof(Complex) * a_size);
double *z1 = (double *)malloc(sizeof(double) * a_size);
double *polyz1 = (double *)malloc(sizeof(double) * (a_size+1));
Complex *ctemp = (Complex *)malloc(sizeof(Complex) * (a_size + 1));
double k0=0;
for (int i = 0; i < a_size; i++)
{
p1[i].real = u_real[i];
p1[i].imaginary = v_imag[i];
}
complex_poly(p1, a_size, ctemp);
for (int j = 0; j < a_size + 1; j++)
{
ab[j] = get_real(ctemp[j]);
}
buttzeros(p1,a_size,type,n,Wn[0],z1,&k0);
vector_poly(z1,a_size,polyz1);
for (size_t i = 0; i < a_size+1; i++)
{
bb[i]=k0*polyz1[i];
}
free(z);
free(p);
free(a);
free(b);
free(u_real);
free(v_imag);
free(ctemp);
free(p1);
free(z1);
free(polyz1);
}
函数中出现比较多的数学计算后面会文章会一一罗列，这里主要介绍巴特沃斯相关的函数
其中buttap为Butterworth filter prototype,创建N阶巴特沃斯滤波器的零极点模型，代码如下：
void buttap(int n, Complex *p, double *k)
{
int size = (n - 1) / 2;
int tempsize=floor(n/2);
Complex *ptemp=(Complex *)malloc(sizeof(Complex)*tempsize);
int index=0;
for (size_t i = 1; i <= n-1; i+=2)
{
double pp=PI*i/(2.0*n)+PI/2.0;
ptemp[index]=complex_multiply_double(create(0,1.0),pp);
ptemp[index]=complex_exponent(ptemp[index]);
index++;
}
for (size_t i = 0; i < n; i++)
{
p[i] = create(0, 1);
}
if (n % 2)
{
p[n - 1] = create(-1, 0);
}
for (size_t k = 0; k < tempsize; k++)
{
p[2*k]=ptemp[k];
p[2*k+1]=complex_conj(ptemp[k]);
}
*k =complex_array_prod(complex_array_multiply_double(p,-1.0,n),n).real;
free(ptemp);
}
其中buttzeros用于计算零点（zeros）和增益(gain)，原代码如下：
function [z,k] = buttzeros(btype,n,Wn,analog,p)
% This internal function computes the zeros and gain of the ZPK
% representation. Wn is scalar (sqrt(Wn(1)*Wn(2)) for bandpass/stop).
if analog
% for lowpass and bandpass, don't include zeros at +Inf or -Inf
switch btype
case 1 % lowpass: H(0)=1
z = zeros(0,1,'like',p);
k = Wn^n; % prod(-p) = Wn^n
case 2 % bandpass: H(1i*Wn) = 1
z = zeros(n,1,'like',p);
k = real(prod(1i*Wn-p)/(1i*Wn)^n);
case 3 % highpass: H(Inf) = 1
z = zeros(n,1,'like',p);
k = 1;
case 4 % bandstop: H(0) = 1
z = 1i*Wn*((-1).^(0:2*n-1)');
k = 1; % prod(p) = prod(z) = Wn^(2n)
otherwise
coder.internal.error('signal:iirchk:BadFilterType','high','stop','low','bandpass');
end
else
Wn = 2*atan2(Wn,4);
switch btype
case 1 % lowpass: H(1)=1
z = -ones(n,1,'like',p);
k = real(prod(1-p))/2^n;
case 2 % bandpass: H(z) = 1 for z=exp(1i*sqrt(Wn(1)*Wn(2)))
z = [ones(n,1,'like',p); -ones(n,1,'like',p)];
zWn = exp(1i*Wn);
k = real(prod(zWn-p)/prod(zWn-z));
case 3 % highpass: H(-1) = 1
z = ones(n,1,'like',p);
k = real(prod(1+p))/2^n;
case 4 % bandstop: H(1) = 1
z = exp(1i*Wn*( (-1).^(0:2*n-1)' ));
k = real(prod(1-p)/prod(1-z));
otherwise
coder.internal.error('signal:iirchk:BadFilterType','high','stop','low','bandpass');
end
end
% Note: codegen complains when z set to both real and complex values above
if ~any(imag(z))
z = real(z);
end
end
C实现代码如下：
void buttzeros(Complex *p,int np,int type,int n,double wn,double *z,double *k)
{
Complex *zt=(Complex*)malloc(sizeof(Complex)*np);
Complex *pt=(Complex*)malloc(sizeof(Complex)*np);
Complex zwn;
wn=2*atan2(wn,4);
switch (type)
{
case 1:
for (size_t i = 0; i < np; i++)
{
z[i]=-1.0;
pt[i] = complex_subtract(create(1.0,0.0),p[i]);
}
*k=get_real(complex_array_prod(pt,np))/pow(2,n);
break;
case 2:
for (size_t i = 0; i < np; i++)
{
z[i] = -1.0;
pt[i] = complex_add(create(1.0, 0.0), p[i]);
}
*k = get_real(complex_array_prod(pt, np)) / pow(2, n);
break;
case 3:
zwn = complex_exponent(complex_multiply_double(create(0, 1.0), wn));
for (size_t i = 0; i < np; i++)
{
if (i < n)
{
z[i] = 1.0;
}
else
{
z[i] = -1.0;
}
zt[i] = complex_subtract_double(zwn, z[i]);
pt[i] = complex_subtract(zwn, p[i]);
}
*k = get_real(complex_divide(complex_array_prod(pt, np), complex_array_prod(zt, np)));
break;
case 4:
for (size_t i = 0; i < np; i++)
{
if (i < n)
{
z[i] = 1.0;
}
else
{
z[i] = -1.0;
}
zt[i] = complex_subtract_double(create(1.0,0.0), z[i]);
pt[i] = complex_subtract(create(1.0,0.0), p[i]);
}
*k = get_real(complex_divide(complex_array_prod(pt, np), complex_array_prod(zt, np)));
break;
default:
break;
}
free(zt);
free(pt);
}
测试
于matlab结果对比，当n>500阶时，会在计算过程中出现个位数差距，实现过程基本无误。
下一篇: matlab cheby1 C实现
扫一扫
1646
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
