設傳遞函式$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 長雙精度版本