C语言结合gnuplot设计巴特沃斯型模拟低通滤波器

设计步骤

在已知通带截频Wp、阻带截频Ws、通带衰减Ap和阻带衰减As等设计指标的前提下，设计BW型模拟低通滤波器的步骤如下：

1. 确定滤波器阶数
   

2. 确定3dB截频Wc
   

3. 求系统函数的极点
   

4. 求系统函数
   

设计一个满足下列指标的巴特沃斯型模拟低通滤波器：
fp=1kHz，fs=2kHz，Ap<=1dB，As>=40dB

matlab代码

Wp = 2*pi*1000; Ws = 2*pi*2000; Ap = 1; As = 40;
[N, Wc] = buttord(Wp, Ws, Ap, As, 's');	%确定滤波器阶数N和3dB截频Wc，'s'表示设计模拟滤波器
fprintf('Order of the filter = %.0f\n',N);	%输出滤波器阶数N
[num, den] = butter(N, Wc, 's');	%求传输函数分子和分母多项式的系数，'s'表示设计模拟滤波器
omega = [Wp Ws]; 
h = freqs(num, den, omega);	%求通带和阻带边界频率这两点的频率衰耗
fprintf('Ap = %.4f\n', -20*log10(abs(h(1))));	%将幅度值转化为dB值
fprintf('As = %.4f\n', -20*log10(abs(h(2))));	%将幅度值转化为dB值
omega = [0 : 200 : 3000*2*pi];
h = freqs(num, den, omega);
gain = 20*log10(abs(h));
plot(omega/(2*pi), gain);
title('巴特沃斯型低通滤波器');
xlabel('Frequency in Hz');
ylabel('Gain in dB');

运行结果：

Ap=0.6167，As=40.0000，满足设计指标。

matlab绘图

C程序

#include<stdio.h>
#include<math.h>
#include<complex.h>
#define pi 3.1416
int main()
{
    double Wp = 2*pi*1000, Ws = 2*pi*2000, Ap = 1, As = 40;
    //确定滤波器阶数N
    int N = ceil(log10((pow(10,0.1*As)-1)/(pow(10,0.1*Ap)-1))/(2*log10(Ws/Wp)));
    //确定3dB截频Wc
    double Wc = Ws / (pow(pow(10,0.1*As)-1,1.0/(2*N)));
    int i;
    double complex s[100];
    //求系统函数的极点
    for(i=1; i<=N; i++)
        s[i] = Wc*(cexp(I*pi*(1.0/2+1.0*(2*i-1)/(2*N))));
    //求系统函数
    int w;
    for(w=0; w<=2*pi*3000; w++)
    {
	double complex H=1, t;
	for(i=1; i<=N; i++)
        {
	    t = (-s[i]) / (I*w-s[i]);
            H = H * t;
        }
        printf("%f\t%f\n", w/(2*pi), 20*log10(cabs(H)));
    }
}

gnuplot绘图