用双线性变换法设计原型低通为椭圆型的数字IIR带通滤波器

首先,把z=ejω，可得 21?e?j?2???s??jtan???j?T1?e?j?T?2?即S平面的虚轴映射到Z平面的单位圆。

（2.9）

其次，将s=σ+jΩ代入式（2.8），得

2 j? ???z?T2 ???j?T因此

2 (2.10)

?2?2???????T?|z|? 2 ?2?2???????T? (2.11)

由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。

2、设计思路

具体实现方法如下：

（1）确定数字低通滤波器的技术指标：通带截止频率ωp、通带衰减Rp、阻带截止频率ωs、阻带衰减Rs。

（2）将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。如果采用脉冲响不变法，边界频率的转换关系为:

???T (2.12)

如果采用双线性变换法，边界频率的转换关系为 21??tg(?) T 2

(2.13)

（3）按照模拟低通滤波器的技术指标设计模拟低通滤波器。

（4）利用双线性变换法将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。 （5）数字低通技术指标为：

ωp=0.4πrad， Rp=1dB； ωs=0.5πrad， Rs=40dB （6）模拟低通的技术指标为：

通带截止频率ωp=(2/T)*tan(ωp/2)=tan(0.392)=0.0068 阻带截止频率ωs=(2/T)*tan(ωs/2)=tan(0.628)=0.011 通带最小衰减Rp=1dB和阻带最大衰减Rs=40dB； 归一化截止角频率ωp=2pi*Fs/Ft； ωs=2pi*Fs/Ft （7）利用模拟椭圆滤波器设计数字滤波器。通带截止频率为：

ωp=0.4*pi；阻带截止频率为：ωs=0.5*pi；通带最大衰减为：Rp=1；阻带最大衰减为：Rs=40；设定周期为1s；模拟低通滤波器的生成：[num1,den1] =ellip(N,Rp,Rs,Wn,'s')；最后实现输入输出、幅频特性、相频特性的图形。

三 程序流程图

数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标 开始 确定数字低通滤波器的各项性能指标

利用双线性变换法将模拟滤波器转变为数字低通滤波器 按照模拟低通滤波器的技术指标设计模拟低通滤波器

利用模拟椭圆滤波器设计数字滤波器 结束 图2 滤波器设计流程方框图

四 程序源代码

MATLAB 程序如下：

clear all close all

Coeff_LF = 1; % 低频信号系数 F1(300 Hz) Coeff_HF = 1; % 高频信号系数 F2(600 Hz)

Coeff_NOISE = 0; % 噪声系数, 如0(没有噪声), 1, sqrt(2), 3等等.

F1 = 300; % Hz F2 = 600; % Hz

FS = 2000; % Hz 抽样频率

% 时间定义

EndTime = 1/FS * 1023; Time = 0:1/FS:EndTime; N = length(Time);

% 产生正弦信号 x1(t),x2(t)和 x(t)

SinSignal_x1= Coeff_LF*sin(2*pi*F1*Time)% x1(t) SinSignal_x2= Coeff_HF*sin(2*pi*F2*Time)% x2(t)

SinSignalSum = SinSignal_x1 + SinSignal_x2+ Coeff_NOISE*randn(1, N); SinSignal = SinSignalSum;% x(t)

% 绘制信号 x1(t) figure subplot(3,1,1)

plot(Time, SinSignal_x1)

xlabel('时间') xlim([0 0.5]) title('时域信号 x1(t)')

subplot(3,1,2)

plot(Time, SinSignal_x1) xlabel('时间') xlim([0 0.05]) title('时域信号 x1(t)')

FFTSinSignal_x1 = 2/N*fft(SinSignal_x1) subplot(3,1,3)

FnSin_x1 = abs(FFTSinSignal_x1);

FW = 2*pi*(Time)/EndTime;

plot(FW(2:length(FW))*FS/pi/2, FnSin_x1(2:length(FW))) hold on

plot((FW(2:length(FW))-2*pi)*FS/pi/2, FnSin_x1(2:length(FW))) xlim([-FS/2 FS/2])

xlabel('频率')

% 绘制 x2(t) figure subplot(3,1,1)

plot(Time, SinSignal_x2) xlabel('时间') xlim([0 0.5])

title('时域信号 x(2)')

subplot(3,1,2)

plot(Time, SinSignal_x2) xlabel('时间') xlim([0 0.05]) title('时域信号 x2(t)')

FFTSinSignal_x2 = 2/N*fft(SinSignal_x2) subplot(3,1,3)

FnSin_x2 = abs(FFTSinSignal_x2);

plot(FW(2:length(FW))*FS/pi/2, FnSin_x2(2:length(FW))) hold on

plot((FW(2:length(FW))-2*pi)*FS/pi/2, FnSin_x2(2:length(FW))) xlim([-FS/2 FS/2])

xlabel('频率')

% 绘制信号 x(t) figure subplot(3,1,1) plot(Time, SinSignal) xlabel('时间') xlim([0 0.5]) title('时域信号')

subplot(3,1,2) plot(Time, SinSignal)