C语言计算数字滤波器的幅频响应和相频响应

数字信号处理C语言程序集，P168

方法简介：

子函数语句：

void gain(double b[], double a[], int m, int n, double x[], double y[], int len, int sign)

系数介绍：

b：长度为m+1，存放滤波器分子多项式的系数b（i）
a：长度为n+1，存放滤波器分母多项式的系数a（i)
m: 滤波器分子多项式的阶数
n: 滤波器分母多项式的阶数
x: 长度为len,sign为0时，存放滤波器频率响应的实部Re[H(w)],当sign=1时，存放滤波器幅度响应|H（w）|；当sign=2时，存放用分贝表示的滤波器幅频响应|H（w）|
y：长度为len，当sign=0时，存放滤波器频率响应的虚部Im[H(w)],当sign=1和2时，存放滤波器的相频响应
len：频率响应的长度
sign：当sign=0时，计算滤波器频率响应的实部Re[H(w)]和虚部Im[H（w)],当sign=1时，计算滤波器的幅频响应|H(w)|和相频响应；当sign=2时，计算滤波器的幅频响应|H（w）|和相频响应。

子程序代码如下：

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);
         }
        }
    }
}

例如：

数字系统的传递函数为：
求该系统的幅频响应和相频响应，并画出相应的图形。

我是利用Qt5.9实现的。列举主要代码。绘图利用的是QCustomPlot。

double a[] = {1.0, 0.0, 0.9};
double b[] = {0.0, -0.1};
double x[300];
double y[300];
gain(b, a, 1, 2, x, y, 300, 1);

数字系统的幅频响应：

数字系统的相频响应：