绘制数字滤波器的频域响应，对比C语言与MATLAB的结果

设传递函数$H(z)$的分子分母系数为：

b=[0.0563   -0.0009   -0.0009    0.0563];%分子
a=[1.0000   -2.1291    1.7834   -0.5435];%分母

MATLAB代码

b=[0.0563   -0.0009   -0.0009    0.0563];%分子
a=[1.0000   -2.1291    1.7834   -0.5435];%分母
for k=1:2002
    w(k)=2*pi*(k-1)/2001;
    z=exp(j*w(k));
    Z=[1 z^-1 z^-2 z^-3];%多项式
    H_z(k)=sum(b.*Z)./sum(a.*Z);%传递函数H(z)
end
figure(1);plot(w/pi,20*log10(abs(H_z)));
figure(2);plot(w/pi,angle(abs(H_z)));

【注】MATLAB代码可以更简单，用freqz(b,a,2001)即可画出幅频、相频特性。但为了与C代码原理上更接近，这里直接将$z=\exp (j\omega )$带入$H(z)$计算。

生成的幅频、相频特性曲线如下：

C代码（H_zdatagen.c）：复数要对实部虚部分别计算。特别注意，多项式用了循环，可扩展到更高阶。计算相角时，要将反正切的结果（$-\pi /2,\pi /2$）扩展到$-\pi$到$\pi$。

#include<stdio.h>
#include<math.h>
#define pi 3.14159265
int main()
{
  double b[]= {0.0563,-0.0009,-0.0009,0.0563};//分子系数
  double a[]= {1.0000,-2.1291,1.7834,-0.5435};//分母系数
  double H_z[2],fenzi[2]= {0,0},fenmu[2]= {0,0};
  double HzAmp,HzAng,tmp;//传递函数幅频、相频
  double w;//角频率
  int i,j;
  for(i=0; i<2001; i++)
    {
      w=2*pi*i/2000.0;
      fenzi[0]=0;//分子实部
      fenzi[1]=0;//分子虚部
      fenmu[0]=0;//分母实部
      fenmu[1]=0;//分母虚部
      for(j=0; j<4; j++) //计算传递函数的分子分母
        {
          fenzi[0]=fenzi[0]+b[j]*cos(-j*w);//real
          fenzi[1]=fenzi[1]+b[j]*sin(-j*w);//imag
          fenmu[0]=fenmu[0]+a[j]*cos(-j*w);//real
          fenmu[1]=fenmu[1]+a[j]*sin(-j*w);//imag
        }
	  //计算分子除以分母
	  tmp=fenmu[0]*fenmu[0]+fenmu[1]*fenmu[1];
      H_z[0]=(fenzi[0]*fenmu[0]+fenzi[1]*fenmu[1])/tmp;
      H_z[1]=(fenzi[1]*fenmu[0]-fenzi[0]*fenmu[1])/tmp;
      //计算幅频、相频响应数据
      HzAmp=sqrt(H_z[0]*H_z[0]+H_z[1]*H_z[1]);
      //相角在-pi到pi
      HzAng=atan(H_z[1]/H_z[0]);
      if(H_z[1]>0 & H_z[0]<0) HzAng=HzAng+pi;
      if(H_z[1]<0 & H_z[0]<0) HzAng=HzAng-pi;
      printf("%e\t%e\t%e\n",w/pi,20*log10(HzAmp),HzAng);
    }
}

编译生成a.exe并执行：输出三列为：归一化角频率、幅度、相角

C:\Users\邵玉斌\Desktop\123>gcc H_zdatagen.c
C:\Users\邵玉斌\Desktop\123>a.exe

0.000000e+000   7.714620e-015   0.000000e+000
1.000000e-003   -1.838922e-004  -1.017884e-002
2.000000e-003   -7.354295e-004  -2.035681e-002
3.000000e-003   -1.654195e-003  -3.053304e-002
4.000000e-003   -2.939493e-003  -4.070666e-002
5.000000e-003   -4.590353e-003  -5.087681e-002
6.000000e-003   -6.605527e-003  -6.104262e-002
7.000000e-003   -8.983494e-003  -7.120324e-002
...共2001行。

用gnuplot作图：

幅度频率响应：

gnuplot> plot [0:2] [-90:10] "<a.exe" u 1:2 w l

相频：

gnuplot> plot [0:2] [-4:4] "<a.exe" u 1:3 w l

从图上结果与MATLAB的一致。

