freqz之C实现例程

参考文章

• 绘制数字滤波器的频域响应
• gnuplot

Octave

pkg load signal

b=[+1.009874 -1.973835 +0.980993];
a=[+1.000000 -1.973835 +0.990867];

Fs = 48000
for k=1:Fs
    w(k)=2*pi*(k-1)/(Fs - 1);
    z=exp(j*w(k));
    Z=[1 z^-1 z^-2];
    H_z(k)=sum(b.*Z)./sum(a.*Z);
end

% Transfer from Rad to Hz.
fz = (w / (2 * pi)) * Fs;

% mag2db
Hf_PEAK = 20*log10(abs(H_z));
Hx_PEAK = angle(H_z);

figure(1);plot(fz(1:4000), Hf_PEAK(1:4000));
figure(2);plot(fz(1:4000), Hx_PEAK(1:4000));

Magnitude Vs Frequency

Phase Vs Frequency

C

安装gnuplot

• msys2环境

pacman -S mingw-w64-x86_64-gnuplot

然后，进入gnuplot互交模式。

gnuplot

注意： 实测无法启动，可能是版本问题。

• Octave
  gnuplot官网介绍来看，Octave也基于gnuplot绘图。所以可以直接调用Octave的软件包即可。Command Prompt下直接进入gnuplot互交模式。

E:\Ethan\NOTEs\dev\C\Effendi>d:\Octave\Octave-5.2.0\mingw64\bin\gnuplot.exe

int print_freqz(uint32_t frequency, uint32_t order, double a[], double b[])
{
    double H_z[2]   = {+0.0};
    double fenzi[2] = {+0,0};
    double fenmu[2] = {+0,0};

    double HzAmp  = +0.0;
    double HzAng  = +0.0;
    double tmp    = +0.0;

    double  w;
    double  hz;
    int     i;
    int     j;

    for (i = 0; i < frequency; i++) {
        w = (2.0 * M_PI * i) / frequency;
        fenzi[0] = +0.0;
        fenzi[1] = +0.0;
        fenmu[0] = +0.0;
        fenmu[1] = +0.0;

        for (j = 0; j <= order; 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]);

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

        hz = (w / (2 * M_PI)) * frequency;
        printf("%e\t%e\t%e\n", hz, 20 * log10(HzAmp), HzAng);
    }

    return 0;
}

./Effendi.exe freqz -o 2 -f 48000 -c +1.000000:-1.973835:+0.990867:+1.009874:-1.973835:+0.980993 > freqz.tx

gnuplot互交模式下。

gnuplot> plot [0:4000] [-2:10] "freqz.txt" u 1:2 w l
gnuplot> plot [0:4000] [-0.8:0.8] "freqz.txt" u 1:3 w l

Magnitude Vs Frequency

Phase Vs Frequency