Horner 算法是以英国数学家 William George Horner 命名的一种多项式求值的快速算法,实际上,这种快速算法在他之前就已经被Paolo Ruffini使用过了。而中国数学家秦九韶提出这种算法要比William George Horner 早600多年。
这个代码用到了C99中对复数支持的新特性,需要编译器支持C99标准。
将输出结果画出了可以得到如下的幅频特性:
P(x) 是一个多项式:
我们希望计算x取某个特殊值x0时多项式的值p(x0).
构造一个序列:
那么这个序列b0的值就是多项式的值了。
用程序实现如下:
double horner(double p[], int n, double x)
{
double sum;
sum = p[--n];
while ( n > 0 )
{
sum = p[--n] + sum * x;
}
return sum;
}
我经常要设计FIR、IIR 滤波器,设计完滤波器后都要验证一下频率特性是否正确。这时就需要计算实系数多项式在x取复数值时的结果。因此就有了下面的代码:
double _Complex horner_C(double p[], int n, double _Complex x)
{
double _Complex sum;
sum = p[--n];
while ( n > 0 )
{
sum = p[--n] + sum * x;
}
return sum;
}
这个代码用到了C99中对复数支持的新特性,需要编译器支持C99标准。
最后给一个测试代码,其中多项式p构成了一个FIR低通滤波器,这个低通滤波器是我用scilab 中的ffilt 函数设计的。
滤波器构造的代码如下:
ffilt('lp',20 , 0.2 , 0.5);
下面是计算这个滤波器的频响特性的代码。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
double horner(double p[], int n, double x);
double _Complex horner_C(double p[], int n, double _Complex x);
int main()
{
int i;
double p[] = {
-0.0196945, -0.0356154,
1.559E-17, 0.0465740,
0.0340178, -0.0415773,
-0.0864945, 1.559E-17,
0.2018205, 0.3741957,
0.3741957, 0.2018205,
1.559E-17, -0.0864945,
-0.0415773, 0.0340178,
0.0465740, 1.559E-17,
-0.0356154, -0.0196945
};
double x[201], mag[201], pha[201];
double _Complex xx;
for(i = 0; i <= 200; i++)
{
x[i] = 0.5 * i / 200.0;
xx = cexp(-2 * M_PI * x[i] * I);
mag[i] = cabs(horner_C(p, 20, xx));
pha[i] = carg(horner_C(p, 20, xx));
printf("%f, %f, %f\n", x[i] , mag[i], pha[i]);
}
return 0;
}
将输出结果画出了可以得到如下的幅频特性:
scilab 计算出的幅频特性如下图:
可以看出两幅图完全相同,说明我们的代码计算结果是正确的。