【注释】在C代码中，反正切和相角计算的三句可用函数atan2(虚部，实部)实现。

即把

HzAng=atan(H_z[1]/H_z[0]);
if(H_z[1]>0 & H_z[0]<0) HzAng=HzAng+pi;
if(H_z[1]<0 & H_z[0]<0) HzAng=HzAng-pi;

写为一句：

HzAng=atan2(H_z[1],H_z[0]);

利用C99中新定义的复数类型（包含complex.h,在）可简化计算代码。改造后代码更简单。

#include<stdio.h>
#include<math.h>
#include<complex.h>
#define pi 3.14159265
int main()
{
  double b[]= {0.0563,-0.0009,-0.0009,0.0563};//分子系数
  double a[]= {1.0000,-2.1291,1.7834,-0.5435};//分母系数
  double w;//角频率
  int i,j;
  double complex fenzi,fenmu,H_z;
  for(i=0; i<2001; i++)
    {
      w=2*pi*i/2000.0;
      fenzi=0;
      fenmu=0;
      for(j=0; j<4; j++) 
        { //计算传递函数的分子分母
          fenzi=fenzi+b[j]*cexp(-I*j*w);
		  fenmu=fenmu+a[j]*cexp(-I*j*w);
        }
	  H_z=fenzi/fenmu;//计算分子除以分母得复传递函数值
      printf("%e\t%e\t%e\n",w/pi,20*log10(cabs(H_z)),carg(H_z));
    }
}

用gcc编译。gnuplot作图，一切同前。结果相同。

C99 中已经引入了复数类型，使用时需要包含<complex.h>

定义复数变量时，用：

double complex v1=3.1+5*I;

除了基本的加减乘除可以直接运算外，还有以下函数可以使用：

• 复三角函数

  • 反余弦
    • cacos 双精度版本
    • cacosf 单精度版本
    • cacosl 长双精度版本
  • 反正弦
    • casin 双精度版本
    • casinf 单精度版本
    • casinl 长双精度版本
  • 反正切
    • catan 双精度版本
    • catanf 单精度版本
    • catanl 长双精度版本
  • 余弦
    • ccos 双精度版本
    • ccosf 单精度版本
    • ccosl 长双精度版本
  • 正弦
    • csin 双精度版本
    • csinf 单精度版本
    • csinl 长双精度版本
  • 正切
    • ctan 双精度版本
    • ctanf 单精度版本
    • ctanl 长双精度版本

• 复双曲函数

  • 反双曲余弦
    • cacosh 双精度版本
    • cacoshf 单精度版本
    • cacoshl 长双精度版本
  • 反双曲正弦
    • casinh 双精度版本
    • casinhf 单精度版本
    • casinhl 长双精度版本
  • 反双曲正切
    • catanh 双精度版本
    • catanhf 单精度版本
    • catanhl 长双精度版本
  • 双曲余弦
    • ccosh 双精度版本
    • ccoshf 单精度版本
    • ccoshl 长双精度版本
  • 双曲正弦
    • csinh 双精度版本
    • csinhf 单精度版本
    • csinhl 长双精度版本
  • 双曲正切
    • ctanh 双精度版本
    • ctanhf 单精度版本
    • ctanhl 长双精度版本

• 指数

  与

  对数

  函数

  • 指数
    • cexp 双精度版本
    • cexpf 单精度版本
    • cexpl 长双精度版本
  • 自然对数
    • clog 双精度版本
    • clogf 单精度版本
    • clogl 长双精度版本

• 幂

  运算和

  绝对值

  • 绝对值
    • cabs 双精度版本
    • cabsf 单精度版本
    • cabsl 长双精度版本
  • 幂运算
    • cpow 双精度版本
    • cpowf 单精度版本
    • cpowl 长双精度版本
  • 平方根
    • csqrt 双精度版本
    • csqrtf 单精度版本
    • csqrtl 长双精度版本

• 操作

  • 相角

    • carg 双精度版本
    • cargf 单精度版本
    • cargl 长双精度版本

  • 虚部

    • cimag 双精度版本
    • cimagf 单精度版本
    • cimagl 长双精度版本

  • 复共轭

    • cong 双精度版本
    • congf 单精度版本
    • congl 长双精度版本

  • 黎曼球面

    投影

    • cproj 双精度版本
    • cprojf 单精度版本
    • cprojl 长双精度版本

  • 实部

    • creal 双精度版本
    • crealf 单精度版本
    • creall 长双精度版本