双线性变换法设计椭圆低通数字IIR带阻滤波器
课程设计名称:数字信号处理
数字信号处理专业课程设计任务书
用双线性变换法设计了一个Butterworth型数字IIR带通滤波器原型。
主要内容
用双线性变换法设计了一个Butterworth型数字IIR带通滤波器原型。
要求通带边界频率为400Hz和500Hz，阻带边界频率分别为350Hz和550Hz，通带最大衰减为1dB，阻带最小衰减为40dB，采样频率为2000Hz。利用MATLAB绘制幅频特性，绘制并分析滤波器传递函数的零极点。
信号通过滤波器，其中450赫兹和600赫兹。滤波器的输出是什么？用Matlab验证你的结论，并给出图表。
任务要求
1，掌握用双线性变换法设计巴特沃兹型原型低通数字IIR带通滤波器的原理和设计方法。
2.求设计滤波器的z变换。
3.用MATLAB绘制幅频特性图。
4.验证设计的滤波器。
参考
1，程培清，数字信号处理教程，清华大学出版社，2001。
2.Sanjit K. Mitra，孙红、余翔宇译，《数字信号处理实验指南》(MATLAB版)，电子工业出版社，2005年6月5438+10月。
3.郭等。，MATLAB 7.x数字信号处理，人民邮电出版社，2006。
4.胡光枢，《数字信号处理的理论算法与实现》，清华大学出版社，2003。
评审意见指导老师签名:
注:本表由指导教师填写，教研室主任审核后发放给入选学生，装订在设计(论文)首页。
填写表格并解释
1.“主题性质”栏:
A.工程设计；
B.工程技术研究；
C.软件工程(如CAI项目等。);
D.文献综述；
E.其他人。
2.“主题来源”栏:
A.部、省、市级以上自然科学基金和科研项目；
B.企事业单位委托的项目；
C.校、院(系、部)级基金项目；
D.自制题目。
1需求分析
本课程的设计需要双线性变换来设计一个巴特沃兹型原型的数字IIR带通滤波器。给定的设计参数分别是400Hz和500Hz的通带边界频率，350Hz和550Hz的阻带边界频率。通带最大衰减为1dB，阻带最小衰减为40dB，采样频率为2000Hz。要求双线性变换方法，使S平面和Z平面有单值一一对应，不存在频谱混淆的问题。由于数字频域和模拟频域的频率不是线性的，这种非线性关系使得通带截止频率和过渡带边缘频率的相对位置非线性失真。在设计巴特沃兹带通滤波器后，让信号通过滤波器，包括450Hz和600Hz，分析滤波器的输出是什么，用Matlab验证结论并给出图形。
2概要设计
1 & gt；双线性变换法设计IIR数字滤波器
脉冲响应不变性方法的主要缺点是频率响应的混叠失真。这是由S平面到Z平面的多值映射关系造成的。为了克服这个缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/t ~ π/t，然后用z=esT将其变换到z平面。也就是说，第一步，将整个S平面压缩映射成S1平面的-π/t ~ π/t的水平带；第二步，通过标准变换关系z=es1T将横带变换到整个Z平面。这样就建立了S平面和Z平面一对一的单值关系，消除了多值变换和频谱混叠。映射关系如图1-3所示。
图1-3双线性变换的映射关系
为了将S平面的整个虚轴J ω压缩到S1平面的J ω 1轴上的-π/T到π/T段，可以实现如下的切线变换。
(1-5)
其中t仍然是采样间隔。
当ω1从-π/T经0变为π/T时，ω从-∞经0变为+∞，即映射整个J ω轴。写出公式(1-5)
解析地将这种关系推广到整个S平面和S1平面，使得J ω = S，J ω 1 = S1，我们得到
然后通过下面的标准变换关系将S1平面映射到Z平面。
z=es1T
这样，S平面和Z平面之间的单值映射关系得到如下:
(1-6)
(1-7)
等式(1-6)和(1-7)是S平面和Z平面之间的单值映射关系。这些变换是两个线性函数的比值，所以称为双线性变换。
公式(1-5)和公式(1-6)的双线性变换满足映射变换的两个要求。
首先，设z = EJω，可以得到
(1-8)
即S平面的虚轴映射到Z平面的单位圆上。
其次，将s = σ+j ω代入公式(1-8)，得到
因此
可以看出，当σ
2 & gt双线性变换法的优缺点
与脉冲响应不变量法相比，双线性变换法的主要优点是避免了频率响应的混叠。这是因为S平面和Z平面具有单值一一对应关系。S平面的整个jω轴对应Z平面单位圆的一个周期，即频率轴是单值变换关系。这种关系如公式(1-8)所示，改写如下:
上式表明，S平面上的ω与Z平面上的ω呈非线性相切关系，如图7-7所示。
从图7-7可以看出，模拟角频率ω与数字频率ω的变换关系在零频率附近接近线性；然而，当ω进一步增加时，ω增长越来越慢。最后，当ω→∞时，ω在折叠频率ω = π处结束，所以双线性变换不会因为高频部分超过折叠频率而与低频部分混淆，从而消除了频率混叠现象。
图1-4双线性变换法的频率变换关系
而双线性变换的这种特性是通过频率的严重非线性关系得到的，如方程(1-8)和图1-4所示。因为这种频率之间的非线性变换关系，产生了新的问题。首先对一个具有线性相位的模拟滤波器进行双线性变换，得到一个具有非线性相位的数字滤波器，它不再保持原来的线性相位；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数，即某一频带的幅频响应近似等于某一常数(这是典型的低通、高通、带通和带阻滤波器的响应特性)，否则变换产生的数字滤波器的幅频响应相对于原模拟滤波器会发生畸变，如图1-5所示。
图1-5双线性变换法非线性映射振幅和相位特性
对于分段常数滤波器，经过双线性变换后，仍然得到了具有分段常数幅频特性的滤波器，但各段边缘的临界频率点发生了畸变，可以通过频率预畸变进行校正。即预先对临界模拟频率进行失真处理，然后经过变换映射到所需的数字频率。
根据上述IIR数字滤波器的设计方法，利用双线性变换法，基于MATLAB设计了一个IIR带通滤波器，通带截止频率为FP 1 = 500 Hz，FP 2 = 400 Hz。阻带截止频率fs1 = 350hz，fs2 = 550hz通带最大衰减αp = 1dB；阻带最小衰减αs=40dB，采样频率为2000HZ。
一、设计步骤:
(1)根据任务确定绩效指标:
通带截止频率FP 1 = 500 Hz，FP2 = 400Hz；
阻带截止频率fs1 = 350hz，fs2 = 550hz
通带最大衰减αp = 1dB；
阻带最小衰减αs = 40dB；
采样频率为2000赫兹。
(2)用ω = (2/t) tan (w/2)对带通数字滤波器H(z)的数字边界频率进行预失真，得到带通模拟滤波器H(s)的边界频率主要是通带截止频率fp1，fp2阻带截止频率fs1，fs2的转换。
为了简化计算，双线性变换法一般是T=2s。
通带截止频率fc 1 =(2/t)* tan(2 * pi * FP 1/2)
fc2=(2/T)*tan(2*pi*fp2/2)
阻带截止频率FR 1 =(2/t)* Tan(2 * PI * FS 1/2)
fr2=(2/T)*tan(2*pi*fs2/2)
阻带最小衰减αs=1dB，通带最大衰减αp = 40dB；
(3)利用低通带通频率转换公式λ = ((f 2)-(f 0 2))/f)将模拟带通滤波器指数转换为模拟低通滤波器指数。
B=fc2-fc1
normwr1=(((fc1*fc2)-(fc1^2))/fc1)normwr2=(((fc1*fc2)-(fc2^2))/fc2)
normwc1=(((fc1*fc2)-(fr1^2))/fr1)
normwc2=(((fc1*fc2)-(fr1^2))/fr1)
模拟低通滤波器的指标:normfc，normfr，αp，α s。
(4)设计模拟低通原型滤波器。借助巴特沃兹滤波器，使用模拟低通滤波器。
模拟低通滤波器的传递函数H(s)通过滤波器设计方法获得。
(5)调用lp2bp函数将模拟低通滤波器转换成模拟带通滤波器。
(6)通过双线性变换将模拟带通滤波器H(s)转换成数字带通滤波器H(z)。
二。程序流程图:
开始
↓
读入数字滤波器的技术指标
↓
将该指数转换成归一化模拟低通滤波器的指数。
↓
设计归一化模拟低通滤波器阶数n和1db截止频率。
↓
模拟域频率变换，将G(P)变换为模拟带通滤波器H(s)
↓
用双线性变换法将H(s)转换成数字带通滤波器H(z)。
↓
输入信号后，显示相关结果。
↓
结束
3操作环境
PC MATLAB
4开发工具和编程语言
MATLAB语言
5详细设计
clc全部清除；
%设计过滤器
要设计的带通滤波器参数的%。
Fres = 2000
ts = 1/Fres；
Omgap1=500
Omgap2=400
Omgas1=350
Omgas2=550
%双线性变换，使得S平面和Z平面是单值一一对应的。
Omgap=(Omgap1-Omgap2)*Ts
om gap 3 = tan(om gap 1 * Ts * 2 * pi/2)
Omgap4=tan(Omgap2*Ts*2*pi/2)
OMG as 3 = tan(OMG as 1 * Ts * 2 * pi/2)
Omgas4=tan(Omgas2*Ts*2*pi/2)
参数标准化百分比
ap1=Omgap3/Omgap
ap2=Omgap4/Omgap
as1=Omgas3/Omgap
as2=Omgas4/Omgap
%将模拟带通滤波器指数转换为模拟低通滤波器指数。
I 1 =(AP 1 * AP2-as 1 * as 1)/as 1
I2=(ap1*ap2-as2*as2)/as2
i3 =(AP 1 * AP2-AP 1 * AP 1)/AP 1
I4=(ap1*ap2-ap2*ap2)/ap2
I5 =最小值(I1，-I2)
alf PP = 1；
alfps = 40
%设计巴特沃兹低通滤波器HL(s)
[N，Wn]= but ord(I4，I5，alfpp，alfps，' s ')
黄油
%将设计的模拟低通滤波器传递函数HL(s)转换成模拟带通滤波器H(s)。
[Num2，Den2]=lp2bp(Num，Den，sqrt(Omgap3*Omgap4)，om gap)；
%将设计的模拟带通滤波器H(s)转换成数字带通滤波器H(Z)。
[Num3，den 3]=双线性(Num2，Den2，0.5)；
%画出H(Z)的幅频特性曲线
w = 0:pi/255:pi；
h=freqz(Num3，Den3，w)；
h = 20 * log 10(ABS(h))；
plot(w，H)；
%设计信号功能
f 1 = 450；f2 = 600
t = 0:1/2000:40/F2；
f = sin(2 * pi * f 1 * t)+sin(2 * pi * F2 * t)；
情节(f)；
%的信号通过带通滤波器。
g=invfreqz(高，宽，40，50)
g1=conv(g，f)
绘图(g1)
6调试分析
设计滤波器时，计算幅频特性。然而，相位-频率特性不是通过频谱来分析的。但是信号和经过滤波器的信号的图形分析是时间和时间分析，没有频域分析，所以不知道用哪个函数进行傅里叶变换。我试了freqz，fft等函数，都没出来。呼叫格式不正确。
本实验采用巴特沃斯低通滤波器的设计方法，也可以采用切比雪夫低通滤波器，或者采用椭圆滤波器。
在这个实验中，我们只使用两个频率的简单叠加作为输入信号，但实际信号是多个频率的组合，还可以添加更多的频率。
现实中应该还有噪音的影响，这个实验不考虑。可以添加噪声信号。
7测试结果
频率=
2000
Ts =
5.0000e-004
Omgap1 =
500
Omgap2 =
名流
Omgas1 =
350
Omgas2 =
550
Omgap =
0.0500
Omgap3 =
1.0000
Omgap4 =
0.7265
Omgas3 =
0.6128
Omgas4 =
1.1708
ap1 =
20.0000
ap2 =
14.5309
as1 =
12.2560
as2 =
23.4170
I1 =
11.4562
I2 =
-11.0065
I3 =
-5.4691
I4 =
5.4691
I5 =
11.0065
alfpp =
1
alfps =
40
N =
八
Wn =
6.1894
数量=
1.0e+006 *
列1到8
0 0 0 0 0 0 0 0
第9栏
2.1538
Den =
1.0e+006 *
列1到8
0.0000 0.0000 0.0005 0.0052 0.0377 0.1984 0.7386 1.7837
第9栏
2.1538
Num2 =
1.0e-004 *
列1到8
0.8413 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000
第9栏
-0.0000
Den2 =
列1到8
1.0000 1.5863 7.0705 8.7151 20.5005 19.9985 32.1353 24.8474
列9到16
29.9185 18.0527 16.9631 7.6698 5.7123 1.7643 1.0400 0.1695
列17
0.0776
Num3 =
1.0e-004 *
列1到8
0.0043 0.0000 -0.0341 0.0000 0.1194 0.0000 -0.2389 0.0000
列9到16
0.2986 0.0000 -0.2389 0.0000 0.1194 0.0000 -0.0341 0.0000
列17
0.0043
Den3 =
列1到8
1.0000 -2.2463 8.3946 -13.4519 27.6744 -33.6694 48.3386 -45.7295
列9到16
49.5687 -36.4140 30.6576 -16.9904 11.1207 -4.2944 2.1332 -0.4524
列17
0.1603
h =
列1到5
-0.0000-0.0000+0.0000 I-0.0000 I-0.0000+0.0000 I-0.0000+0.0000 I
……
第251至256列
-0.0000-0.0000 I-0.0000 I-0.0000 I-0.0000 I-0.0000 I-0.0000 I-0.0000 I-0.0000 I-0.0000 I-0.0000 I-0.0000 I
H =
列1到8
-279.2737 -279.2732 -279.2723 -279.2813 -279.3855 -280.0020 -283.0382 -292.4507
……
第249至256列
-281.0032 -280.3352 -280.1410 -280.0979 -280.0943 -280.0973 -280.0996 -280.1004
f1 =
450
f2 =
600
t =
列1到8
0 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035
……
列129到134
0.0640 0.0645 0.0650 0.0655 0.0660 0.0665
f =
列1到8
0 1.9387 -0.2788 -1.4788 0.3633 0.7071 -0.1420 0.1338
……
列129到134
-0.3633 -0.7946 1.0000 1.1075 -1.5388 -1.0418
g =
列1到8
0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0001 -0.0002
……
第33列到41列
-0.0002 0.0000 0.0001 -0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000
g1 =
列1到8
0 0.0000 0.0000 -0.0000 -0.0000 0.0001 0.0001 -0.0003
……
列169到174
0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000
综合以上分析，数据基本正确。
上图为带通滤波器的幅频响应。
、
上图显示了传递函数的零极点分析。
上图是输入信号的时域。
上图显示了输出波形的时域分析。
参考
[1]吴大正《信号与线性系统分析》，第四版高等教育出版社
[2]郑·《信号与系统》第二版高等教育出版社。
[3] Sanjit K. Mitra《数字信号处理——基于计算机的方法》第三版清华大学出版社。
[4]余《数字信号处理与MATLAB实现》清华大学出版社。
[5]周·《数字信号处理基础》北京邮电大学出版社
[6]莱昂的《数字信号处理》英文第二版机械工业出版社。
感受和经历
................0.0