|
|
相关文章推荐
/*********************************************************************************************************
* 函数名称： Gain
* 函数功能： 数字滤波器的频率响应
* 输入参数： b ：双精度实型一维数组，长度为（m+1）。存放滤波器分子多项式的系数 b(i)。
* a ：双精度实型一维数组，长度为（n+1）。存放滤波器分母多项式的系数 a(i)。
* m ：整型变量。滤波器分子多项式的阶数。
* n ：整形变量。滤波器分母多项式的阶数。
* x ：双精度实型一维数组，长度为 len。
* 当 sign = 0 时，存放滤波器频率响应的实部 Rc[H(w)]；
* 当 sign = 1 时，存放滤波器幅频响应 |H(w)|；
* 当 sign = 2 时，存放分贝表示的滤波器幅频响应 |H(w)|。
* y ：双精度实型一维数组，长度为 len。
* 当 sign = 0 时，存放滤波器频率响应的虚部 Im[H(w)]；
* 当 sign = 1 或 2 时，存放滤波器的相频响应 φ(w)。
* sign：整形变量。
* 当 sign = 0 时，计算滤波器频率响应的实部 Rc[H(w)] 和 虚部 Im[H(w)]；
* 当 sign = 1 时，计算滤波器的幅频响应 |H(w)| 和相频响应 φ(w)；
* 当 sign = 2 时，计算滤波器的幅频响应 |H(w)| （用 dB 表示）和相频响应 φ(w)。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年09月29日
* 注 意：
*********************************************************************************************************/
void Gain(double* b, double* a, int m, int n, double* x, double* y, int len, int sign)
int i, k;
double ar, ai, br, bi, zr, zi, im, re, den, numr, numi, freq, temp;
for (k = 0; k < len; k++)
freq = k * 0.5 / (len - 1);
zr = cos(-8.0 * atan(1.0) * freq);
zi = sin(-8.0 * atan(1.0) * freq);
br = 0.0;
bi = 0.0;
for (i = m; i > 0; i--)
re = br;
im = bi;
br = (re + b[i]) * zr - im * zi;
bi = (re + b[i]) * zi + im * zr;
ar = 0.0;
ai = 0.0;
for (i = n; i > 0; i--)
re = ar;
im = ai;
ar = (re + a[i]) * zr - im * zi;
ai = (re + a[i]) * zi + im * zr;
br = br + b[0];
ar = ar + 1.0;
numr = ar * br + ai * bi;
numi = ar * bi - ai * br;
den = ar * ar + ai * ai;
x[k] = numr / den;
y[k] = numi / den;
switch (sign)
case 1:
temp = sqrt(x[k] * x[k] + y[k] * y[k]);
y[k] = atan2(y[k], x[k]);
x[k] = temp;
break;
case 2:
temp = x[k] * x[k] + y[k] * y[k];
y[k] = atan2(y[k], x[k]);
x[k] = 10.0 * log10(temp);
数字系统的传递函数为
\(H(z) = \frac{-0.1z^{-1}}{1 + 0.9z^{-2}}\)
求该系统的幅频响应和相频响应，并画出相应的波形。
int main(void)
int i;
static double a[] = { 1.0, 0.0, 0.9 }, b[] = { 0.0, -0.1 };
static double f, x[300], y[300];
FILE* fp;
Gain(b, a, 1, 2, x, y, 300, 1);
fp = fopen("gainam.csv", "w");
for (i = 0; i < 300; i++)
f = i * 0.5 / 299;
fprintf(fp, "%lf,%lf\n", f, x[i]);
fclose(fp);
fp = fopen("gainph.csv", "w");
for (i = 0; i < 300; i++)
f = i * 0.5 / 299;
fprintf(fp, "%lf,%lf\n", f, y[i]);
fclose(fp);
return 0;
输出的幅频响应曲线如下图所示.
/*********************************************************************************************************
* 函数名称： Gainc
* 函数功能： 级联型数字滤波器的频率响应
* 输入参数： b ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分子多项式的系数，
* b[j][i] 表示第 j 个 n 阶节的分子多项式的第 i 个系数。
* a ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分母多项式的系数，
* a[j][i] 表示第 j 个 n 阶节的分母多项式的第 i 个系数。
* n ：整型变量。级联型滤波器每节的阶数。
* ns ：整形变量。级联型滤波器的 n 阶节数 L。
* x ：双精度实型一维数组，长度为 len。
* 当 sign = 0 时，存放滤波器频率响应的实部 Rc[H(w)]；
* 当 sign = 1 时，存放滤波器幅频响应 |H(w)|；
* 当 sign = 2 时，存放分贝表示的滤波器幅频响应 |H(w)|。
* y ：双精度实型一维数组，长度为 len。
* 当 sign = 0 时，存放滤波器频率响应的虚部 Im[H(w)]；
* 当 sign = 1 或 2 时，存放滤波器的相频响应 φ(w)。
* sign：整形变量。
* 当 sign = 0 时，计算滤波器频率响应的实部 Rc[H(w)] 和 虚部 Im[H(w)]；
* 当 sign = 1 时，计算滤波器的幅频响应 |H(w)| 和相频响应 φ(w)；
* 当 sign = 2 时，计算滤波器的幅频响应 |H(w)| （用 dB 表示）和相频响应 φ(w)。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年09月29日
* 注 意：
*********************************************************************************************************/
void Gainc(double* b, double* a, int n, int ns, double* x, double* y, int len, int sign)
int i, j, k, n1;
double ar, ai, br, bi, zr, zi, im, re, den, numr, numi, freq, temp;
double hr, hi, tr, ti;
n1 = n + 1;
for (k = 0; k < len; k++)
freq = k * 0.5 / (len - 1);
zr = cos(-8.0 * atan(1.0) * freq);
zi = sin(-8.0 * atan(1.0) * freq);
x[k] = 1.0;
y[k] = 0.0;
for (j = 0; j < ns; j++)
br = 0.0;
bi = 0.0;
for (i = n; i > 0; i--)
re = br;
im = bi;
br = (re + b[j * n1 + i]) * zr - im * zi;
bi = (re + b[j * n1 + i]) * zi + im * zr;
ar = 0.0;
ai = 0.0;
for (i = n; i > 0; i--)
re = ar;
im = ai;
ar = (re + a[j * n1 + i]) * zr - im * zi;
ai = (re + a[j * n1 + i]) * zi + im * zr;
br = br + b[j * n1 + 0];
ar = ar + 1.0;
numr = ar * br + ai * bi;
numi = ar * bi - ai * br;
den = ar * ar + ai * ai;
hr = numr / den;
hi = numi / den;
tr = x[k] * hr - y[k] * hi;
ti = x[k] * hi + y[k] * hr;
x[k] = tr;
y[k] = ti;
switch(sign)
case 1:
temp = sqrt(x[k] * x[k] + y[k] * y[k]);
if (temp != 0.0)
y[k] = atan2(y[k], x[k]);
y[k] = 0.0;
x[k] = temp;
break;
case 2:
temp = x[k] * x[k] + y[k] * y[k];
if (temp != 0.0)
y[k] = atan2(y[k], x[k]);
temp = 1.0e-40;
y[k] = 0.0;
x[k] = 10.0 * log10(temp);
数字滤波器由两个子系统 \(H_{1}(z)\) 和 \(H_{2}(z)\) 级联而成，其中 \(H_{1}(z)\) 为
\(H_{1}(z) = \frac{0.2}{1 – 0.8z^{-1}}\)
\(H_{2}(z)\) 为
\(H_{2}(z) = \frac{0.05 – 0.1z^{-1}}{1 – 0.95z^{-1} + 0.9z^{-2}}\)
选取参数 n = 2，ns = 2，len = 300，求该滤波器的幅频响应。
测试代码如下：
int main(void)
int i, n, ns;
double f;
static double a[2][3] = { {1.0, -0.8, 0.0}, {1.0, -0.95, 0.9} };
static double b[2][3] = { {0.2, 0.0, 0.0}, {0.05, -0.1, 0.0} };
static double x[300] = { 0 };
static double y[300] = { 0 };
FILE* fp;
n = 2;
ns = 2;
Gainc(b, a, n, ns, x, y, 300, 1);
fp = fopen("gaincam.csv", "w");
for (i = 0; i < 300; i++)
f = i * 0.5 / 299;
fprintf(fp, "%lf,%lf\n", f, x[i]);
fclose(fp);
return 0;
输出的幅频响应曲线如下所示。
/*********************************************************************************************************
* 函数名称： Resp
* 函数功能： 数字滤波器的时域响应
* 输入参数： x ：双精度实型一维数组，长度为 lx。存放滤波器的输入序列。
* y ：双精度实型一维数组，长度为 ly。存放滤波器的输出序列。
* lx：整型变量。输入序列的长度。
* ly：整形变量。输出序列的长度。
* b ：双精度实型一维数组，长度为（m+1）。存放滤波器分子多项式的系数 b(i)。
* a ：双精度实型一维数组，长度为（n+1）。存放滤波器分母多项式的系数 a(i)。
* m ：整形变量。滤波器分子多项式的阶数。
* n ：整形变量。滤波器分母多项式的阶数。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年09月29日
* 注 意：
*********************************************************************************************************/
void Resp(double* x, double* y, int lx, int ly, double* b, double* a, int m, int n)
int k, i, i1;
double sum;
for (k = 0; k < ly; k++)
sum = 0.0;
for (i = 0; i <= m; i++)
if ((k - i) >= 0)
i1 = ((k - i) < lx) ? (k - i) : (lx - 1);
sum = sum + b[i] * x[i1];
for (i = 1; i <= n; i++)
if ((k - i) >= 0)
sum = sum - a[i] * y[k - i];
y[k] = sum;
数字滤波器由两个子系统 \(H_{1}(z)\) 和 \(H_{2}(z)\) 级联而成，其中 \(H_{1}(z)\) 为
\(H_{1}(z) = \frac{0.05 – 0.1z^{-1}}{1 – 0.95z^{-1} + 0.9z^{-2}}\)
\(H_{1}(z)\) 为
\(H_{2}(z) = \frac{0.2}{1 – 0.8z^{-1}}\)
选取参数 lx = 2，ly = 100，试求高滤波器的单位阶跃响应。
测试代码如下：
int main(void)
int i, m1, m2, n1, n2, lx, ly;
static double y[100];
static double z[100];
static double x[2] = { 1.0, 1.0 };
static double b[2] = { 0.05, -0.1 };
static double a[3] = { 1.0, -0.95, 0.9 };
static double d[1] = { 0.2 };
static double c[2] = { 1.0, -0.8 };
FILE* fp;
lx = 2;
ly = 100;
m1 = 1;
n1 = 2;
m2 = 0;
n2 = 1;
fp = fopen("resp.csv", "w");
Resp(x, y, lx, ly, b, a, m1, n1);
Resp(y, z, ly, ly, d, c, m2, n2);
for (i = 0; i < ly; i++)
fprintf(fp, "%4d,%lf\n", i, z[i]);
fclose(fp);
for (i = 0; i < 32; i++)
printf(" %10.7lf", z[i]);
if ((i % 4) == 3)
printf("\n");
return 0;
终端输出的单位阶跃响应如下所示。
0.0100000 0.0075000 -0.0134750 -0.0388313
-0.0501862 -0.0430680 -0.0300183 -0.0271734
-0.0387318 -0.0542862 -0.0602707 -0.0532457
-0.0422166 -0.0388861 -0.0463080 -0.0568841
-0.0606738 -0.0550934 -0.0466514 -0.0438701
-0.0489986 -0.0565122 -0.0591451 -0.0549727
-0.0487101 -0.0465725 -0.0502234 -0.0556519
-0.0575522 -0.0544950 -0.0498990 -0.0482991
输出的单位阶跃响应图如下所示。
/*********************************************************************************************************
* 函数名称： Filter
* 函数功能： 直接性 IIR 数字滤波（一）
* 输入参数： b ：双精度实型一维数组，长度为（m+1）。存放滤波器分子多项式的系数 b(i)。
* a ：双精度实型一维数组，长度为（n+1）。存放滤波器分母多项式的系数 a(i)。
* m ：整形变量。滤波器分子多项式的阶数。
* n ：整形变量。滤波器分母多项式的阶数。
* x ：双精度实型一维数组，长度为 len。开始时存放滤波器的输入序列，最后存放滤波器的输出序列；
* 分块处理时，它用于表示当前块内的滤波器的输入序列与输出序列。
* len：整型变量。输入序列与输出序列的长度；在分块处理时，它用于表示块的长度
* px ：双精度实型一维数组，长度为（m+1）。在分块处理时，它用于保存前一块滤波时的（m+1）个输入序列值。
* py ：双精度实型一维数组，长度为（n+1）。在分块处理时，它用于保存前一块滤波时的 n 个输出序列值。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年10月01日
* 注 意： 当输入序列 x(k) 很长时，由于计算机内存的限制，常将其分成彼此衔接的若干块进行处理。数组 px 与
* py 就是专为分块处理而设置的。
* px 用于保存前一块滤波时的（m+1）个输入序列值，即 px[] = {x(k), x(k-1), ..., x(k-m)}；
* py 用于保存前一块滤波时的 n 个输出序列值。即 py[] = {y(k-1), y(k-2), ..., y(k-n)}；
* 通常，我们假定滤波器的初始条件为零，因此数组 px[] 与 py[] 在滤波前都要初始化为零。
*********************************************************************************************************/
void Filter(double* b, double* a, int m, int n, double* x, int len, double* px, double* py)
int k, i;
for (k = 0; k < len; k++)
px[0] = x[k];
x[k] = 0.0;
for (i = 0; i <= m; i++)
x[k] = x[k] + b[i] * px[i];
for (i = 1; i <= n; i++)
x[k] = x[k] - a[i] * py[i];
if (fabs(x[k]) > 1.0e10)
return; //This is an unstable filter!
for (i = m; i >= 1; i--)
px[i] = px[i - 1];
for (i = n; i >= 2; i--)
py[i] = py[i - 1];
py[1] = x[k];
4 阶切比雪夫低通数字滤波器的传递函数为
\(H(z) = \frac{0.001836 + 0.007344z^{-1} + 0.011016z^{-2} + 0.007344z^{-3} + 0.001836z^{-4}}{1.0 – 3.0544z^{-1} + 3.8291z^{-2} – 2.2925z^{-3} + 0.55075z^{-4}}\)
选取参数 m = 4，n = 4，块长 len = 25，求该滤波器的单位冲击响应，并画出相应波形。
测试代码如下：
int main(void)
int i, k, m, n, len, nblk;
static double b[5] = { 0.001836, 0.007344, 0.011016, 0.007344, 0.001836 };
static double a[5] = { 1.0, -3.0544, 3.8291, -2.2925, 0.55075 };
static double px[5] = { 0 };
static double py[5] = { 0 };
static double x[25] = {0};
static double data[100] = {0};
FILE* fp;
m = 4;
n = 4;
len = 25;
nblk = 4;
x[0] = 1.0;
for (k = 1; k < len; k++)
x[k] = 0.0;
for (i = 0; i < nblk; i++)
Filter(b, a, m, n, x, len, px, py);
for (k = 0; k < len; k++)
data[i * len + k] = x[k];
x[k] = 0.0;
printf("Unit Impulse Respomse\n");
for (i = 0; i < 16; i++)
printf(" %10.7lf", data[i]);
if (i % 4 == 3)
printf("\n");
fp = fopen("filter.csv", "w");
for (i = 0; i < 100; i++)
fprintf(fp, "%2d,%lf\n", i, data[i]);
fclose(fp);
return 0;
终端输出的单位冲激响应如下所示。
Unit Impulse Respomse
0.0018360 0.0129519 0.0435460 0.0949659
0.1538388 0.1989473 0.2123269 0.1871151
0.1298634 0.0573615 -0.0100326 -0.0556283
-0.0715163 -0.0600246 -0.0314986 0.0003172
输出的单位冲激响应如下图所示。
/*********************************************************************************************************
* 函数名称： Filtery
* 函数功能： 直接性 IIR 数字滤波（二）
* 输入参数： b ：双精度实型一维数组，长度为（m+1）。存放滤波器分子多项式的系数 b(i)。
* a ：双精度实型一维数组，长度为（n+1）。存放滤波器分母多项式的系数 a(i)。
* m ：整形变量。滤波器分子多项式的阶数。
* n ：整形变量。滤波器分母多项式的阶数。
* x ：双精度实型一维数组，长度为 len。开始时存放滤波器的输入序列；
* 分块处理时，它用于表示当前块内的滤波器的输入序列。
* y ：双精度实型一维数组，长度为 len。开始时存放滤波器的输出序列；
* 分块处理时，它用于表示当前块内的滤波器的输出序列。在滤波前必须将其初始化为零。
* len：整型变量。输入序列与输出序列的长度；在分块处理时，它用于表示块的长度
* px ：双精度实型一维数组，长度为（m+1）。在分块处理时，它用于保存前一块滤波时的（m+1）个输入序列值。
* py ：双精度实型一维数组，长度为（n+1）。在分块处理时，它用于保存前一块滤波时的 n 个输出序列值。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年10月01日
* 注 意： 当输入序列 x(k) 很长时，由于计算机内存的限制，常将其分成彼此衔接的若干块进行处理。数组 px 与
* py 就是专为分块处理而设置的。
* px 用于保存前一块滤波时的（m+1）个输入序列值，即 px[] = {x(k), x(k-1), ..., x(k-m)}；
* py 用于保存前一块滤波时的 n 个输出序列值。即 py[] = {y(k-1), y(k-2), ..., y(k-n)}；
* 通常，我们假定滤波器的初始条件为零，因此数组 px[] 与 py[] 在滤波前都要初始化为零。
*********************************************************************************************************/
void Filtery(double* b, double* a, int m, int n, double* x, double* y, int len, double* px, double* py)
int k, i;
double sum;
for (k = 0; k < len; k++)
px[0] = x[k];
sum = 0.0;
for (i = 0; i <= m; i++)
sum = sum + b[i] * px[i];
for (i = 1; i <= n; i++)
sum = sum - a[i] * py[i];
if (fabs(x[k]) > 1.0e10)
return; //This is an unstable filter!
for (i = m; i >= 1; i--)
px[i] = px[i - 1];
for (i = n; i >= 2; i--)
py[i] = py[i - 1];
py[1] = sum;
y[k] = y[k] + sum;
4 阶切比雪夫低通数字滤波器的传递函数为
\(H(z) = \frac{0.001836 + 0.007344z^{-1} + 0.011016z^{-2} + 0.007344z^{-3} + 0.001836z^{-4}}{1.0 – 3.0544z^{-1} + 3.8291z^{-2} – 2.2925z^{-3} + 0.55075z^{-4}}\)
选取参数 m = 4，n = 4，块长 len = 100，求该滤波器的单位冲击响应，并画出相应波形。
测试代码如下：
int main(void)
int i, m, n, len;
static double b[5] = { 0.001836, 0.007344, 0.011016, 0.007344, 0.001836 };
static double a[5] = { 1.0, -3.0544, 3.8291, -2.2925, 0.55075 };
static double px[5] = { 0 };
static double py[5] = { 0 };
static double x[100] = { 0 };
static double y[100] = { 0 };
FILE* fp;
m = 4;
n = 4;
len = 100;
for (i = 0; i < len; i++)
x[i] = 0.0;
y[i] = 0.0;
x[0] = 1.0;
Filtery(b, a, m, n, x, y, len, px, py);
printf("Unit Impulse Respomse\n");
for (i = 0; i < 16; i++)
printf(" %10.7lf", y[i]);
if (i % 4 == 3)
printf("\n");
fp = fopen("filtery.csv", "w");
for (i = 0; i < 100; i++)
fprintf(fp, "%2d,%lf\n", i, y[i]);
fclose(fp);
return 0;
终端输出的单位冲激响应如下所示。
Unit Impulse Respomse
0.0018360 0.0129519 0.0435460 0.0949659
0.1538388 0.1989473 0.2123269 0.1871151
0.1298634 0.0573615 -0.0100326 -0.0556283
-0.0715163 -0.0600246 -0.0314986 0.0003172
输出的单位冲激响应如下图所示。
/*********************************************************************************************************
* 函数名称： Filterc
* 函数功能： 级联型 IIR 数字滤波器
* 输入参数： b ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分子多项式的系数，
* b[j][i] 表示第 j 个 n 阶节的分子多项式的第 i 个系数。
* a ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分母多项式的系数，
* a[j][i] 表示第 j 个 n 阶节的分母多项式的第 i 个系数。
* n ：整型变量。级联型滤波器每节的阶数。
* ns ：整形变量。级联型滤波器的 n 阶节数 L。
* x ：双精度实型一维数组，长度为 len。开始时存放滤波器的输入序列，最后存放滤波器的输出序列；
* 分块处理时，它用于表示当前块内的滤波器的输入序列与输出序列。
* len：整型变量。输入序列与输出序列的长度；在分块处理时，它用于表示块的长度
* px ：双精度实型二维数组，体积为 ns*（n+1）。在分块处理时，它用于保存前一块滤波时的（n+1）个输入序列值。
* py ：双精度实型二维数组，体积为 ns*（n+1）。在分块处理时，它用于保存前一块滤波时的 n 个输出序列值。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年10月01日
* 注 意： 当输入序列 x(k) 很长时，由于计算机内存的限制，常将其分成彼此衔接的若干块进行处理。数组 px 与
* py 就是专为分块处理而设置的。
* px 用于保存前一块滤波时的（n+1）个输入序列值，即 px[j][] = {x(k), x(k-1), ..., x(k-m)}；
* py 用于保存前一块滤波时的 n 个输出序列值。即 py[j][] = {y(k-1), y(k-2), ..., y(k-n)}；
* 通常，我们假定滤波器的初始条件为零，因此数组 px[][] 与 py[][] 在滤波前都要初始化为零。
*********************************************************************************************************/
void Filterc(double* b, double* a, int n, int ns, double* x, int len, double* px, double* py)
int i, j, k, n1;
n1 = n + 1;
for (j = 0; j < ns; j++)
for (k = 0; k < len; k++)
px[j * n1 + 0] = x[k];
x[k] = b[j * n1 + 0] * px[j * n1 + 0];
for (i = 1; i <= n; i++)
x[k] += b[j * n1 + i] * px[j * n1 + i] - a[j * n1 + i] * py[j * n1 + i];
if (fabs(x[k]) > 1.0e10)
return; //This is an unstable filter!
for (i = n; i >= 2; i--)
px[j * n1 + i] = px[j * n1 + i - 1];
py[j * n1 + i] = py[j * n1 + i - 1];
px[j * n1 + 1] = px[j * n1 + 0];
py[j * n1 + 1] = x[k];
4 阶切比雪夫低通数字滤波器的传递函数为
\(H(z) = \frac{0.001836(1+z^{-1})^{4}}{(1 – 1.4996z^{-1} + 0.8482z^{-2})(1 – 1.5548z^{-1} + 0.6493z^{-2})}\)
它由两个 2 阶节级联而成。选取参数 n = 2，ns = 2，共分 4 块进行处理，每块长度 len = 25。求该滤波器的单位冲激响应，并画出相应图形。
测试代码如下：
int main(void)
int i, j, n, ns, len, nblk;
static double b[2][3] = { {0.001836, 0.003672, 0.001836}, {1.0, 2.0, 1.0} };
static double a[2][3] = { {1.0, -1.4996, 0.8482}, {1.0, -1.5548, 0.6493} };
static double px[2][3] = { 0.0 };
static double py[2][3] = { 0.0 };
static double x[25] = { 0.0 };
static double data[100] = { 0.0 };
FILE* fp;
n = 2;
ns = 2;
len = 25;
nblk = 4;
x[0] = 1.0;
for (i = 1; i < len; i++)
x[i] = 0.0;
for (i = 0; i < ns; i++)
for (j = 0; j <= n; j++)
px[i][j] = 0.0;
py[i][j] = 0.0;
for (j = 0; j < nblk; j++)
Filterc(b, a, n, ns, x, len, px, py);
for (i = 0; i < len; i++)
data[j * len + i] = x[i];
x[i] = 0.0;
printf("Unit Impulse Respomse\n");
for (i = 0; i < 16; i++)
printf(" %10.7lf", data[i]);
if (i % 4 == 3)
printf("\n");
fp = fopen("filterc.csv", "w");
for (i = 0; i < 100; i++)
fprintf(fp, "%2d,%lf\n", i, data[i]);
fclose(fp);
return 0;
终端输出的单位冲激响应如下所示。
Unit Impulse Respomse
0.0018360 0.0129519 0.0435460 0.0949662
0.1538403 0.1989518 0.2123368 0.1871326
0.1298898 0.0573960 -0.0099931 -0.0555889
-0.0714825 -0.0600010 -0.0314873 0.0003168
输出的单位冲激响应如下图所示。
/*********************************************************************************************************
* 函数名称： Filterp
* 函数功能： 并联型 IIR 数字滤波器
* 输入参数： b ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分子多项式的系数，
* b[j][i] 表示第 j 个 n 阶节的分子多项式的第 i 个系数。
* a ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分母多项式的系数，
* a[j][i] 表示第 j 个 n 阶节的分母多项式的第 i 个系数。
* n ：整型变量。级联型滤波器每节的阶数。
* ns ：整形变量。级联型滤波器的 n 阶节数 L。
* x ：双精度实型一维数组，长度为 len。存放滤波器的输入序列；
* 分块处理时，它用于表示当前块内的滤波器的输入序列。
* y ：双精度实型一维数组，长度为 len。存放滤波器的输出序列；
* 分块处理时，它用于表示当前块内的滤波器的输出序列。
* len：整型变量。输入序列与输出序列的长度；在分块处理时，它用于表示块的长度
* px ：双精度实型二维数组，体积为 ns*（n+1）。在分块处理时，它用于保存前一块滤波时的（n+1）个输入序列值。
* py ：双精度实型二维数组，体积为 ns*（n+1）。在分块处理时，它用于保存前一块滤波时的 n 个输出序列值。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年10月01日
* 注 意： 当输入序列 x(k) 很长时，由于计算机内存的限制，常将其分成彼此衔接的若干块进行处理。数组 px 与
* py 就是专为分块处理而设置的。
* px 用于保存前一块滤波时的（n+1）个输入序列值，即 px[j][] = {x(k), x(k-1), ..., x(k-m)}；
* py 用于保存前一块滤波时的 n 个输出序列值。即 py[j][] = {y(k-1), y(k-2), ..., y(k-n)}；
* 通常，我们假定滤波器的初始条件为零，因此数组 px[][] 与 py[][] 在滤波前都要初始化为零。
*********************************************************************************************************/
void Filterp(double* b, double* a, int n, int ns, double* x, double* y, int len, double* px, double* py)
int i, j, k, n1;
double sum;
n1 = n + 1;
for (k = 0; k < len; k++)
y[k] = 0.0;
for (j = 0; j < ns; j++)
for (k = 0; k < len; k++)
px[j * n1 + 0] = x[k];
sum = b[j * n1 + 0] * px[j * n1 + 0];
for (i = 1; i <= n; i++)
sum += b[j * n1 + i] * px[j * n1 + i] - a[j * n1 + i] * py[j * n1 + i];
if (fabs(sum) > 1.0e10)
return; //This is an unstable filter!
for (i = n; i >= 2; i--)
px[j * n1 + i] = px[j * n1 + i - 1];
py[j * n1 + i] = py[j * n1 + i - 1];
px[j * n1 + 1] = px[j * n1 + 0];
py[j * n1 + 1] = sum;
y[k] = y[k] + sum;
4 阶切比雪夫低通数字滤波器的传递函数为
\(H(z) = \frac{0.08372 + 0.0239z^{-1}}{1 – 1.5658z^{-1} + 0.6549z^{-2}} – \frac{0.08327 + 0.0246z^{-1}}{1 – 1.4934z^{-1} + 0.8392z^{-2}}\)
它由两个 2 阶节并联而成。选取参数 n = 2，ns = 2；共分 4 块进行处理，每块长度 len = 25。求该滤波器的单位冲击相应，并画出相应波形。
测试代码如下：
int main(void)
int i, j, n, ns, len, nblk;
static double b[2][3] = { {0.08327, 0.0239, 0.0}, {-0.08327, -0.0246, 0.0} };
static double a[2][3] = { {1.0, -1.5658, 0.6549}, {1.0, -1.4934, 0.8392} };
static double px[2][3] = { 0 };
static double py[2][3] = { 0 };
static double x[25] = { 0 };
static double y[25] = { 0 };
static double data[100] = { 0 };
FILE* fp;
n = 2;
ns = 2;
len = 25;
nblk = 4;
x[0] = 1.0;
for (i = 1; i < len; i++)
x[i] = 0.0;
for (i = 0; i < ns; i++)
for (j = 0; j <= n; j++)
px[i][j] = 0.0;
py[i][j] = 0.0;
for (j = 0; j < nblk; j++)
Filterp(b, a, n, ns, x, y, len, px, py);
for (i = 0; i < len; i++)
data[j * len + i] = y[i];
x[i] = 0.0;
printf("Unit Impulse Respomse\n");
for (i = 0; i < 16; i++)
printf(" %10.7lf", data[i]);
if (i % 4 == 3)
printf("\n");
fp = fopen("filterp.csv", "w");
for (i = 0; i < 100; i++)
fprintf(fp, "%2d,%lf\n", i, data[i]);
fclose(fp);
return 0;
Unit Impulse Respomse
0.0000000 0.0053287 0.0344748 0.0889894
0.1523266 0.2010379 0.2162496 0.1913117
0.1335624 0.0605640 -0.0069497 -0.0523818
-0.0682015 -0.0571050 -0.0295577 0.0008975
输出的单位冲激响应如下图所示。
/*********************************************************************************************************
* 函数名称： iirbcf
* 函数功能： 巴特沃兹和切比雪夫数字滤波器的设计
* 输入参数： ifilt：整形变量。滤波器的类型。取值为 1、2 和 3，分别对应切比雪夫、逆切比雪夫和巴特沃兹滤波器。
* band ：整形变量。滤波器的通带形式，取值为 1、2、3 和 4，分别对应低通、高通、带通和带阻。
* ns ：整形变量。滤波器的 n 阶节数。
* n ：整形变量。滤波器每节的阶数。对于低通和高通滤波器，n = 2；对于带通和带阻滤波器，n = 4。
* f1-f4：双精度实型变量。
* 对于巴特沃兹滤波器：
* 低通时，f1 是通带边界频率，f2 = f3 = f4 = 0；
* 高通时，f2 是通带边界频率，f1 = f3 = f4 = 0；
* 带通时，f2 是带通下边界频率，f3 是带通上边界频率，f1 = f4 = 0；
* 带阻时，f1 是带通下边界频率，f4 是带通上边界频率，f2 = f3 = 0。
* 对于切比雪夫滤波器：
* 低通时，f1 是通带边界频率，f2 是阻带边界频率，f3 = f4 = 0；
* 高通时，f2 是通带边界频率，f1 是阻带边界频率，f3 = f4 = 0；
* 带通时，f2 是通带下边界频率，f3 是通带上边界频率，f1 是阻带下边界频率，f4 是阻带上边界频率；
* 带阻时，f1 是通带下边界频率，f4 是通带上边界频率，f2 是阻带下边界频率，f3 是阻带上边界频率。
* db ：双精度实型变量。滤波器的阻带衰减（用 dB 表示）。
* b ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分子多项式的系数，
* b[j][i] 表示第 j 个 n 阶节的分子多项式的第 i 个系数。
* a ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分母多项式的系数，
* a[j][i] 表示第 j 个 n 阶节的分母多项式的第 i 个系数。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年10月01日
* 注 意：
*********************************************************************************************************/
void iirbcf(int ifilt, int band, int ns, int n, double f1, double f2, double f3, double f4, double db, double* b, double* a)
int k;
double omega, lamda, epslon, fl, fh;
double d[5], c[5];
void chebyi(int ln, int k, int n, double ep, double* d, double* c);
void chebyii(int ln, int k, int n, double ws, double att, double* d, double* c);
void bwtf(int ln, int k, int n, double* d, double* c); double cosh1(double x);
void fblt(double* d, double* c, int n, int band, double fln, double fhn, double* b, double* a);
double warp(double f); double bpsub(double om, double fh, double fl); double omin(double om1, double om2);
if ((band == 1) || (band == 4)) { fl = f1; }
if ((band == 2) || (band == 3)) { fl = f2; }
if (band <= 3) { fh = f3; }
if (band == 4) { fh = f4; }
if (ifilt < 3)
switch (band)
case 1:
case 2: omega = warp(f2) / warp(f1); break;
case 3: omega = omin(bpsub(warp(f1), fh, fl), bpsub(warp(f4), fh, fl)); break;
case 4: omega = omin(1.0 / bpsub(warp(f2), fh, fl), 1.0 / bpsub(warp(f3), fh, fl)); break;
lamda = pow(10.0, (db / 20.0));
epslon = lamda / cosh(2 * ns * cosh1(omega));
for (k = 0; k < ns; k++)
switch (ifilt)
case 1: chebyi(2 * ns, k, 4, epslon, d, c); break;
case 2: chebyii(2 * ns, k, 4, omega, lamda, d, c); break;
case 3: bwtf(2 * ns, k, 4, d, c); break;
fblt(d, c, n, band, fl, fh, &b[k * (n + 1) + 0], &a[k * (n + 1) + 0]);
static double cosh1(double x)
double z;
z = log(x + sqrt(x * x - 1.0));
return (z);
static double warp(double f)
double pi, z;
pi = 4.0 * atan(1.0);
z = tan(pi * f);
return (z);
static double bpsub(double om, double fh, double fl)
double z;
z = (om * om - warp(fh) * warp(fl)) / ((warp(fh) - warp(fl)) * om);
return (z);
static double omin(double om1, double om2)
double z, z1, z2;
z1 = fabs(om1);
z2 = fabs(om2);
z = (z1 < z2) ? z1 : z2;
return (z);
static void bwtf(int ln, int k, int n, double* d, double* c)
int i;
double pi, tmp;
pi = 4.0 * atan(1.0);
d[0] = 1.0;
c[0] = 1.0;
for (i = 1; i <= n; i++)
d[i] = 0.0;
c[i] = 0.0;
tmp = (k + 1) - (ln + 1.0) / 2.0;
if (tmp == 0.0)
c[1] = 1.0;
c[1] = -2.0 * cos((2 * (k + 1) + ln - 1) * pi / (2 * ln));
c[2] = 1.0;
static void chebyi(int ln, int k, int n, double ep, double* d, double* c)
int i;
double pi, gam, omega, sigma;
pi = 4.0 * atan(1.0);
gam = pow(((1.0 + sqrt(1.0 + ep * ep)) / ep), 1.0 / ln);
sigma = 0.5 * (1.0 / gam - gam) * sin((2 * (k + 1) - 1) * pi / (2 * ln));
omega = 0.5 * (1.0 / gam + gam) * cos((2 * (k + 1) - 1) * pi / (2 * ln));
for (i = 0; i <= n; i++)
d[i] = 0.0;
c[i] = 0.0;
if (((ln % 2) == 1) && ((k + 1) == (ln + 1) / 2))
d[0] = -sigma;
c[0] = d[0];
c[1] = 1.0;
c[0] = sigma * sigma + omega * omega;
c[1] = -2.0 * sigma;
c[2] = 1.0;
d[0] = c[0];
if (((ln % 2) == 0) && (k == 0))
d[0] = d[0] / sqrt(1.0 + ep * ep);
static void chebyii(int ln, int k, int n, double ws, double att, double* d, double* c)
int i;
double pi, gam, alpha, beta, sigma, omega, scln, scld;
pi = 4.0 * atan(1.0);
gam = pow((att + sqrt(att * att - 1.0)), 1.0 / ln);
alpha = 0.5 * (1.0 / gam - gam) * sin((2 * (k + 1) - 1) * pi / (2 * ln));
beta = 0.5 * (1.0 / gam + gam) * cos((2 * (k + 1) - 1) * pi / (2 * ln));
sigma = ws * alpha / (alpha * alpha + beta * beta);
omega = -1.0 * ws * beta / (alpha * alpha + beta * beta);
for (i = 0; i <= n; i++)
d[i] = 0.0;
c[i] = 0.0;
if (((ln % 2) == 1) && ((k + 1) == (ln + 1) / 2))
d[0] = -1.0 * sigma;
c[0] = d[0];
c[1] = 1.0;
scln = sigma * sigma + omega * omega;
scld = pow((ws / cos((2 * (k + 1) - 1) * pi / (2 * ln))), 2);
d[0] = scln * scld;
d[2] = scln;
c[0] = d[0];
c[1] = -2.0 * sigma * scld;
c[2] = scld;
static void fblt(double* d, double* c, int n, int band, double fln, double fhn, double* b, double* a)
int i, k, m, n1, n2, ls;
double pi, w, w0, w1, w2, tmp, tmpd, tmpc, *work;
double combin(int i1, int i2);
void bilinear(double* d, double* c, double* b, double* a, int n);
pi = 4.0 * atan(1.0);
w1 = tan(pi * fln);
for (i = n; i >= 0; i--)
if ((c[i] != 0.0) || (d[i] != 0.0))
break;
m = i;
switch (band)
case 1:
case 2:
n2 = m;
n1 = n2 + 1;
if (band == 2)
for (i = 0; i <= m / 2; i++)
tmp = d[i];
d[i] = d[m - i];
d[m - i] = tmp;
tmp = c[i];
c[i] = c[m - i];
c[m - i] = tmp;
for (i = 0; i <= m; i++)
d[i] = d[i] / pow(w1, i);
c[i] = c[i] / pow(w1, i);
break;
case 3:
case 4:
n2 = 2 * m;
n1 = n2 + 1;
work = malloc(n1 * n1 * sizeof(double));
w2 = tan(pi * fhn);
w = w2 - w1;
w0 = w1 * w2;
if (band == 4)
for (i = 0; i <= m / 2; i++)
tmp = d[i];
d[i] = d[m - i];
d[m - i] = tmp;
tmp = c[i];
c[i] = c[m - i];
c[m - i] = tmp;
for (i = 0; i <= n2; i++)
work[0 * n1 + i] = 0.0;
work[1 * n1 + i] = 0.0;
for (i = 0; i <= m; i++)
tmpd = d[i] * pow(w, (m - i));
tmpc = c[i] * pow(w, (m - i));
for (k = 0; k <= i; k++)
ls = m + i - 2 * k;
tmp = combin(i, i) / (combin(k, k) * combin(i - k, i - k));
work[0 * n1 + ls] += tmpd * pow(w0, k) * tmp;
work[1 * n1 + ls] += tmpc * pow(w0, k) * tmp;
for (i = 0; i <= n2; i++)
d[i] = work[0 * n1 + i];
c[i] = work[1 * n1 + i];
free(work);
bilinear(d, c, b, a, n);
static double combin(int i1, int i2)
int i;
double s;
s = 1.0;
if (i2 == 0)
return (s);
for (i = i1; i > (i1 - i2); i--)
s *= i;
return (s);
static void bilinear(double* d, double* c, double* b, double* a, int n)
int i, j, n1;
double sum, atmp, scale, *temp;
n1 = n + 1;
temp = malloc(n1 * n1 * sizeof(double));
for (j = 0; j <= n; j++)
temp[j * n1 + 0] = 1.0;
sum = 1.0;
for (i = 1; i <= n; i++)
sum = sum * (double)(n - i + 1) / (double)i;
temp[0 * n1 + i] = sum;
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
temp[j * n1 + i] = temp[(j - 1) * n1 + i] - temp[j * n1 + i - 1] - temp[(j - 1) * n1 + i - 1];
for (i = n; i >= 0; i--)
b[i] = 0.0;
atmp = 0.0;
for (j = 0; j <= n; j++)
b[i] = b[i] + temp[j * n1 + i] * d[j];
atmp = atmp + temp[j * n1 + i] * c[j];
scale = atmp;
if (i != 0)
a[i] = atmp;
for (i = 0; i <= n; i++)
b[i] = b[i] / scale;
a[i] = a[i] / scale;
a[0] = 1.0;
free(temp);
下面给出主函数程序，它调用 iirbcf.c 子函数。通过人机对话方式输入参数后，它可以设计巴特沃兹、切比雪夫和逆切比雪夫滤波器，每种滤波器都具有低通、高通、带通、带阻这 4 种形式。下面对输入参数进行说明：
ifilt：滤波器的类型；band：滤波器的通带；ns：滤波器的 n 阶节数；fs：采样频率；db：阻带衰减；fname：幅频响应的文件名。
对于低通和高通滤波器，fc：通带边界频率；fr：阻带边界频率。
对于带通和带阻滤波器，flc：通带下边界频率；fhc：通带上边界频率；fls：阻带下边界频率；fhs：阻带上边界频率。
测试代码如下：
int main(void)
int i, k, n, ns, ifilt, band;
static double a[50] = { 0 };
static double b[50] = { 0 };
static double x[300] = { 0 };
static double y[300] = { 0 };
double f1, f2, f3, f4, fc, fr, fs, flc, fls, fhc, fhs, freq, db;
static char fname[40] = {0};
FILE* fp;
db = 1;
printf("Enter 1 for Chebyshev I, 2 for Chebyshev II, 3 for Butterworth\n");
scanf("%d", &ifilt);
printf("Enter 1 for lowpass, 2 for highpass, 3 for bandpass, 4 for bandstop\n");
scanf("%d", &band);
n = (band <= 2) ? 2 : 4;
printf("Enter the number of filter section ns\n");
scanf("%d", &ns);
printf("Enter sample frequence fs\n");
scanf("%lf", &fs);
if (ifilt <= 2)
switch (band)
case 1:
case 2:
printf("Enter passband edge frequency fc\n");
scanf("%lf", &fc);
printf("Enter stopband edge frequency fr\n");
scanf("%lf", &fr);
if (band == 1)
f1 = fc;
f2 = fr;
f1 = fr;
f2 = fc;
f3 = f4 = 0.0;
break;
case 3:
case 4:
printf("Enter the lower passband edge frequency flc\n");
scanf("%lf", &flc);
printf("Enter the higher passband edge frequency fhc\n");
scanf("%lf", &fhc);
printf("Enter the lower stopband edge frequency fls\n");
scanf("%lf", &fls);
printf("Enter the higher stopband edge frequency fhs\n");
scanf("%lf", &fhs);
if (band == 3)
f1 = fls;
f2 = flc;
f3 = fhc;
f4 = fhs;
f1 = flc;
f2 = fls;
f3 = fhs;
f4 = fhc;
printf("Enter stopband attenuation (dB)\n");
scanf("%lf", &db);
switch (band)
case 1:
case 2:
printf("Enter passband edge frequency fc\n");
scanf("%lf", &fc);
f1 = f2 = f3 = f4 = 0.0;
if (band == 1)
f1 = fc;
f2 = fc;
break;
case 3:
case 4:
printf("Enter the lower passband edge frequency flc\n");
scanf("%lf", &flc);
printf("Enter the higher passband edge frequency fhc\n");
scanf("%lf", &fhc);
f1 = f2 = f3 = f4 = 0.0;
if (band == 3)
f2 = flc;
f3 = fhc;
f1 = flc;
f4 = fhc;
f1 = f1 / fs;
f2 = f2 / fs;
f3 = f3 / fs;
f4 = f4 / fs;
iirbcf(ifilt, band, ns, n, f1, f2, f3, f4, db, b, a);
for (k = 0; k < ns; k++)
printf("\nSection %d\n\n", k + 1);
for (i = 0; i <= n; i++)
printf(" b[%d][%d] = %10.7lf", k, i, b[k * (n + 1) + i]);
if (((i % 2) == 0) && (i != 0))
printf("\n");
printf("\n");
for (i = 0; i <= n; i++)
printf(" a[%d][%d] = %10.7lf", k, i, a[k * (n + 1) + i]);
if (((i % 2) == 0) && (i != 0))
printf("\n");
printf("\n");
printf("\nEnter file name of magnitude response\n");
scanf("%s", fname);
fp = fopen(fname, "w");
Gainc(b, a, n, ns, x, y, 300, 2);
for (i = 0; i < 300; i++)
freq = i * 0.5 / 300.0;
fprintf(fp, "%lf,%lf\n", freq, x[i]);
fclose(fp);
return 0;
使用 IIR 滤波器滤除心电信号工频干扰
首先，使用上一节提供的人机接口程序设计 IIR 滤波器。心电信号采样频率为 2kHz，在文章的末尾给出。使用 4 级巴特沃兹低通数字滤波器，通带边界频率为 25Hz。即输入参数 ifilt = 3，band = 1，ns = 4， fc = 25，fs = 2000。将幅频特性曲线输出到 result.csv。终端的输入输出如下所示。
Enter 1 for Chebyshev I, 2 for Chebyshev II, 3 for Butterworth
Enter 1 for lowpass, 2 for highpass, 3 for bandpass, 4 for bandstop
Enter the number of filter section ns
Enter sample frequence fs
Enter passband edge frequency fc
Section 1
b[0][0] = 0.0015181 b[0][1] = 0.0030362 b[0][2] = 0.0015181
a[0][0] = 1.0000000 a[0][1] = -1.9637759 a[0][2] = 0.9698483
Section 2
b[1][0] = 0.0014770 b[1][1] = 0.0029539 b[1][2] = 0.0014770
a[1][0] = 1.0000000 a[1][1] = -1.9105545 a[1][2] = 0.9164623
Section 3
b[2][0] = 0.0014469 b[2][1] = 0.0028939 b[2][2] = 0.0014469
a[2][0] = 1.0000000 a[2][1] = -1.8717298 a[2][2] = 0.8775176
Section 4
b[3][0] = 0.0014312 b[3][1] = 0.0028624 b[3][2] = 0.0014312
a[3][0] = 1.0000000 a[3][1] = -1.8513690 a[3][2] = 0.8570938
Enter file name of magnitude response
result.csv
输出的滤波器系数为 2 阶级联型的，将滤波器系数保存到数组后，效果如下所示。
static const double s_arrB[4][3] = {{0.0015181, 0.0030362, 0.0015181}, {0.0014770, 0.0029539, 0.0014770},
{0.0014469, 0.0028939, 0.0014469}, {0.0014312, 0.0028624, 0.0014312}};
static const double s_arrA[4][3] = {{1.0000000, -1.9637759, 0.9698483}, {1.0000000, -1.9105545, 0.9164623},
{1.0000000, -1.8717298, 0.8775176}, {1.0000000, -1.8513690, 0.8570938}};
滤波器的幅频特性曲线如下所示。纵坐标的单位是 dB，横坐标是归一化频率，乘以采样率，即 2000Hz，即可转换到频率。滤波器在 0.013333，即 26.67Hz 处的增益为 -5.983596dB，小于 -3dB。
extern const double g_arrEcgWave[4000];
static double s_arrInput[4000] = {0};
static double s_arrOutput[4000] = {0};
static const double s_arrB[4][3] = {{0.0015181, 0.0030362, 0.0015181}, {0.0014770, 0.0029539, 0.0014770},
{0.0014469, 0.0028939, 0.0014469}, {0.0014312, 0.0028624, 0.0014312}};
static const double s_arrA[4][3] = {{1.0000000, -1.9637759, 0.9698483}, {1.0000000, -1.9105545, 0.9164623},
{1.0000000, -1.8717298, 0.8775176}, {1.0000000, -1.8513690, 0.8570938}};
static double s_arrPx[12] = {0};
static double s_arrPy[12] = {0};
int i, len;
long seed;
double pi, time;
FILE* fp;
//定义数据量
len = 4000;
//获取输入数据，加 50Hz 工频干扰
seed = 13579l;
pi = 3.1415926535;
time = 0;
for (i = 0; i < len; i++)
s_arrInput[i] = g_arrEcgWave[i] + 0.1 * sin(time * 2 * pi / (1.0 / 50.0));
time = time + 0.0005;
s_arrOutput[i] = s_arrInput[i];
//IIR 滤波
Filterc(s_arrB, s_arrA, 2, 4, s_arrInput, len, s_arrOutput, s_arrPx, s_arrPy);
//输出波形数据，按照原始信号、参考信号、滤波后信号顺序
fp = fopen("IIRFilter.csv", "w");
for (i = 0; i < len; i++)
fprintf(fp, "%d,%10.7lf,%10.7lf,%10.7lf\n", i, s_arrInput[i], g_arrEcgWave[i], s_arrOutput[i]);
return 0;
最终实验结果如下所示。蓝色为参考信号，黄色为输入信号，绿色为输出信号。
/*********************************************************************************************************
* 函数名称： IIRFilterc
* 函数功能： 级联型 IIR 数字滤波器
* 输入参数： b ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分子多项式的系数，
* b[j][i] 表示第 j 个 n 阶节的分子多项式的第 i 个系数。
* a ：双精度实型二维数组，体积为 ns*（n+1）。存放滤波器分母多项式的系数，
* a[j][i] 表示第 j 个 n 阶节的分母多项式的第 i 个系数。
* n ：整型变量。级联型滤波器每节的阶数。
* ns ：整形变量。级联型滤波器的 n 阶节数 L。
* px ：双精度实型二维数组，体积为 ns*（n+1）。存放历史输入数据，初始状态必须为零。
* py ：双精度实型二维数组，体积为 ns*（n+1）。存放历史输入数据，初始状态必须为零。
* x ：双精度实型变量。新的采样值；
* 输出参数： void
* 返 回 值： 滤波后的采样值
* 创建日期： 2023年10月10日
* 注 意：
*********************************************************************************************************/
double IIRFilterc(double* b, double* a, int n, int ns, double* px, double* py, double x)
int i, j, n1;
n1 = n + 1;
for (j = 0; j < ns; j++)
px[j * n1 + 0] = x;
x = b[j * n1 + 0] * px[j * n1 + 0];
for (i = 1; i <= n; i++)
x += b[j * n1 + i] * px[j * n1 + i] - a[j * n1 + i] * py[j * n1 + i];
for (i = n; i >= 2; i--)
px[j * n1 + i] = px[j * n1 + i - 1];
py[j * n1 + i] = py[j * n1 + i - 1];
px[j * n1 + 1] = px[j * n1 + 0];
py[j * n1 + 1] = x;
return x;
应用示例代码如下，同样可以直接滤除掉心电信号中的 50Hz 工频干扰。
int main(void)
extern const double g_arrEcgWave[4000];
static double s_arrInput[4000] = {0};
static double s_arrOutput[4000] = {0};
static const double s_arrB[4][3] = {{0.0015181, 0.0030362, 0.0015181}, {0.0014770, 0.0029539, 0.0014770},
{0.0014469, 0.0028939, 0.0014469}, {0.0014312, 0.0028624, 0.0014312}};
static const double s_arrA[4][3] = {{1.0000000, -1.9637759, 0.9698483}, {1.0000000, -1.9105545, 0.9164623},
{1.0000000, -1.8717298, 0.8775176}, {1.0000000, -1.8513690, 0.8570938}};
static double s_arrPx[12] = {0};
static double s_arrPy[12] = {0};
int i, len;
long seed;
double pi, time;
FILE* fp;
//定义数据量
len = 4000;
//获取输入数据，加 50Hz 工频干扰
seed = 13579l;
pi = 3.1415926535;
time = 0;
for (i = 0; i < len; i++)
s_arrInput[i] = g_arrEcgWave[i] + 0.1 * sin(time * 2 * pi / (1.0 / 50.0));
time = time + 0.0005;
//IIR 滤波
for (i = 0; i < len; i++)
s_arrOutput[i] = IIRFilterc(s_arrB, s_arrA, 2, 4, s_arrPx, s_arrPy, s_arrInput[i]);
//输出波形数据，按照原始信号、参考信号、滤波后信号顺序
fp = fopen("IIRFilter.csv", "w");
for (i = 0; i < len; i++)
fprintf(fp, "%d,%10.7lf,%10.7lf,%10.7lf\n", i, s_arrInput[i], g_arrEcgWave[i], s_arrOutput[i]);
return 0;
本文中使用到的心电数据如下所示。
* 心电波形数据
* 采样率 2kHz
* 脉率值 60BPM
* 单位 mV
const double g_arrEcgWave[4000] =
-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,
-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.05,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,
-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,
-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,
-0.05,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.05,-0.06,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,
-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,
-0.06,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,
-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,
-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.06,-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,
-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.04,-0.04,-0.04,-0.04,-0.03,-0.03,-0.02,-0.02,-0.03,
-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.01,-0.01,-0.01,-0.01,0,0,0,0,0,0,0,0,0.01,0.01,0.01,0.01,0.01,0.01,
0.01,0.01,0.01,0.01,0.01,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.03,0.03,0.02,0.03,0.03,0.03,0.03,0.03,
0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,
0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.02,0.02,0.03,0.03,0.03,0.03,0.02,
0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.02,0.02,0.02,0.02,0.01,0.01,0.01,0.01,0.01,0.01,0,0,0.01,0,0,0.01,0,0,0,
0,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.02,-0.02,-0.02,-0.02,-0.03,-0.03,-0.02,-0.03,-0.03,-0.03,-0.03,-0.04,
-0.04,-0.04,-0.04,-0.04,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,
-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,
-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.06,-0.07,
-0.08,-0.08,-0.09,-0.1,-0.1,-0.1,-0.11,-0.12,-0.12,-0.12,-0.12,-0.14,-0.14,-0.14,-0.15,-0.16,-0.17,-0.17,-0.17,
-0.19,-0.19,-0.19,-0.17,-0.13,-0.12,-0.12,-0.09,-0.05,-0.05,-0.05,-0.04,0,0.02,0.02,0.03,0.08,0.09,0.09,0.11,0.15,
0.16,0.16,0.19,0.23,0.23,0.23,0.26,0.3,0.3,0.3,0.31,0.35,0.37,0.37,0.39,0.43,0.44,0.44,0.45,0.49,0.51,0.51,0.53,
0.57,0.58,0.58,0.61,0.64,0.64,0.64,0.65,0.64,0.64,0.64,0.65,0.65,0.65,0.64,0.63,0.58,0.57,0.56,0.53,0.49,0.48,0.49,
0.47,0.42,0.4,0.4,0.38,0.34,0.32,0.32,0.3,0.26,0.24,0.24,0.22,0.17,0.16,0.15,0.12,0.08,0.07,0.07,0.04,0,-0.01,-0.01,
-0.02,-0.07,-0.09,-0.09,-0.11,-0.16,-0.17,-0.17,-0.18,-0.24,-0.25,-0.25,-0.28,-0.33,-0.33,-0.33,-0.32,-0.29,-0.28,
-0.29,-0.28,-0.25,-0.23,-0.24,-0.23,-0.2,-0.2,-0.19,-0.18,-0.16,-0.15,-0.14,-0.13,-0.1,-0.1,-0.1,-0.08,-0.05,-0.06,
-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,
-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,
-0.05,-0.05,-0.05,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.05,-0.05,-0.05,-0.04,-0.04,
-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,
-0.04,-0.04,-0.04,-0.04,-0.04,-0.03,-0.03,-0.04,-0.04,-0.04,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,
-0.03,-0.03,-0.03,-0.02,-0.02,-0.02,-0.03,-0.03,-0.03,-0.02,-0.02,-0.02,-0.03,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,
-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.01,-0.01,-0.02,-0.02,-0.02,-0.01,-0.01,-0.02,-0.02,-0.01,
-0.02,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,0,-0.01,-0.01,-0.01,-0.01,-0.01,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.02,0.02,0.02,0.02,
0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.03,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.03,0.03,0.03,0.03,
0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.03,0.04,0.04,
0.04,0.04,0.05,0.05,0.05,0.04,0.04,0.04,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,
0.05,0.06,0.06,0.06,0.05,0.06,0.06,0.06,0.05,0.05,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,
0.06,0.06,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.08,0.07,0.07,0.07,0.07,0.08,0.08,0.08,0.08,0.08,
0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.09,0.09,0.08,0.09,0.08,0.08,0.08,0.08,
0.08,0.08,0.08,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.1,0.1,0.09,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.11,0.11,0.11,
0.11,0.11,0.11,0.11,0.11,0.12,0.12,0.12,0.12,0.12,0.13,0.13,0.13,0.13,0.13,0.14,0.14,0.14,0.14,0.14,0.14,0.14,
0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.15,0.16,0.16,0.16,0.16,0.16,0.16,0.16,
0.16,0.16,0.16,0.16,0.16,0.16,0.17,0.17,0.17,0.17,0.16,0.16,0.16,0.16,0.16,0.16,0.16,0.16,0.15,0.15,0.15,0.15,
0.15,0.15,0.15,0.15,0.14,0.14,0.14,0.14,0.13,0.14,0.13,0.13,0.13,0.13,0.13,0.13,0.12,0.12,0.12,0.12,0.12,0.12,
0.12,0.12,0.12,0.11,0.12,0.12,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.11,
0.1,0.1,0.1,0.1,0.1,0.09,0.1,0.09,0.08,0.08,0.08,0.08,0.07,0.07,0.06,0.06,0.06,0.05,0.05,0.05,0.04,0.04,0.04,
0.03,0.03,0.03,0.03,0.02,0.01,0.01,0.01,0.01,0,0,0,0,-0.01,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,
-0.03,-0.03,-0.03,-0.03,-0.03,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.05,-0.05,-0.04,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.05,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,
-0.06,-0.06,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,
-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,
-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.05,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.06,
-0.06,-0.06,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,
-0.06,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,-0.06,-0.06,-0.06,-0.05,-0.06,-0.06,
-0.06,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,
-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.06,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.04,-0.04,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.04,-0.04,-0.04,-0.04,-0.03,-0.03,
-0.03,-0.03,-0.03,-0.03,-0.02,-0.03,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.01,-0.01,-0.01,-0.01,0,0,0,0,0.01,0.01,
0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,
0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.04,0.04,0.04,0.03,0.03,0.03,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.03,0.03,0.04,
0.04,0.04,0.04,0.04,0.04,0.03,0.03,0.04,0.04,0.04,0.04,0.04,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.04,0.03,0.03,0.03,
0.03,0.03,0.03,0.03,0.03,0.02,0.02,0.02,0.03,0.02,0.02,0.03,0.03,0.02,0.02,0.02,0.02,0.02,0.01,0.02,0.01,0.01,0.01,
0.01,0.01,0.01,0.01,0.01,0.01,0,0,0,0,0,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.02,-0.02,-0.02,-0.02,-0.02,
-0.02,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.04,-0.04,-0.04,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.04,-0.04,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.04,-0.05,-0.04,-0.05,-0.05,-0.05,-0.04,
-0.04,-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.07,-0.07,-0.07,-0.08,-0.09,-0.09,
-0.1,-0.1,-0.11,-0.12,-0.12,-0.12,-0.13,-0.14,-0.14,-0.14,-0.15,-0.16,-0.17,-0.17,-0.19,-0.19,-0.19,-0.17,-0.13,
-0.12,-0.12,-0.11,-0.06,-0.04,-0.04,-0.03,0,0.02,0.02,0.03,0.06,0.09,0.1,0.11,0.15,0.16,0.16,0.18,0.22,0.23,0.23,
0.24,0.28,0.3,0.3,0.31,0.35,0.37,0.37,0.38,0.42,0.44,0.44,0.46,0.5,0.51,0.51,0.54,0.57,0.58,0.58,0.59,0.63,0.65,
0.65,0.65,0.65,0.65,0.65,0.65,0.65,0.65,0.65,0.63,0.58,0.57,0.57,0.54,0.5,0.48,0.49,0.48,0.44,0.41,0.4,0.39,0.35,
0.33,0.33,0.32,0.27,0.24,0.24,0.22,0.18,0.16,0.16,0.14,0.09,0.08,0.08,0.07,0.02,-0.01,0,-0.02,-0.06,-0.09,-0.09,
-0.1,-0.14,-0.17,-0.17,-0.19,-0.24,-0.25,-0.25,-0.28,-0.33,-0.33,-0.33,-0.32,-0.3,-0.29,-0.28,-0.28,-0.25,-0.24,
-0.24,-0.24,-0.21,-0.19,-0.19,-0.18,-0.15,-0.14,-0.14,-0.13,-0.1,-0.1,-0.09,-0.09,-0.06,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.04,-0.04,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.04,-0.05,-0.05,-0.05,-0.05,-0.04,-0.04,-0.05,-0.05,-0.05,-0.05,-0.04,-0.04,-0.04,-0.04,
-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.04,-0.03,-0.03,-0.04,-0.04,
-0.04,-0.04,-0.04,-0.03,-0.03,-0.04,-0.04,-0.04,-0.03,-0.03,-0.03,-0.03,-0.03,-0.04,-0.04,-0.03,-0.04,-0.04,-0.03,
-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.02,-0.02,-0.02,-0.03,
-0.03,-0.02,-0.03,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,
-0.02,-0.02,-0.02,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,-0.01,0,0,0,
-0.01,0,0,0,-0.01,0,0,0,0,0,0,0,0,0,0.01,0.01,0,0,0,0.01,0,0.01,0.01,0.01,0,0.01,0.01,0.01,0.01,0.01,0.01,0.01,
0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.02,
0.02,0.02,0.02,0.02,0.02,0.02,0.02,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.03,0.04,0.04,0.04,
0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.04,0.05,0.05,0.05,0.04,
0.04,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.06,0.05,0.05,0.05,0.06,0.06,0.06,0.06,0.06,0.06,0.05,0.06,
0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.07,0.07,0.06,0.07,0.06,0.06,0.07,0.07,0.07,0.07,
0.07,0.07,0.08,0.07,0.07,0.07,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.08,0.09,0.08,0.08,0.08,0.08,0.09,0.08,0.08,
0.08,0.09,0.09,0.09,0.08,0.08,0.09,0.08,0.08,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.09,0.1,
0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.11,0.11,0.11,0.11,0.11,0.11,0.11,0.12,0.12,0.12,0.13,0.12,
0.12,0.13,0.13,0.13,0.13,0.13,0.13,0.13,0.13,0.13,0.14,0.14,0.14,0.14,0.15,0.15,0.15,0.15,0.15,0.14,0.15,0.15,
0.16,0.16,0.16,0.15,0.15,0.16,0.16,0.16,0.16,0.16,0.16,0.16,0.16,0.16,0.16,0.17,0.17,0.17,0.17,0.17,0.16,0.16,
0.17,0.17,0.17,0.17,0.17,0.17,0.17,0.17,0.16,0.16,0.16,0.16,0.16,0.15,0.16,0.15,0.15,0.15,0.15,0.15,0.15,0.14,
0.14,0.14,0.14,0.14,0.13,0.13,0.14,0.14,0.13,0.13,0.13,0.13,0.13,0.12,0.12,0.12,0.12,0.12,0.12,0.12,0.12,0.12,
0.12,0.12,0.12,0.12,0.12,0.12,0.12,0.11,0.11,0.12,0.11,0.11,0.11,0.11,0.11,0.11,0.1,0.1,0.11,0.1,0.11,0.11,0.1,
0.1,0.1,0.09,0.09,0.09,0.08,0.07,0.07,0.07,0.07,0.06,0.05,0.06,0.06,0.04,0.04,0.04,0.04,0.04,0.03,0.03,0.03,0.02,
0.02,0.01,0.01,0.01,0.01,0.01,0,0,-0.01,-0.01,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.02,-0.03,-0.03,-0.03,-0.03,
-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.03,-0.04,-0.04,-0.04,-0.05,-0.05,-0.04,-0.04,-0.04,-0.05,-0.05,-0.05,-0.06,
-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,
-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,
-0.05,-0.06,-0.05,-0.05,-0.06,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.06,
-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,
-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.06,-0.05,-0.06,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.06,-0.06,-0.06,-0.05,-0.05,-0.05,-0.06,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.06,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,-0.05,
统计
- 1
- 113
- 49,189
归档
2023年10月
2023年9月
2023年8月
2023年7月
2023年5月
2023年4月
2023年3月
2023年2月
2022年12月
推荐文章
近期评论