哈工大数字信号处理实验
哈尔滨工程大学数字信号处理实验四
t/T x (n )k X (k)t/T x (n)k X (k )一、 实验原理DFT 的快速算法FFT 利用了WN^(nk)的三个固有特性：(1)对称性，(WN^(nk))*=WN^(-nk),(2)周期性，WN^(nk)=WN^(n+N)k=WN^n(k+K),(3)可约性WN^(nk)=WmN^(nmk)和WN(nk)=WN/m^(nk/m)。
FFT 算法基本上可以分为两大类，即按时间抽选法和按频率抽选法。
MATLAB 中提供了进行快速傅里叶变换的fft 函数，X=fft(x),基2时间抽取FFT 算法，x 是表示离散信号的向量；X 是系数向量； X=fft(x ，N)，补零或截断的N 点DFT ，当x 的长度小于N 时，对x 补零使其长度为N ，当x 的长度大于N 时，对x 截断使其长度为N 。
Ifft 函数计算IDFT ，其调用格式与fft 函数相同，参考help 文件。
例3.1程序及图形文件 >> k=8;>> n1=[0:19];>> xa1=sin(2*pi*n1/k); >> subplot(221) >> stem(n1,xa1)>> xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xa1);>> xk1=abs(xk1);>> subplot(222)>> stem(n1,xk1)>> xlabel('k');ylabel('X(k)'); >> n2=[0:1:15]; >> xa2=sin(2*pi*n2/k); >> subplot(223) >> stem(n2,xa2)>> xlabel('t/T');ylabel('x(n)');>> xk2=fft(xa2);>> xk2=abs(xk2);>> subplot(224)>> stem(n2,xk2)>> xlabel('k');ylabel('X(k)');上两个图为N=20是的截取信号和DFT 结果，由于截取了两个半周期，频谱出现泄漏；下面的两个图为N=16时的截取信号和DFT 结果，由于截取了两个整周期，得到单一谱线的频谱。
哈工程数字信号处理实验二
实验二离散时间傅立叶变换一：实验原理经由正、负离散时间傅立叶变换表达式是信号分析的一个关键部分。
当LTI系统用于滤波的时候，作为冲激响应离散时间傅立叶的频率响应，提供了LTI系统间接的描述。
离散时间傅立叶变换X（）是w的周期复值函数，周期总是2π，并且基周期通常选在区间[-π，π]上。
对离散时间傅立叶变换DTFT来说有两个问题：1.DTFT的定义对无限长信号是有效的。
2.DTFT是连续变量w的函数在MATLAB中，任何信号（向量）必须是有限长度的，所以DTFT 无法用MATLAB直接计算。
当能从变换定义式推导出解析式并计算它时，可以用MATLAB计算。
第二个问题是频率抽样问题。
MATLAB 擅长在有线网络点上计算DTFT，通常选足够多的频率以使绘出的图平滑，逼近真实的DTFT。
二：实验内容1.矩形信号的DTFT设矩形脉冲r[n]由下式定义。
（3.12）A．证明r[n]的DTFT可由式（3.13）得出B．使用DTFT函数计算12点脉冲信号的DTFT。
绘出在区间上对w的DTFT。
把实部和虚部分开绘出，但是注意这些图不是很有用。
另绘出DTFT的幅度。
绘图时，注意要正确地标注频率坐标轴的变量。
C．注意asinc函数零点的位置是规律分布的。
对奇数长脉冲，比如L=15的脉冲重复进行DTFT计算并绘出幅度；同样再次检验零点位置，注意峰值高度。
D．对asinc函数零点的间距与ASINC函数的直流值，确定出通用规则。
4.指数信号对于信号x[n], 使用freqz函数计算其DTFT X(eⁿ)。
A..对w在区间-π<=w<π上绘出幅度与相位特性。
这需要从freqz返回的[X,W] 向量的移位。
解释为什么幅度特性是w的偶函数，而相位特性是w的奇函数。
B.推算一阶系统的幅度特性与相位特性的表达式。
C.直接以这些表达式来计算幅度特性和相位特性，并与用freqz函数计算出的结果相对比。
三．实验程序1.脉冲信号的DTFTA.format compact,subplot(111)L=10;a=ones(1,L);n=0:(L-1); xn=a.^n;[X,W]=dtft(xn,128);subplot(211),plot(W/2/pi,abs(X));grid,title('Test2_1_a_1');xlabel('NORMALIZED FREQUENCY'),ylabel('|H(w)|') subplot(212),plot(W/2/pi,180/pi*angle(X));gridxlabel('NORMALZEDFREQUENCY'),ylabel('DEGEREES')title('Test2_1_a_2')B.format compact,subplot(111)L=12;a=ones(1,L);n=0:(L-1); xn=a.^n;[X,W]=dtft(xn,128);subplot(311),plot(W/2/pi,real(X),'b');grid,title('Test2_1_b_1');xlabel('NORMALIZEDFREQUENCY/w'),ylabel('REAL(X)')subplot(312),plot(W/2/pi,imag(X),'b');grid,title('Test2_1_b_2');xlabel('NORMALIZEDFREQUENCY/w'),ylabel('IMAG(X)')subplot(313),plot(W/2/pi,abs(X));grid,title('Test2_1_a_1');xlabel('NORMALIZED FREQUENCY/w'),ylabel('|H(w)|')C.format compact,subplot(111)L=15;a=ones(1,L);n=0:(L-1); xn=a.^n;[X,W]=dtft(xn,128);subplot(111),plot(W/2/pi,abs(X));grid,title('Test2_1_c_1');xlabel('NORMALIZED FREQUENCY'),ylabel('|H(w)|')D.4.指数信号A.B.C.四．结果分析1.脉冲信号的DTFTA.由dtft函数计算出r[n]的dtft如下图，与标准sinc函数图象一致，故证明得证。
哈工大数字信号处理实验2011
实验一 离散傅里叶变换的性质一、 实验目的1、 掌握离散傅里叶变换的性质，包括线性特性、时移特性、频移特性、对称性和循环卷积等性质；2、 通过编程验证傅里叶变换的性质，加强对傅里叶变换性质的认识。
二、 实验原理和方法 1. 线性特性1212DFT[()()]()()ax n bx n aX k bX k +=+2. 时移特性DFT[()]()DFT[()]()km kmx n m W X k x n m W X k -+=-=3. 频移特性()()nl N IDFT X k l IDFT X k W +=⎡⎤⎡⎤⎣⎦⎣⎦4. 对称性设由x(n) 延拓成的周期序列为 ()x n 则()()()e o xn x n x n =+ 共轭对称序列()()()*12e xn x n x N n ⎡⎤=+-⎣⎦ 共轭反对称序列()()()*12o x n x n x N n ⎡⎤=--⎣⎦ 将()e xn 和()o x n 截取主周期，分别得 ()()()ep e N x n x n R n = ()()()o p o N x n x n R n= 则()()()()()N ep op x n xn R n x n x n ==+ x(n)序列的实部和虚部的离散立叶变换(){}()Re ep DFT x n X k =⎡⎤⎣⎦ (){}()Im op DFT j x n X k =⎡⎤⎣⎦当x(n)为实数序列[][]()(())()()(())()()()(())()(())()()(())()(())()()()arg ()arg ()N N N N R R N N R N N I I N N I N N X k X k R k X k X N k R k X N k X k X k R k X N k R k X k X k R k X N k R k X k X N k X k X k *=-≅-=-≅-=-=-=--=--=-=--5. 循环卷积()312312()()()()()x n x n x n X k X k X k =⊗⇒=有限长序列线性卷积与循环卷积的关系 x1(n)和x2(n)的线性卷积：111212()()()()()N l m m x n x m x n m x m x n m -∞=-∞==-=-∑∑1120()()N m x m x n m -==-∑将x1(n)和x2(n)延拓成以N 为周期的周期序列11()()r xn x n rN ∞=-∞=+∑ 22()()q xn x n qN ∞=-∞=+∑ 则它们的周期卷积为14120()()()N p m x n xm x n m -==-∑ 1120()()N m x m xn m -==-∑ 1120()()N m q x m x n m qN -∞==-∞=-+∑∑1120()()N q m x m x n qN m ∞-=-∞=⎡⎤=+-⎢⎥⎣⎦∑∑ ()lq x n qN ∞=-∞=+∑x1(n)和x2(n)周期延拓后的周期卷积等于他们的线性卷积的的周期延拓。
哈尔滨工程大学数组信号处理实验三
1020304050607080一、实验原理时域抽样定理给出了连续信号抽样过程中信号不失真的约束条件：对于基带信号，信号抽样频率fsam 大于等于2倍的信号最高频率m ，即fsam>=2fm 。
时域抽样是把连续信号x(t)变成适于数字系统处理的离散信号x[k] ；信号重建是将离散信号x[k]转换为连续时间信号x(t)。
1. 信号的时域抽样若x[k]=x(kT)|t=kT ，则信号x(t)与x[k]的频谱之间存在：其中：x(t)的频谱为X(jw)，x[k]的频谱为X(ejW) 可见，信号时域抽样导致信号频谱的周期化。
2. 信号的频域抽样非周期离散序列x[k]的频谱X(ejW)是以2p 为周期的连续函数。
频域抽样是将X(ejW)离散化以便于数值计算。
频域抽样与时域抽样形成对偶关系。
在[0,2p]内对X(ejW) 进行N 点均匀抽样，引起时域序列x[k]以N 点为周期进行周期延拓。
频域抽样定理给出了频域抽样过程中时域不发生混叠的约束条件：若序列x[k]的长度L ，则应有NL 。
二、实验内容（一）抽样引起的混叠正弦信号混叠：可以按抽样频率fs=1/Ts 对x(t)抽样来获得离散时间信号。
a.抽样频率fs=8kHz ，令正弦波频率为300Hz ，然后在10ms 长间隔上抽样，相位任意指定。
使用stem 绘出产生的离散时间信号。
程序如下： >> nn=0:80;>> x=sin(2*pi*300/8000*nn+pi/2);>> stem(nn,x);%时间长度为10ms ，抽样间隔为1/8000s %抽取样本数N=10ms/(1/8000s)=80 %规定范围为0到80%使用stem 绘出离散时间信号b.使用plot 绘图，将点用直线连接起来。
>>plot(nn,x);%使用plot 绘图10203040506070801020304050607080102030405060708010203040506070801020304050607080102030405060708010203040506070801020304050607080c.把正弦的频率从100Hz 变至475Hz ，每次增加125Hz ，绘出一系列相应的图，并用subplot 指令把四个图放在同一张图上。
哈尔滨工业大学威海校区_《数字信号处理》实验一
哈尔滨⼯业⼤学威海校区_《数字信号处理》实验⼀数字信号处理实验报告实验名称：实验⼀离散傅⾥叶变换的性质实验⽇期：2011．11．16姓名：尤伟学号：090240328哈尔滨⼯业⼤学（威海）实验⼀离散傅⾥叶变换的性质⼀、实验⽬的1、掌握离散傅⾥叶变换的性质，包括线性特性、时移特性、频移特性、对称性和循环卷积等性质；2、通过编程验证傅⾥叶变换的性质，加强对傅⾥叶变换性质的认识。
⼆、实验原理和⽅法1.线性特性1212D FT [()()]()()ax n bx n aX k bX k +=+ 2.时移特性DFT[()]()DFT[()]()km kmx n m W X k x n m WX k -+=-=3.频移特性()()nlN IDFT X k l IDFT X k W +=4. 对称性设由x(n)开拓成的周期序列为 ()p x n 则()()()p pe po x n x n x n =+ 偶序列()()()*12pe p p x n x n x N n ??=+-?奇序列()()()*12pop p x n x n x N n ??=--?? 将()pe x n 和()po x n 截取主周期，分别得()()()pet pe N x n x n R n = ()()()p o tp oN x n x n R n =则()()()()()p N pet pot x n x n R n x n x n ==+ x(n)序列的实部和虚部的离散⽴叶变换(){}()R e petD FT x n X k = (){}()Im potj x n Xk =[][]()()()()()()()()()()()arg ()arg ()R R R I I I X k X k X N k X k X k X N k X k X k X N k X k X N k X k X k * =-=-=-=-=--=--=-=-- 5.循环卷积()3123121()()()()()x n x n x n X k X k X k N=?=有限长序列线性卷积与循环卷积的关系 X1(n)和x2(n)的线性卷积：11312120()()()()()N m m x n x m x n m x m x n m -∞=-∞==-=-∑∑112()()N m x m xn m -==-∑将X1(n)和x2(n)开拓成以N 为周期的周期序列11()()p r x n x n rN ∞=-∞=+∑22()()p q x n x n qN ∞=-∞=+∑则它们的周期卷积为14120()()()N p p p m x n xm x n m -==-∑12()()N p m x m xn m -==-∑1120()()N m q x m x n m qN -∞==-∞=-+∑∑1120()()N q m x m x n qN m ∞-=-∞=??=+-∑∑ 3()q x n qN ∞=-∞=+∑X1(n)和x2(n)周期开拓后的周期卷积等于他们的线性卷积的的周期开拓。
哈工程数字信号处理实验六离散时间滤波器
实验六 离散时间滤波器一、实验名称离散时间滤波器设计二、实验目的：1、掌握利用脉冲响应不变法设计IIR 数字滤波器的原理及具体方法。
2、加深理解数字滤波器与连续时间滤波器之间的技术转化。
3、掌握脉冲响应不变法设计IIR 数字滤波器的优缺点及使用范围。
4、掌握利用双线性变换法设计IIR 数字滤波器的原理及具体方法。
5、深入理解利用双线性变换法设计IIR 数字滤波器的优缺点及使用范围。
三、实验原理：1、脉冲响应不变法变换原理脉冲响应不变法将模拟滤波器的s 平面变换成数字滤波器的z 平面，从而将模拟滤波器映射成数字滤波器。
IIR 滤波器的系数函数为1-z （或z ）的有理分式，即 ∑∑=-=--=Nk kk M k k kz a z b z H 101)( 一般满足N M ≤。
⑴转换思路：)(H )()(h )(h )(z n h nT t s H z a a −−→−=−−−→−−−−−→−变换时域采样拉普拉斯逆变换若模拟滤波器的系统函数H （s ）只有单阶极点，且假定分母的阶次大于分子的阶次，表达式：∑=--=Nk T s k z e TA z H k 111)( ⑵s 平面与z 平面之间的映射关系。
Ω+==j s re z j σω→=→=→ΩT j T jw sT e e re e z σTe r T Ω==ωσ IIR 数字滤波器设计的重要环节式模拟低通滤波器的设计，典型的模拟低通滤波器有巴特沃兹和切比雪夫等滤波器。
由模拟低通滤波器经过相应的复频率转换为H （s ），由H （s ）经过脉冲响应不变法就得到所需要的IIR 数字滤波器H （z ）。
Matlab 信号处理工具箱中提供了IIR 滤波器设计的函数，常用的函数：IIR 滤波器阶数选择Buttord--巴特沃兹滤波器阶数选择。
Cheb1ord--切比雪夫I 型滤波器阶数选择。
Cheb2ord--切比雪夫II 型滤波器阶数选择。
IIR 滤波器设计Butter--巴特沃兹滤波器设计。
哈尔滨工程大学数字信号处理实验八
一、实验原理1、频率（周期）检测对周期信号来说，可以用时域波形分析来确定信号的周期，也就是计算相邻的两个信号波峰的时间差或过零点的时间差。
这里采用过零点（ti）的时间差T（周期）。
频率即为f=1/T,由于能够求得多个T的值（ti有多个），故采用它们的平均值作为周期的估计值。
2、幅值检测在一个周期内，求出信号最大值ymax与最小值ymin的差的一半，即A=（ymax-ymin）/2,同样，也会求出多个A值，但第一个A值对应的ymax和ymin不是在一个周期内搜索得到的，故以除以第一个以外的A值的平均作为幅值的估计值。
3、相位检测采用过零法，即通过判断与同频零相位信号过零点时刻，计算其时间差，然后换成相应的相位差。
φ=2п（１－ｔｉ／Ｔ），同样，以φ的平均值作为相位的估计值。
频率、幅值和相位估计的流程图见实验讲义64页图3.13。
4、数字信号统计量估计（1）峰值p的估计在样本数据x中找出最大值与最小值，其差值为双峰值，双峰值的一半即为峰值。
（2）均值估计（3）均方值估计(4) 差估估计5、频谱分析原理时域分析只能反映信号的幅值随时间的变化情况，除单频谱分量的简单波形外，很难明确提示信号的频率组成的各频率分量的大小，而频域分析能很好的解决此问题。
(1)DFT与FFT对于给定的时域信号y，可以通过Fourier变换得到频域信息Y。
采样信号是一个不连续的频谱，不能算出所有点的值，所以采用离散Fourier变换DFT，但其计算效率较低，故实际中常用快速傅立叶变换。
(2)频率、周期的估计对于Y(k△f)，如果当k△f=f~时，Y(k△f)取最大值，则f为频率的估计值，由于采样间隔的误差，f~也存在误差，其误差最大为△f/2。
周期T=1/f，从原理上可以看出，如果在标准信号中混有噪声，用上述方法仍能够精确地估计出原标准信号的频率和周期。
(3)频谱图为了直观地表示信号的频率特性，常常将傅立叶变换的结果用图形的方式表示，及频谱图。
哈尔滨工程大学数字信号处理实验四离散傅里叶变换
哈尔滨工程大学数字信号处理实验四离散傅里叶变换实验四离散傅里叶变换一、实验目的:加深对DFT性质的理解，拓展它们在DSP中的使用;二(实验原理:)对称性，，(2)周期性，， DFT的快速算法FFT利用了的三个固有特性:(1(3)可约性，和。
FFT算法基本上可以分为两大类，即按时间抽选法(DIT,Decimation-In-Time)和按频率抽选法(DIF,Decimation-In-frequency)。
MATLAB中提供了进行快速傅里叶变换的fft函数:X=fft(x)，基2时间抽取FFT算法，x是表示离散信号的向量;X是系数向量;X=fft(x,N),补零或截断的N点DFT，当x得长度小于N时，对补零使其长度为N，当x的长度大于N时，对x截断使其长度为N。
三、实验内容:1.(1) 确定信号的基频w和基本周期Tp，以及分析时采用的采样点数N;w=2pi,Tp=1,N=length(n);(2)当分析长度取0.5Tp和1.5Tp时，对x(t)采样，利用FFT计算其幅度谱，对所得结果进行比较，总结应如何选择分析长度。
0.5T程序>>k=50;>>n=[0:1:24]>>x=cos(10*pi*n/k)+2*sin(18*pi*n/k);>>subplot(211);>>stem(n,x);>>xlabel('n'),ylabel('x(n)');>>xk=fft(x);>>subplot(212);>>stem(n,abs(xk));>>xlabel('k'),ylabel('X(k)');1.5T时程序>>k=50;>>n=[0:1:74]>>x=cos(10*pi*n/k)+2*sin(18*pi*n/k); >>subplot(211);>>stem(n,x);>>xlabel('n'),ylabel('x(n)');>>xk=fft(x);>>subplot(212);>>stem(n,abs(xk));>>xlabel('k'),ylabel('X(k)');jw2、设x(n)=R(n)，分别计算X(e)在[0,2pi]上的32点和64点等间隔采样，并绘制幅频和相8频特性图。
哈工大试验方法和数字信号分析处理作业一
题目:（1）给定数字信号：x(t)=sin(20*pi*t)+sin(100*pi*t)+sin(400*pi*t);即该信号由10HZ,50HZ,200HZ。
三个正弦信号合成。
要求：绘出上述给定数字信号的曲线x(t)。
低通滤波练习：分别用FIR、IIR滤波器滤去50Hz、200Hz信号，保留10Hz信号；绘出滤波前和滤波后的信号曲线，并做对比；滤波过程中的问题讨论。
带通滤波练习：用FIR滤波器滤去10Hz、200Hz信号，保留50Hz信号；绘出滤波前和滤波后的信号曲线，并做对比；滤波过程中的问题讨论。
（2）给定数字信号：X（t）=sin(2*pi*10*t)+sin(2*pi*50*t)+sin(2*pi*200*t)+0.6*randn(1,N)即在原信号上叠加上一个白噪声信号。
要求：绘出上述给定数字信号的曲线x(t)。
分别用低通滤波器和带通滤波器（FIR、IIR任选）滤波、绘曲线对比、讨论。
注：本次作业要求使用我们课上（§3-3、§3-4）所推导的滤波器（公式）滤波；不许使用MATLAB中的滤波函数。
1.数字信号为：x(t)=sin(20*pi*t)+sin(100*pi*t)+sin(400*pi*t);时因为，最大频率为200HZ，故由采样定理dt<=1/2*f max，可得dt<=0.0025s，取dt=0.0003s，满足采样定理。
（1）绘出x(t)图像：Matlab代码：clear allt=0:0.0005:0.6;t1=0.0005;F=15;N=1201;x=sin(2*pi*10*t)+sin(2*pi*50*t)+sin(2*pi*200*t);x1=sin(2*pi*10*t);plot(t,x,'b')；图形如下：图1 原始信号图像（2）低通滤波练习：1.FIR滤波器：Matlab代码：clear allt=0:0.0005:0.6;t1=0.0005;F=15;x=sin(2*pi*10*t)+sin(2*pi*50*t)+sin(2*pi*200*t);x1=sin(2*pi*10*t);y(1201)=0;for k=50:1100for i=-20:20if i==0fi=2*F*t1;elsefi=sin(2*pi*F*i*t1)/pi/i;endy(k)=y(k)+fi*x(k-i);endendplot(t,x1,'k',t,x,'b',t,y,'r');图像如下：图2 FIR低通滤波信号图像图3 FIR低通滤波信号图像i=-30:30,k=70:1100时分析讨论：由图可以看出，原始图像有正弦信号叠加后十分混乱，滤波后基本滤出了10HZ的信号，设计滤波器时，通过改变N1和N2以及采样的数量来生成不同的滤波后图像，最终选择了如上代码中的数值。
哈尔滨工程大学数字信号处理实验五 谱分析
实验五 谱分析一.实验原理信号是无限长的，而在进行信号处理是只能采用有限长信号，所以需要将信号“截断”。
在信号处理中，“截断”被看成是用一个有限长的“窗口”看无限长的信号，或者从分析的角度是无限长的信号乘以有限长的窗函数。
二.实验内容1、用matlab 编程绘制各种窗函数的形状。
2、用matlab 编程绘制各种窗函数的幅频响应。
矩形窗N=20;n=0:(N-1);w=boxcar(N);subplot(211);stem(n,w);title('形状');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('幅频响应');0246810121416182000.51矩形窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.505101520矩形窗幅频响应汉宁窗N=20;n=0:(N-1);w=hanning(N);subplot(211);stem(n,w);title('形状');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('幅频响应');0246810121416182000.51汉宁窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5051015汉宁窗幅频响应汉明窗N=20;n=0:(N-1);w=hamming(N);subplot(211);stem(n,w);title('形状');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('幅频响应');0246810121416182000.51汉明窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5051015汉明窗幅频响应巴特利特窗N=20;n=0:(N-1);w=bartlett(N);subplot(211);stem(n,w);title('形状');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('幅频响应');0246810121416182000.51巴特利特窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50510巴特利特窗幅频响应布莱克曼窗N=20;n=0:(N-1);w=blackman(N);subplot(211);stem(n,w);title('形状');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('幅频响应');0246810121416182000.51布莱克曼窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.502468布莱克曼窗幅频响应Triang 窗N=20;n=0:(N-1);w=triang(N);subplot(211);stem(n,w);title('形状');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('幅频响应');0246810121416182000.51triang 窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50510triang 窗幅频响应Kaiser 窗N=20;n=0:(N-1);w=kaiser(N);subplot(211);stem(n,w);title('kaiser´°ÐÎ×´');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('kaiser´°·ùƵÏìÓ¦');0246810121416182000.51kaiser 窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.505101520kaiser 窗幅频响应切比雪夫窗N=20;n=0:(N-1);w=chebwin(N);subplot(211);stem(n,w);title('Æõ±ÈÑ©·ò´°ÐÎ×´');[H,W]=dtft(w,1024);subplot(212);plot(W/2/pi,abs(H));title('Æõ±ÈÑ©·ò´°·ùƵÏìÓ¦');0246810121416182000.51契比雪夫窗形状-0.5-0.4-0.3-0.2-0.100.10.20.30.40.502468契比雪夫窗幅频响应3、绘制矩形窗的幅频响应，窗长度分别为：10.20,50,100.N=10时N=10;n=0:(N-1);w=boxcar(N);[H,W]=dtft(w,1024);plot(W/2/pi,abs(H));title('¾ØÐδ°·ùƵÏìÓ¦');-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5012345678910矩形窗幅频响应N=20时N=20;n=0:(N-1);w=boxcar(N);[H,W]=dtft(w,1024);plot(W/2/pi,abs(H));title('¾ØÐδ°·ùƵÏìÓ¦');-0.5-0.4-0.3-0.2-0.100.10.20.30.40.502468101214161820矩形窗幅频响应N=50时N=50;n=0:(N-1);w=boxcar(N);[H,W]=dtft(w,1024);plot(W/2/pi,abs(H));title('¾ØÐδ°·ùƵÏìÓ¦');-0.5-0.4-0.3-0.2-0.100.10.20.30.40.505101520253035404550矩形窗幅频响应N=100时N=100;n=0:(N-1);w=boxcar(N);[H,W]=dtft(w,1024);plot(W/2/pi,abs(H));title('¾ØÐδ°·ùƵÏìÓ¦');-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50102030405060708090100矩形窗幅频响应4、已知周期信号，若截取时间长度分别为信号周期的0.9和1.1倍，试绘制和比较采用下面窗函数提取的频谱。
哈工程数字信号处理实验报告1
数字信号处理实验实验一：基本信号班级：姓名：学号：指导教师：2012年10月实验一：基本信号一：实验原理：本节专注于用MATLAB产生一些基本离散信号的问题。
主要是有那个MATLAB内部向量程序来产生信号。
用MATLAB的stem指令会出离散时间信号。
依据MATLAB的编址约定，标号n=0必须对应nn（1）；必须给指定向量的第一个参数以得到正确的n轴。
二：实验内容：1.冲击信号产生并绘出下面的序列。
在每种情况下，水平n轴应该只在指定的区间上展开并应该相应标注。
使用stem指令使每个序列显示成离散时间信号。
x[n]=0.9δ[n-5] 1<=n<=20x[n]=0.8δ[n] -15<=n<=15x[n]=1.5δ[n-333] 300<=n<=350x[n]=4.5δ[n+7] -10<=n<=0L=20;nn=1:(L);imp=zeros(L,1);imp(5)=0.9;stem(nn,imp))L=31;nn=-15:(L-16);imp=zeros(L,1);imp(16)=0.8;stem(nn,imp))L=51;nn=300:350;imp=[zeros(L,1)]'; imp(34)=1.5 stem(nn,imp)L=11;nn=-10:(L-11);imp=zeros(L,1);imp(4)=4.5;stem(nn,imp)实验分析：所得4个图形均符合题目要求3、指数信号衰减的指数信号是数字信号是数字信号处理的基本信号。
因为它是线性常系数差分方程的解。
A.使用函数在区间n=0,1,2，。
，20上绘出信号x[n]=(0.9)ⁿ。
B.在许多推导中，指数信号序列aⁿu[n]须在有限区间上求和。
使用（a）中的函数产生一个指数信号然后对其求和并比较结果。
C.指数序列在信号处理中常常出现的一个原因是，时移并不改变其信号特征。
证明一有限长指数信号满足移位关系：y[n]=ay[n-1], 1<=n<=L-1比较向量y（2：L）和a*y(1:L-1)。
哈工大 数字信号处理实验报告
实验一： 用FFT 作谱分析实验目的：(1) 进一步加深DFT 算法原理和基本性质的理解(因为FFT 只是DFT 的一种快速算法， 所以FFT 的运算结果必然满足DFT 的基本性质)。
(2) 熟悉FFT 算法原理和FFT 子程序的应用。
(3) 学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT 。
实验原理： DFT 的运算量：一次完整的DFT 运算总共需要2N 次复数乘法和(1)N N -复数加法运算，因而直接计算DFT 时，乘法次数和加法次数都和2N 成正比，当N 很大时，运算量很客观的。
例如，当N=8时，DFT 运算需64位复数乘法，当N=1024时，DFT 运算需1048576次复数乘法。
而N 的取值可能会很大，因而寻找运算量的途径是很必要的。
FFT 算法原理：大多数减少离散傅里叶变换运算次数的方法都是基于nk N W 的对称性和周期性。
（1）对称性()*()k N n kn knNN NW W W --==（2）周期性kn Nkn n N kn k NNN NNW W W W ++===由此可得()()/2(/2)1n N k N n k nk N N N N N k N k N N W W W W W W ---+⎧==⎪=-⎨⎪=-⎩这样：1.利用第三个方程的这些特性，DFT 运算中有些项可以合并；2.利用nk N W 的对称性和周期性，可以将长序列的DFT 分解为短序列的DFT 。
前面已经说过，DFT 的运算量是与2N 成正比的，所以N 越小对计算越有利，因而小点数序列的DFT 比大点数序列的DFT 运算量要小。
快速傅里叶变换算法正是基于这样的基本思路而发展起来的，她的算法基本上可分成两大类，即按时间抽取法和按频率抽取法。
我们最常用的是2M N =的情况，该情况下的变换成为基2快速傅里叶变换。
完成一次完整的FFT 计算总共需要2log 2N N次复数乘法运算和2log N N 次复数加法运算。
哈尔滨工程大学数字信号处理实验六
哈尔滨工程大学数字信号处理实验六实验原理IIR数字滤波器一、脉冲响应不变法变换原理将模拟滤波器的s平面变换成数字滤波器的平面，而将模拟滤波器映射成数字滤波器。
MATLAB 信号处理工具箱中提供了IIR数字滤波器设计的函数，常用的函数：IIR滤波器的阶数选择Buttord----巴特沃斯滤波器阶数选择cheb1ord-----切比雪夫I型滤波器阶数选择cheb2ord ----切比雪夫Ⅱ型滤波器阶数选择IIR滤波器的设计Buttrer-----巴特沃斯滤波器设计cheby1-----切比雪夫I型滤波器设计cheby2 ----切比雪夫Ⅱ型滤波器设计maxflat----通用的巴特沃斯低通滤波器设计二、巴特沃斯滤波器设计巴特沃斯滤波器是通带，阻带都单调衰减的滤波器。
(1)调用buttord函数确定巴特沃斯滤波器的阶数，格式[N,Wc]=buttord(Wp,Ws,Ap,As)(2)调用butter函数设计巴特沃斯滤波器，格式[b,a]=butter(N,Wc,options)利用以上两个函数可以设计出模拟滤波器，格式为[N,Wc]=buttord(Wp,Ws,Ap,As,’s’)[b,a]=butter(N,Wc,options.’s’)三、切比雪夫I型滤波器的设计切比雪夫I型滤波器为通带波纹控制器；在通带呈现波纹特性，阻带单调衰减。
[N,Wc]= cheb1ord (Wp,Ws,Ap,As)[b,a]= cheby1 (N,Ap,Wc,options)四、切比雪夫Ⅱ型滤波器的设计切比雪夫Ⅱ型滤波器为阻带波纹控制器；在阻带呈现波纹特性，通带单调衰减。
[N,Wc]= cheb2ord (Wp,Ws,Ap,As)[b,a]= cheby2 (N,As,Wc,options)已知模拟滤波器。
可以利用脉冲响应不变法转换函数impinvar将其变为数字滤波器，调用格式为[bz,az]=impinvar(b,a,Fs).五、双线性变换原理采用非线性频率压缩方法，克服了脉冲响应不变法产生频率响应的混叠失真，使是s平面与z平面建立了一一对应的单值关系，消除多值变换性，频谱混叠现象。
哈工程 数字信号处理实验 实验一
实验一 基本信号一：实验原理使用MATLAB 内部向量程序来产生信号。
用MATLAB 的stem 指令绘出离散时间信号.用MATLAB 的stem 指令会出离散时间信号。
依据MATLAB 的编址约定，标号n=0必须对应nn （1）；必须给指定向量的第一个参数以得到正确的n 轴。
二：实验内容2.正弦信号X 【n 】=Acos （ωn+ρ）使用MATLAB 的向量功能求解此问题，将向量赋予余弦或正弦函数，再利用一个函数调用。
在每种指定区间上展开并标注水平轴n 轴。
使用stem 指令显示每个序列。
A ．X ₁【n 】=sin （πn ／17） 0≤n ≤25 B ．X ₂【n 】=sin （πn ／17） -15≤n ≤25 C ．X ₃【n 】=sin （3πn ﹢π／2） -10≤n ≤10 D ． X ₄【n 】=sin （πn ／√23） 0≤n ≤503.指数信号衰减的指数信号是数字信号处理中的基本信号，因为它是线性常系数方程的解。
A.利用functions 研究下面的MATLAB 函数，看他如何产生离线指数信号。
然后是用函数在区间n=0，1，2，3，. . . ，20上绘出指数信号x 【n 】=（0.9）ⁿ。
B.指数信号序列a ⁿu 【n 】须在有限区间上求和。
这个和以下面闭合时表示: a a a L nL n --=∑-=1110C.指数序列在信号处理中常常出现的一个原因是，时移并不改变其信号特征。
Y【n】=ay【n-1】，1≤n≤L-1D.产生指数信号另外的方法是使用查分方程给出的递归表示式。
当输入x【n】是一个冲击信号的时候，信号y【n】=aⁿu【n】是下面查分方程的解：y【n】-a【n-1】=x【n】，初始条件y【-1】=0三．实验程序2.正弦函数A.n=0:0.01:25y=sin（pi*n／17）plot(n,y)ylabel(‘y= sin（pi*n／17）’)gridgtext(‘n’)B.n=-15:0.01:25y=sin（pi*n／17）plot(n,y)ylabel(‘y= sin（pi*n／17）’)gridgtext(‘n’)C. n=-10:0.01:10y=sin（3*pi*n＋pi／2）ylabel(‘y= sin（3*pi*n＋pi／2）’)gridgtext(‘n’)D. n=0:0.01:50y=cos（pi*n/(23^0.5)）plot(n,y)ylabel(‘cos（pi*n/(23^0.5)’)gridgtext(‘n’)3.指数信号A function y = genexp( b, n0, L)%GENEXP generate an exponential signal:b^n% usage:Y = genexp(B,N0,L)% B input scalar giving ratio between terns% N0 starting index (integer)% L length of geberated signal% Y output signal Y(1:L)if( L <= 0 )error('GENEXP:length not positive')endnn = n0 + [1 : L]'-1; %---vector of indices y =b .^ nn;endn=0:9;x1=genexp(0.9,0,20)stem(n,x1,'b')B. function y=signal(a,n0,L)if(L<=0)error('SIGNAL:length not positive') endnn=n0+[1:L]'-1y=(1-a.^nn)/(1-a)endu1=genexp(0.9,0,21)ss(1)=0;for i=1:19ss(i+1)=u1(i)+ss(i);endss(19)nn=[1:21]'-1;a=0.9y=(1-a.^nn)/(1-a)C. format compact, subplot(111)n=0:20;m=1:21;x1=genexp(0.9,0,21)subplot(211)stem(n,x1,'b')sum(x1(:))grid,title('Test1_3_3_1')x2=0.9*genexp(0.9,0,21)subplot(212)stem(m,x2,'b')grid,title('Test1_3_3_2')四．结果分析1.正弦函数A.BCD3.A.B.x1 =1.00000.90000.81000.72900.65610.59050.53140.47830.43050.38740.34870.31380.28240.25420.22880.20590.18530.16680.15010.13510.1216ans =8.9058用题目所给公式所求结果为8.9058，结果一致。
哈工大数字信号处理大作业 (2)(word文档良心出品)
数字信号处理上机实验报告学号：姓名:实验题目一1. 实验要求：序列卷积计算（1）编写序列基本运算函数，序列相加、相乘、翻转、求和；（2）使用自定义函数计算序列线性卷积，并与直接计算结果相比较。
两个序列分别为：() 1,05 0,others n nx n≤≤⎧=⎨⎩，()2,030,othersn nx n≤≤⎧=⎨⎩2. 实验过程和步骤：包含题目分析，实验程序和流程图（程序要有必要的注释）3. 实验结果和分析：包含程序运行结果图，结果分析和讨论(一)基本运算函数1.原序列2.序列相加序列相加程序function [y,n]=sigadd(x1,n1,x2,n2)%implements y(n)=x1(n)+x2(n)%---------------------------------------------% [y,n] = sigadd(x1,n1,x2,n2)% y = sum sequence over n, which includes n1 and n2% x1 = first sequence over n1% x2 = second sequence over n2 (n2 can be different from n1)%n=min(min(n1),min(n2)):max(max(n1),max(n2)); %duration of y(n) y1=zeros(1,length(n));y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; %x1 with duration of y y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; %x2 with duration of y y=y1+y2; %sequence addition3.序列相乘序列相乘程序function [y,n]=sigmult(x1,n1,x2,n2)%implements y(n)=x1(n)*x2(n)%---------------------------------------------% [y,n] = sigmult(x1,n1,x2,n2)% y = product sequence over n, which includes n1 and n2% x1 = first sequence over n1% x2 = second sequence over n2 (n2 can be different from n1)%n=min(min(n1),min(n2)):m(min(n1),min(n2)) %duration of y(n)y1=zeros(1,length(n));y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; %x1 with duration of y y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; %x2 with duration of y y=y1.*y2; %sequence multiplication4.序列翻转序列翻转程序function [y,n]=sigfold(x, n)%implements y(n)=x(-n)%--------------------------------------------- % [y,n] = sigfold(x,n)%y=fliplr(x);n=-fliplr(n);5.序列移位序列移位程序function [y,n]=sigshift(x,m,n0)%implements y(n)=x(n-n0)%--------------------------------------------- % [y,n] = sigshift(x,m,n0)%n=m+n0;y=x;主程序x1=[0:5];x2=[0,1,2,3];n1=0:5;n2=0:3;%N=n1+n2-1;figure(1)subplot(211)stem(x1)xlabel('x1')subplot(212)stem(x2)xlabel('x2')title('原序列')x= sigadd(x1,n1,x2,n2);figure(2)stem(x)xlabel('x1+x2')title('序列相加')figure(3)[x,n] = sigfold(x1,n1);stem(n,x)xlabel('x1(-n)')title('序列翻转')[x,n] = sigshift(x,n,2);figure(4)stem(n,x)xlabel('x1(-n+2)')title('序列移位')x= sigmult(x1,n1,x2,n2);figure(5)stem(x)title('序列相乘')xlabel('x1*x2')(二)自定义函数计算线性卷积1.题目分析使用上一题中的序列相乘、翻转和求和子函数计算线性卷积，并与这直接用conv 函数计算的线性卷积结果相比较。
哈工大数字信号处理实验报告
哈工大数字信号处理实验报告数字信号处理实验报告姓名:于淼学院:航天学院班级:21201 学号:1112120216 指导教师:李杨1实验一用FFT作谱分析一、实验目的1．进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质)。
2．熟悉FFT算法原理和FFT子程序的应用。
3．学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。
二、实验步骤1．复习DFT的定义、性质和用DFT作谱分析的有关内容。
复习FFT算法原理与编程思想，并对照DIT—FFT运算流图和程序框图， 2．读懂本实验提供的FFT子程序。
3．编制信号产生子程序，产生以下典型信号供谱分析：x1?n??R4?n?n?1， 0?n?3x2?n?? 8?n， 4?n?70 ，其它n 4?n， 0?n?3x3?n?? n?3， 4?n?7，其它nx4?n??cosn4 x5?n??sin??8nx6?t??cos8?t?cos16?t?cos20?t应当注意，如果给出的是连续信号xa?t?，则首先要根据其最高频率确定采样速率fs以及由频率分辨率选择采样点数N，然后对其进行软件采样(即计算x?n??xa?nT?，0?n?N?1)，产生对应序列x?n?。
对信号x6?t?，频率分辨率的选择要以能分辨开其中2的三个频率对应的谱线为准则。
对周期序列，最好截取周期的整数倍进行谱分析，否则有可能产生较大的分析误差。
4．编写主程序下图给出了主程序框图，供参考。
本实验提供FFT子程序和通用绘图子程序。
开始读入长度N 调用信号产生子程序产生实验信号调用绘图子程序(函数)绘制时间序列波形图调用FFT子程序(函数)计算信号的DFT 调用绘图子程序(函数)绘制X?k?曲线结束主程序框图三、实验结果直接运行程序，按照实验内容及程序提示键入1~8，分别对x1?n?~x6?n?及x7?n??x4?n??x5?n?、x8?n??x4?n??jx5?n?进行谱分析。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数字信号处理实验报告Part1 分别用脉冲响应不变法和双线性变换法将图中的RC 低通滤波器转换为数字滤波器，RC=0.001。
1(),a a H s a a s RC==+ 脉冲响应不变法: 1()1aTa H z ez--=-双线性变换法：111121212112(1)()()|12,22a z s T z a z H z H s a z aT aT a a aT aT -----=++==+-==++利用matlab 绘图,分别取T=0.0005s 、T=0.001s 、T=0.002s1. 模拟滤波器的频响曲线()sa as H a +=程序及结果如下： RC=0.001; a=1/RC;f=0:10:1000-10; Hk=a./(a+j*2*pi*f); Hk=abs(Hk); plot(f/1000,Hk); grid;2. 分别用脉冲响应不变法和双线性变换法得到数字滤波器的系统函数H(z)；冲激响应不变法得到函数：1()1aTa H z ez--=-双线性变换法得到函数：11H(z)=(2)(2)aT aTz aT aT z--+++- 3..当采样间隔T=0.0005s 、T=0.001s 、T=0.002s ，分别画出两种数字滤波器的频响曲线，示例程序如下：脉冲法: %模拟RC=0.001; a=1/RC;f=0:10:1000-10; Hk=a./(a+j*2*pi*f); Hk=abs(Hk); f=f/max(f); plot(f,Hk) hold on%数字 a=1000; T=0.0005; B=[a]; A=[1,-exp(-a*T)]; [H,W]=freqz(B,A,100); n=0:99; n=n/max(n); H=abs(H); H=H/max(H); plot(n,H,'r') axis([0 1 0 1]) grid;双线性法：%模拟RC=0.001; a=1/RC; f=0:10:1000; Hk=a./(a+j*2*pi*f); Hk=abs(Hk); f=f/max(f); plot(f,Hk) hold on;%数字a=1000;T=0.0005;B=[a*T a*T];A=[2+a*T ,a*T-2];[H,W]=freqz(B,A,100);n=0:99;n=n/max(n);H=abs(H);H=H/max(H);plot(n,H,'r')axis([0 1 0 1]);grid;4.将模拟滤波器的频响曲线放置在同一坐标中(蓝色曲线)，T不同时需要模拟频率f的取值范围也不相同；曲线如下，其中红色曲线为数字滤波器，蓝色曲线为模拟滤波器。
冲激相应不变法：双线形变换法：比较两种方法设计的特点：在通带范围内，这两种设计方法的响应曲线与模拟滤波器的响应曲线基本吻合。
但在阻带范围内，这两种设计方法的响应曲线与模拟滤波器的响应曲线有较大偏差。
脉冲响应不变法在高频时性能不如模拟滤波器，即阻带衰减慢；而双线性变换法在高频时性能比模拟滤波器好，阻带衰减快！脉冲响应不变法: 频率坐标变换是线形的，由于由s 到z 是多值映射，造成高频混叠。
在w = 附近的高频混叠。
从图中可以看出，减小采样时间可以减小频率混叠。
双线性变换法: 频率坐标变换是非线性的，由s 到z 的映射是单值的，采用非线性频率压缩的方法，因此消除了频谱高频混叠。
但同时也会造成频响曲线形状的偏离。
采样时间越小，频响曲线形状偏离减小。
Part2 用巴特沃斯模拟滤波器设计一个IIR 数字低通滤波器，滤波器的技术要求为：0.2,1,0.3,15,1p p s s rad a dB rad a dB T sωπωπ=====1、利用脉冲响应不变法设计： 源程序为： clearT=1;Wp=0.2*pi/T;Rp=1;Ws=0.3*pi/T;Rs=15;Wc=Wp/((10^(0.1*abs(Rp))-1)^(1/(2*N)));[z,p,k]=buttap(N);[Bp,Ap]=zp2tf(z,p,k);[Bs,As]=lp2lp(Bp,Ap,Wc);[Bz,Az]=impinvar(Bs,As,T);[H,W]=freqz(Bz,Az,100);w=W/pi;h=abs(H);subplot(211);plot(w,h);grid on;hh=20*log10(h);subplot(212);plot(w,hh);grid on;所求得的H（z）为：画出的数字滤波器的幅度响应和增益曲线如下：2、利用双线性变换法设计：程序为：clear;T=1;Wp=2*tan(0.1*pi)/T;Rp=1;Ws=2*tan(0.15*pi)/T;Rs=15;[z,p,k]=buttap(N); [Bp,Ap]=zp2tf(z,p,k); [Bs,As]=lp2lp(Bp,Ap,Wc); [Bz,Az]=bilinear(Bs,As,T); [H,W]=freqz(Bz,Az,100); w=W/pi;h=abs(H); subplot(211); plot(w,h); grid on;axis([0 1 0 1]) hh=20*log10(h); subplot(212); plot(w,hh); grid on;axis([0 1 -80 0])所求得的H （z ）为：1234561234560.00070.00440.01110.01480.01110.00440.0007() 1.0000 3.1836 4.6222 3.7795 1.81360.48000.0544z z z z z z H z z z z z z z ------------++++++=-+-+-+画出的数字滤波器的幅度响应和增益曲线如下：比较两种设计方法的特点：在通带范围内，频响特性曲线无太大区别；在过渡带和阻带内，差别很大。
脉冲响应不变法: 频率坐标变换是线形的，有高频混叠。
双线性变换法: 由于Ω与ω之间的非线性关系，数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变。
Part3IIR和FIR滤波效果的比较，分别用上述实验中的IIR和FIR低通滤波器对输入信号x n=x1+x2+x3进行滤波，滤除高频干扰信号x3，输出得到低频信号y n=x1+x2。
1. IIR 低通滤波器利用上面双线性变换法设计的IIR 低通滤波器对xn进行滤波。
clear;T=1;Wp=(2/T)*tan(0.5*0.2*pi/T);Rp=1;Ws=(2/T)*tan(0.5*0.3*pi/T);Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs,'s');[Bs,As]=butter(N,Wc,'s');[Bz,Az]=bilinear(Bs,As,T);N=127;n=0:N-1;fs=64;x1=sin(pi*n*6.4/fs);x2=sin(pi*n*9.6/fs);x3=sin(pi*n*28.8/fs);xn=x1+x2+x3;yn=filter(Bz,Az,xn);subplot(3,1,1)plot(n,xn);axis([0 N-1 -3 3]);subplot(3,1,2)plot(n,x1+x2);axis([0 N-1 -2 2]);subplot(3,1,3)plot(n,yn);axis([0 N-1 -2 2]);实验结果：Part4 滤波器的实现 利用差分方程1()()()M Ni i i i y n b x n i a y n i ===---∑∑对上面双线性变换法设计的IIR 滤波器，自己编写滤波程序（只允许包含加减乘除运算），对Part4 中的信号xn 滤波，结果应与filter(Bz, Az, xn)函数结果基本一致。
解：利用双线性变换法设计的IIR 滤波器：1234561234560.00070.00440.01110.01480.01110.00440.0007() 1.0000 3.1836 4.6222 3.7795 1.81360.48000.0544z z z z z z H z z z z z z z ------------++++++=-+-+-+差分方程如下：() 3.1816(1) 4.6222(2) 3.7795(3) 1.8136(4)0.4800(5)0.0544(6)0.0007()0.0044(1)0.0111(2)0.0148(3)0.0111(4)0.0044(5)0.0007(6)y n y n y n y n y n y n y n x n x n x n x n x n x n x n =---+---+---++-+-+-+-+-+-源程序如下：clearN=127;n=0:N-1;fs=64;x1=sin(pi*n*6.4/fs);x2=sin(pi*n*9.6/fs);x3=sin(pi*n*28.8/fs);xn=x1+x2+x3;y(1:6)=0;for i=7:127y(i)=0.0007*xn(i)+0.0044*xn(i-1)+0.0111*xn(i-2)+0.0148*xn (i-3)+0.0111*xn(i-4)+0.0044*xn(i-5)+0.0007*xn(i-6)...+3.1836*y(i-1)-4.6222*y(i-2)+3.7795*y(i-3)-1.8136*y(i-4)+ 0.48*y(i-5)-0.0544*y(i-6);endN=length(y);n=1:1:N;plot(n,y);axis([0 N-1 -2 2]);对比实验结果可以看出,实验结果与 part3 完全相同.。