参考文章
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