数字信号处理第四章
数字信号处理第四章 模拟滤波器频率变换、冲激响应不变法、双线性变换法
4.4 冲激响应不变法
一、基本原理
＝
x(t)
y(t)
取样
取样
x(n) = x(nT)
?
y(n) = y(nT)
?
＝
响应不变
4.4 冲激响应不变法
一、基本原理
其中
取样
其中
另，根据数字系统响应
冲激响应不变原则！
4.4 冲激响应不变法
一、基本原理
模拟滤波器：
（M<N）
部分分式分解
冲激响应不变准则：
数字滤波器：
因此，双线性变换不改变系统稳定性
4.4 双线性变换法
4、频率预畸变
0
高频进行压缩
无混叠，有畸变
频率越高，畸变越大
预畸变
预畸变公式：
根据数字滤波器设计指标，求对应模拟滤波器设计指标时，需预先进行畸变
4.4 双线性变换法
5、双线性变换法设计滤波器步骤
（1）确定数字滤波器技术指标
（Hz表示）
（弧度表示）或
1）带通：计算几何中心
0
若
，则
代替
若
，则
代替
若
，则令
4.2.4 模拟滤波器的频率变换
带通带阻滤波器衰减参数选择
几何对称：
若实际给出的指标不满足几何对称，如何应对？
2）带阻：计算几何中心
0
若
，则
代替
若
，则
代替
若
，则令
固定靠近
的两个值
以让过渡带更窄为选择标准（靠近中心，指标更严）
模拟转数字滤波器
已知一个模拟滤波器H(s)，如何得到数字滤波器H(z)？
3）设计归一化低通滤波器，得到传输函数
数字信号处理DSP第4章
k 0,1, , N 1
2
13
4.2 按时间抽取（DIT）的基2–FFT算法
将系数统一为 WNk 2 WN2k ,则可得
x[0]
N 4点
x[4]
DFT
G[0]
X [0]
G[1]
X [1]
x[2]
N 4点
WN0
x[6]
DFT
WN2
G[2]
1 G[3]
1
X [2] X [3]
x[1]
N 4点
X m1[i] WNr X m1[ j] , X m1[i] WNr X m1[ j]
m 1, 2 ,
每一个蝶形需要一次复数乘法和两次复数加法。
17
4.2 按时间抽取（DIT）的基2–FFT算法
N点的DIT-FFT计算量为
复数乘法：
1
N 2
log2
N
N 2
复数加法：
2
N 2
log2
N
N
例： 如果每次复数乘法需要100us，每次复数加法需要20us，来 计算N=1024点DFT，则需要
12
4.2 按时间抽取（DIT）的基2–FFT算法
同理
( N 4)1
( N 4)1
G[k] DFT[g[r]]
g[2l]WN2lk2
g[2l 1]WN(22l1)k
l 0
l 0
( N 4)1
( N 4)1
g[2l]WNlk 4 WNk 2
g[2l 1]WNlk 4 ,
l 0
l 0
k 0,1,
(3) WN0 WN4 WN8 WN12 WN16 WN20 WN24 WN28
或 WN4i i 0,1, 2, 3, 4, 5, 6, 7 (dm 1)
数字信号处理 第4章 FFT基本思想和2种基本的FFT
= −W
W的对称性
W的可约性
2 rk WN rk = WN / 2
长序列变成短序列 若N → 2个N / 2
2 则N 2次复述乘法 →（N / 2）= N 2 / 2次复数乘法 2
从信号的特殊性上考虑
– 如奇、偶、虚、实性
W 0 X (0) X (1) W 0 = X (2) W 0 0 X (3) W
对 N = 2M ， 共可分 M 次，即 m = 0,1,L , M − 1，
8点FFT时间抽取算法信号流图
每一级有 N/2 个如下的“蝶形”单元：
xm ( p )
xm +1 ( p )
W
r N
xm (q)
−1
xm +1 (q )
算法讨论（ “级”的概念、碟形单元、 “组” 的概念、旋转因子的分布、码位倒置）
r =2l ,r =2l +1
A(k ), B(k )
C(k) = D(k) =
N / 4−1 l =0
∑x(4l)W
l =0
lk N/4
, k = 0,1,..., N / 4 −1
N / 4−1
lk x(4l + 2)WN / 4 , k = 0,1,..., N / 4 −1 ∑
k A(k) = C(k) +WN / 2 D(k), k = 0,1,..., N / 4 −1 k A(k + N / 4) = C(k) −WN / 2 D(k), k = 0,1,..., N / 4 −1
x(6)
n N
N n = 0,1,L , 2
由此得到基本 运算单元
g (0) g (1) g (2) g (3)
数字信号处理课件第四章资料
5、时间抽取蝶形运算流图符号
X1(k)
X1(k) WNk X 2 (k)
X 2 (k )
WNk
1 X1(k) WNk X 2 (k)
返回DIF 返回例题
设 N 23 8
X1(k)
X 2 (k )
WNk
k 0
W80
1
W81
2
W82
3
W83
X (k)
k 0,1,,7
l0
l 0
X1(k) X 3(k) WNk X 4 (k)
2
X1(k
N 4
)
X 3 (k )
W Nk
2
X
4
(k)
k 0,1,..., N 1 4
x2(r)也进行同样的分解：
x5 (l) x2 (2l)
x6 (l) x2 (2l 1)
l 0,1,..., N 1 4
)
N
/ 21
x1(r)WNrk/ 2
X1(k)
r 0
r 0
X2(k N / 2) X2(k) X (k) X1(k) WNk X 2 (k)
W (kN N
/
2)
WNkWNN
/
2
WNk
N点X(k)可以表示成前 N点和后 点N 两部分：
2
2
前半部分X(k)：
X (k) X1(k) WNk X 2 (k)
N 1
X (k) x(n)WNnk k = 0, 1, …, N-1
n0
x(n)
1 N
N 1
X (k )WNnk
k 0
n = 0, 1, …, N-1
二者的差别只在于WN 的指数符号不同，以及差一 个常数因子1/N，所以IDFT与DFT具有相同的运算量。
数字信号处理第四章
第四章线性时不变离散时间系统的频域分析一、传输函数和频率响应clear;M = input('Enter the filter length M: ');w = 0:2*pi/1023:2*pi;num = (1/M)*ones(1,M);den = [1];h = freqz(num, den, w);subplot(2,1,1)plot(w/pi,abs(h));gridtitle('Magnitude Spectrum |H(e^{j\omega})|')xlabel('\omega /\pi');ylabel('Amplitude');subplot(2,1,2)plot(w/pi,angle(h));gridtitle('Phase Spectrum arg[H(e^{j\omega})]')xlabel('\omega /\pi');ylabel('Phase in radians');M=2 M=10 M=15幅度谱为偶对称，相位谱为奇对称，这是一个低通滤波器。
M越大，通带越窄且过渡带越陡峭。
使用修改后的程序P3.1，计算并画出当w=[0,pi]时传输函数的因果线性时不变离散时间系统的频率响应。
它表示哪种类型的滤波器？w = 0:pi/511:pi;num = [0.15 0 -0.15];den = [1 -0.5 0.7];如下列图1这是一个带通滤波器。
图1 图2Q4.3对下面的传输函数重做习题Q4.2：，式〔4.36〕和式〔4.37〕给出的两个滤波器之间的区别是什么？你将选择哪一个滤波器来滤波，为什么？w = 0:pi/511:pi;num = [0.15 0 -0.15];den = [0.7 -0.5 1];如上图2也是一个带通滤波器，这两个滤波器的幅度谱是一样的，相位谱不太一样，我会选择第一个带通滤波器，因为它的相位谱更加平滑，相位失真小。
《数字信号处理》第四章 相关分析
r12 ( )
x1
(t
)
x2
(t
)dt
x1 (t
)x2 (t)dt
r21( )
x1
(t
)x2 (t)dt
x1
(t
)
x2
(t
)dt
以及
rxx ( )
x(t)x (t )dt x (t)x(t )dt
一、自相关函数的性质
1、自相关函数rxx(τ)的极大值在τ=0处，是实数。
rxx ( ) rxx (0)
证明：
rxx ( )
x(t)x (t )dt
x2 (t)dt x2 (t )dt
x2 (t)dt
x2 (u)du
y2 (t)dt
将其代入均方误差，得到这种近似的最小均方误差为
min
xe2 (t)
x2 (t)dt
x(t
)
y(t
)dt
2
y2 (t)dt
第一节 相关
式中右边第一项 x2 (t)dt表示了原信号x(t)的能量。
若将上式用原信号能量归一化成为相对误差，则有
第二节 相关函数的性质
这是由于：
① r(τ)完全由它的能量谱或功率谱P(f )来决定； ② P(f ) =∣X(f )∣2
具有相同的振幅谱而不同相位谱的信号，可以 有相同的自相关函数。
第一节 相关
此时，相关函数r(τ)具有如下性质：
r12 ( ) r21( )
数字信号处理_程佩青_PPT第四章
主要内容
DIT-FFT算法 DIF-FFT算法 IFFT算法 Chirp-z算法 线性卷积的FFT算法
§4.0 引言
FFT: Fast Fourier Transform
1965年，Cooley&Turky 发表文章《机器计算傅 里叶级数的一种算法》，提出FFT算法，解决 DFT运算量太大，在实际使用中受限制的问题。 FFT的应用。频谱分析、滤波器实现、实时信 号处理等。 DSP芯片实现。TI公司的TMS 320c30，10MHz 时钟，基2-FFT1024点FFT时间15ms。
又WN
k
N 2
W
N /2 N
W W
k N
k N
k X (k ) X1 (k ) WN X 2 (k )，k 0,1,2,...N / 2 1 (2) X ( N k ) X ( N k ) W ( N / 2 k ) X ( N k ) 1 N 2 2 2 2 k X1 (k ) WN X 2 (k )，k 0,1,2,...N / 2 1
n为偶
n为奇
N / 2 1
rk k rk x ( r ) W W x ( r ) W 1 N /2 N 2 N /2 r 0 r 0 X1 ( k )
N / 2 1
2 rk rk （这一步利用： WN WN /2
） r , k 0,1,...N / 2 1
N为2的整数幂的FFT算法称基-2FFT算法。
将序列x(n)按n的奇偶分成两组：
x1 (r ) x(2r ) ，r 0， 1， 2， ...N/ 2 1 x2 (r ) x(2r 1)
数字信号处理第四章-数字滤波器的结构
3).H (z)
Y (z) X (z)
(1 bz1) (1 az1)
y(n) ay(n 1) x(n) bx(n 1)
9
10
11
w w
12
转置流图：
w(n) y(n)
原流图：
w(n) ay(n 1) x(n) bx(n 1) 两边作Z变换：
w(n) x(n) aw(n 1) y(n) w(n) bw(n 1) 两边作Z变换：
乘法系数为复数，运算量增加； 系统的稳定性依赖于零、极点相互抵消，对实
现的精度要求很高。在存在有限字长效应的情 况下，有可能造成系统不稳定。
54
确保所有零点、极点在单位圆内。 55
（h(n)为实数）
第k对 极点， 即第k 个与第 N-k个 谐振器 合并
56
谐振频 率不变
还有两点需要注意：（存在实根） 57
1
前言
线性时不变系统用单位冲击响应来表示 系统函数实际上单位冲击响应的Z变换 系统函数反映线性时不变系统的特性 大多数的信号处理可看成是对信号的滤波操作 数字滤波器实际上就是线性时不变系统
因此数字滤波器可以表示为：
2
前言
M
bk zk
H(z) Y(z) / X (z)
k 0 N
1 ak zk
从信号流图中：
可以清楚地看到系统中的运算步骤和运 算结构。FFT时用到了该特点。
运算结构可以直观反映所需的存储单元 和运算次数。由于是数字实现，必然存 在系统误差，运算结构同时也可以反映 系统误差的累积问题。 下面讨论的IIR和FIR滤波器结构将涉及 上述问题。
14
1
15
无限冲击响应滤波器的特点
82
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第四章线性时不变离散时间系统的频域分析一、传输函数和频率响应例4.1传输函数分析Ｑ4.1ｃｌear;M = input('Enter the filteｒ lｅｎgtｈＭ: ')；w ＝ 0:2*pi/１０23:2*ｐi；nuｍ = (１/M)*ｏｎｅs(1，Ｍ)；dｅn ＝ [1];h = freｑz(nuｍ，den, w）;ｓｕｂplot（２,１,１)plot(w／pi,aｂs（ｈ));gridtiｔｌｅ（'Mａgｎiｔude Sｐectrum ｜Ｈ（e^{j\ｏmeｇa})｜')xlａbel('\omeｇa ／\pｉ＇);yｌabel('Amｐlitude');sｕbplｏt(２,１,２)plｏt(w/pi,aｎgｌe(h))；ｇriｄtitle('Pｈaｓe Speｃｔrum arg［H(e^{j\omega})]')xｌａbｅl('\omeｇａ／\pi');ylabel('Phase in rａdianｓ＇）;Ｍ=２M=1０M=15幅度谱为偶对称，相位谱为奇对称,这是一个低通滤波器。
M越大,通带越窄且过渡带越陡峭。
Q4．2使用修改后的程序Ｐ3.1,计算并画出当w=［0,pｉ]时传输函数的因果线性时不变离散时间系统的频率响应。
它表示哪种类型的滤波器？w = 0:pi／５1１:pi;nuｍ =[0.１５ 0 －0.１5];deｎ = ［１ -0.５０．７]；如下图１这是一个带通滤波器。
图1图2Q4.3对下面的传输函数重做习题Q4.２：，式(4.36）和式(４.37)给出的两个滤波器之间的区别是什么？你将选择哪一个滤波器来滤波,为什么？w = 0：pi／5１1:pi；nuｍ = ［０.15 0 -０.15];den = [０.7 －0.５ 1]；如上图2也是一个带通滤波器,这两个滤波器的幅度谱是一样的，相位谱不太一样,我会选择第一个带通滤波器,因为它的相位谱更加平滑，相位失真小。
Q４.4 使用MAＴLＡB计算并画出当ｗ=［0,pi]时因果线性时不变离散时间系统的群延迟。
系统的传输函数为。
ｃlf;w ＝０：ｐi/511:ｐｉ;num = ［1 -1.２1］;den = [1 －1.3 1.04 -0.２2２]；h= ｇrpｄelay(nｕm，den,w);plｏt(w/ｐｉ,ｈ);xｌabｅl('ｗ／ｐｉ＇）;ylａｂel('群延迟'）；Q4．５使用Ｑ３．5０中编写的程序，分别计算并画出式(4.36)和式（4．37)确定的两个滤波器的冲激响应中的前一百个样本。
讨论你的结果。
clｆ;nuｍ = [0.15 0 -0.1５];den ＝ [０.7 -０.5 1］;L = iｎｐut（'输入样本数 L: ＇);［ｇ t］＝ impz(ｎum,den,L）;stem（t,ｇ);tｉtle(［'前 ',num2sｔｒ(L),' 脉冲响应的样本＇]);xlaｂｅl('时间序号ｎ＇）;ylabel('h[n］'）;(4.36）式（4.37）式由图可知：这些情节由impz给生成的因果的脉冲响应实现的H(ｚ)。
我们观察到Q4.3因果滤波器与H(z）在(４.36）稳定，这意味着Ｈ[n]是绝对可和,我们看到交替和指数衰减的脉冲响应。
在另一方面,因果编档人员与H（ｚ)在(4.37)极点以外的单位圆,是不稳定的。
不足为奇的是,相应的ｈ[n］上图显示与n指数增长。
Q4．6 传输函数的极零点图同样能分析线性时不变离散时间系统的性质。
使用命令zplａｎe可以很容易地得到系统的极零点图。
使用zpｌane分别生成式（4.36)和式(4.3７）确定的两个滤波器的极零点图。
讨论你的结果。
clf;ｎｕm = ［0．1５ 0 -0.１5］；deｎ＝［1 -0.５0．7];[z ｐ k］＝ tｆ2zpｋ(num，den）;ｄｉsp(＇Zeros：＇)；diｓp（z);ｄiｓp('Ｐｏｌes：＇);dｉsp(p);inpｕt('Hit <rｅturn＞ｔｏ conｔiｎue．.．');［sｏs ｋ] = ｚｐ2soｓ(ｚ,ｐ,k）input('Hit <reｔｕrn> tｏｃontiｎue.．．'）;zｐｌａｎe（ｚ,p);式(４．36）式(4.37)由图可知：过滤器在(4.３６)在单位圆和两极因此它的因果实现稳定；较低的图显示过滤器(4.3７)极点在单位圆外,其因果关系的实现是不稳定的。
二、传输函数的类型例4.２滤波器Ｑ4．7ｃlｆ;fc = ０.25;ｎ= ［-6.5:1:6.５］;y = 2＊fc*siｎｃ（2*fｃ＊n)；k = n+６.5;stｅm(k,y）；ｔitle('Ｎ＝ 1４');axis(［０13 -０.2 0.６]);ｘlabel(＇Time ｉndｅx n');ylａbel（＇Aｍｐｌituｄe');gｒid;图1 图2如图１低通有限冲激滤波器的长度为14,决定滤波器长度的语句为ｎ = ［-6．5:1:６．5]，而控制截止频率的参数是ｆｃ = 0.25。
Ｑ4.８fc = 0．45;n ＝[-９.５：1:9.5］;y ＝ 2*fc*ｓｉｎc(2*fc*n）;k ＝ｎ+9.5；ｓteｍ(ｋ，y）；tｉtｌe('N = ２０');axiｓ（[0 １9 -０．2 0.7］）;xｌabeｌ(＇Ｔime iｎdｅx n＇)；ylａbel(＇Ampliｔｕde＇);grid;修改参数fc和ｎ,得到如上图2,可知低通有限冲激滤波器的长度变为2０．Ｑ４.９clf;fc = 0.６５;ｎ = [-7.０:１:7.0]；ｙ = 2*ｆc*sinc(2*fｃ*ｎ)；k = n＋７．０;ｓtｅm（ｋ,y)；titlｅ(＇N = 1４＇);ａxis（[0 14 －0.4 1.4]）；xｌabｅl（'Tｉme ｉndex ｎ');yｌaｂel（'Amplｉtｕde＇）;griｄ;Q4.１0cｌear;Ｎ= inpｕt('Ｅntｅｒｔhｅ filtｅｒtimｅｓｈift N: ');No２ = N/2；fc = ０．25;n = [-No2:1:No2］；y = 2*ｆc*sinｃ(2*fc*n);w = 0:ｐi／5１1:pi；h = freqｚ（y, ［1], ｗ); ｐlot(w／ｐi,abs(ｈ));grid；tiｔle（strcat('｜H（ｅ^{j\omega｝)｜, N＝',nｕm2str(N))）;xｌａｂｅｌ('\ｏmega /＼pi'）；yｌabｅl(＇Ａmpｌiｔuｄe');上图依次分别为N=5,1０,30,１00的四幅图，从这四幅图可以看出随着阶数Ｎ的增大,低通滤波器的过渡带越来越窄，阻带衰减越来越快，滤波器越来越接近理想低通滤波器。
Q4．11ｃｌf;M = 2；ｎｕｍ = ｏｎeｓ（１,M)/Ｍ;［ｇ,ｗ] ＝ gaｉn(nuｍ,1);plｏt（w/pｉ,g）;ｇriｄaxｉs([0 1 -５0 ０.5])xlａbel('\oｍega ／\pｉ');yｌabel（'Ｇain in ｄB');ｔitlｅ(['Ｍ = ',ｎuｍ２ｓtr(M)])可以验证3ｄB截止频率在π/２处。
Q4．１2clear;K = iｎpｕt（'Ｅnter tｈe ｎuｍber of ｓeｃtionｓ K: ＇）；Ｈｚ = [1］;for i＝1：K;Hz ＝ｃoｎv(Hz,[1 1]);end;Hｚ＝（0.5）^K * Hz；[g，w］= gain(Hz,1);ＴhreedＢ＝ -3*oｎes(１，lｅngth(g)）；t1 = 2＊acoｓ（(０.5)^（1／（２*K）))＊ｏｎes(1，5１２)/ｐi;t2 = -5０:50.5/51１:0．５;ｐlot（w／pi,g,w/pi,TｈrｅedB,t１,ｔ2）;griｄ；ａxiｓ([０１ -50 0．5]）xlabel（'＼ｏmega /\pi'）;yｌabel(＇Gａin in dB');ｔiｔle(［'K = ',nｕm２sｔｒ(K),'; Ｔheoreｔical \ｏｍega_{ｃ} = ＇,num2ｓtr （t１(1）)]）;Ｑ4.13ｃlear；M = iｎｐuｔ('Enteｒ the filteｒｌength M: ')；n = 0:M-１;ｎum = （-1）.＾n .* oｎes（1,M）/Ｍ;[ｇ,w] = ｇain(nｕm,1);plot（w/pｉ,g)；grid;axiｓ(［0 1 －50 0.5]);ｘlａbｅl('\ｏmega /＼pｉ＇);ｙlabel('Ｇaｉｎ in dB＇）;tiｔle（['M = '， num２sｔｒ（M)］);其3ｄB截止频率约为0.82piQ4.14 设计一个在０.４５ｐｉ处具有3ｄＢ截止频率ｗc的一阶无限冲激响应低通滤波器和一阶无限冲激响应高通滤波器。
用MＡTＬAB计算并画出它们的增益响应，验证设计的滤波器是否满足指标。
用ＭATLAＢ证明两个滤波器是全通互补和功率互补的。
Q４．１５级联10个式(4.15)所示一阶无限冲激响应低通滤波器,设计一个在０.3pi处具有3dB截止频率ｗc的无限冲激响应低通滤波器。
把它与一个具有相同截止频率的一阶无限冲激响应低通滤波器的增益响应作比较。
Q4.16 设计一个中心频率ｗｏ在0.61ｐi处、3dB带宽为0.5１pi的二阶无限冲激响应带通滤波器。
由于式(４.20）是α的二次方程，为了产生相同的3dＢ带宽,参数α将有两个数值，得到的传输函数HBP(z）也会有两个不同的表达式。
使用函数zplane可产生两个传输函数的极零点图,从中可以选择一个稳定的传输函数。
用MATLＡB计算并画出你所设计的滤波器的增益响应,并验证它确实满足给定的条件。
用设计的稳定无限冲激响应带通滤波器的传输函数的参数α和β,生成一个二阶无限冲激响应带阻滤波器的传输函数HBＳ(ｚ）。
用MＡTLAB证明HＢＰ（ｚ)和ＨBＳ(z)都是全通互补和功率互补的。
Q４.17 用MＡTLAB计算并画出一个梳状滤波器的幅度响应，该梳状滤波器是在Ｌ取不同值的情况下，由式（4.40)给出的原型有限冲激响应低通滤波器得到的。