参考文章
- 绘制数字滤波器的频域响应
- gnuplot
Octave
pkg load signalb=[+1.009874 -1.973835 +0.980993];
a=[+1.000000 -1.973835 +0.990867];Fs = 48000
for k=1:Fsw(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); // realfenzi[1] = fenzi[1] + b[j] * sin(-j * w); // imagfenmu[0] = fenmu[0] + a[j] * cos(-j * w); // realfenmu[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