数字信号处理实验

数字信号处理实验报告姓名班级学号实验一信号、系统及系统响应一．实验目的1.掌握线性卷积定义及意义；2.掌握连续卷积与离散卷积的联系，利用离散卷积计算连续信号的卷积结果。二．实验要求1.编写所有程序；2.根据信号系统理论，计算系统的阶跃响应、方波响应，并画出结果；3.编写程序计算系统的阶跃响应、方波响应，并与计算结果进行比较；4.编写实验报告。三．实验内容1.复习采样理论、离散信号与系统、线性卷积有关内容。2.编制实验用主程序及相应子程序，包括（1）方波，阶跃信号产生子程序，体会计算机表示方波，阶跃信号的异同。：closeall;clc;clear;A=1;w0=pi;t=0:0.001:4;y=A*square(w0*t);plot(t,y);title('方波产生信号');axis([-1,5,-1.5,1.5]);阶跃信号子程序：closeall;clc;clear;t=0:0.01:4;plot(t,rectpuls(t,3));title('阶跃信号产生程序');axis([0,3,0,1.5]);（2）系统单位脉冲响应序列产生子程序。本实验要用到两种系统。a、ha(t)=u(t)-u(t-1)产生子程序：closeall;clc;clear;n=-20:20;y1=[zeros(1,20),ones(1,21)];y2=[zeros(1,21),ones(1,20)];subplot(311)stem(n,y1,'fill');xlabel('n');ylabel('u(t)');subplot(312)stem(n,y2,'fill');xlabel('n');ylabel('u(t-1)');subplot(313);y=y1-y2;stem(n,y,'fill');xlabel('n');ylabel('u(t)-u(t-1)');b、ha(t)=e-2tu(t)子程序设计：closeall;clear;clc;t=0:0.01:2;a=exp(-2*t);stem(t,a);axis([0,1,0,1.5]);xlabel('t');ylabel('exp(-2t)u(t)');holdon;（3）有限长序列线性卷积子程序，可以直接调用MATLAB语言中卷积函数conv。closeall;clc;clear;a=[10123];b=[01321];c=conv(a,b);m=length(c)-1;n=0:1:m;stem(n,c,'fill');title('有限序列卷积');gridon;3.通过编程实现方波信号输入系统，观测输出信号，与实际卷积计算结果对比，验证实验结果的正确性。验证方波卷积的程序如下：closeall;clc;clear;Ts=0.2;t=0:Ts:2;a=rectpuls(t,2);b=conv(a,a)*Ts;m=(length(b)-1)*Ts;n=0:Ts:m;t1=0:0.01:2-0.01;c=square(pi*t1,50);t2=0:0.01:1-0.01;y1=t2;y2=1-t2+0.01;y=[y1,y2];subplot(221);stem(t,a);xlabel('t');ylabel('x(t)/h(t)');title('输入方波');subplot(222);stem(n,b);xlabel('t');ylabel('y(t)');title('卷积结果（实际）');subplot(223);plot(t1,c,'r');xlabel('t');ylabel('x(t)/h(t)');title('输入方波');subplot(224);plot(t1,y,'r');axis([0,4,0,1]);xlabel('t');ylabel('y(t)');title('卷积结果（理论）');验证e指数函数卷积的程序如下：closeall;clc;clear;t=0:0.1:10;a=exp(-2*t).*(t>=0);b=[ones(1,100)];c=conv(a,b)*0.1;m=(length(c)-1)*0.1;n=0:0.1:m;t1=0:0.1:10;y1=0.5*(1-exp(-2*t1));t2=10:0.1:20;y2=-0.5*(exp(-20)-exp(-2*(t2-10)));subplot(211);stem(t,a);xlabel('t');ylabel('h(t)');title('系统');axis([0,10,0,1]);subplot(212);stem(n,c);holdon;plot(t1,y1,'r',t2,y2,'r');xlabel('t');ylabel('y(t)');title('卷积结果');axis([0,20,0,1]);实验二用FFT作谱分析下图给出了主程序框图，供参考。本实验提供FFT子程序和通用绘图子程序。开始开始读入长度N读入长度N调用信号产生子程序产生实验信号调用信号产生子程序产生实验信号调用绘图子程序(函数)绘制时间序列波形图调用绘图子程序(函数)绘制时间序列波形图调用FFT子程序(函数)计算信号的DFT调用FFT子程序(函数)计算信号的DFT调用绘图子程序(函数)绘制调用绘图子程序(函数)绘制曲线结束结束主程序框图源程序%x6(n)=cos(8пt)+cos(16пt)+cos(20пt)fs=64kHz,N=16,32,64m=input('FFT点数=');n=0:(m-1);subplot(2,2,1);x6=cos(pi*n/8)+cos(pi*n/4)+cos(pi*n*5/16);stem(n,x6,'.');xlabel('n');ylabel('x6(n)');title('x6(n)=cos(8пt)+cos(16пt)+cos(20пt)的函数');subplot(2,2,2);xa=fft(x6,16);i=0:15;stem(i,abs(xa),'.');xlabel('k');ylabel('x6(n)');title('x6(n)的16点FFT');subplot(2,2,3);xb=fft(x6,32);i=0:31;stem(i,abs(xb),'.');xlabel('k');ylabel('x6(n)');title('x6(n)的32点FFT');subplot(2,2,4);xc=fft(x6,64);i=0:63;stem(i,abs(xc),'.');xlabel('k');ylabel('x6(n)');title('x6(n)的64点FFT');直接运行程序，按照实验内容及程序提示键入1~8，分别对~及、进行谱分析。输出的波形及其8点DFT和16点DFT，的16点、32点和64点采样序列及其DFT。1、及其8点和16点DFT2、及其8点和16点DFT3、及其8点和16点DFT4、的8点和16点波形及其DFT5、的8点和16点波形及其DFT6、的16点、32点和64点采样序列波形及其DFT选7时，=1\*GB3①计算并图示和及其DFT。=2\*GB3②程序自动计算并绘图验证DFT的共轭对称性。当N=16时，，。即为的共轭对称分量，而是的共轭反对称分量。根据DFT的共轭对称性，应有以下结果：的8点和16点波形及其DFT绘出和的模。它们正是图中16点的和。选8时，=1\*GB3①计算并图示和及其DFT。=2\*GB3②程序自动计算并绘图验证DFT的共轭对称性的第二种形式：如果，，则，。其中，。的8点和16点DFT程序计算结果如下：及，正好与图中的16点及相同。及，正好与图中16点的及相同。实验三双线性变换法设计IIR数字滤波器实验目的掌握利用双线性变换法设计IIR数字滤波器的原理和具体方法。深入理解利用双线性变换法设计IIR数字滤波器的优缺点及使用范围。实验内容利用巴特沃思模拟滤波器,通过双线性变换法设计巴特沃思数字滤波器,数字滤波器的技术指标为0.90H采样周期为T=1。程序代码T=2; %设置采样周期为2fs=1/T; %采样频率为周期倒数Wp=0.25*pi/T;Ws=0.35*pi/T; %设置归一化通带和阻带截止频率Ap=20*log10(1/0.9);As=20*log10(1/0.18); %设置通带最大和最小衰减[N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %调用butter函数确定巴特沃斯滤波器阶数[B,A]=butter(N,Wc,'s'); %调用butter函数设计巴特沃斯滤波器W=linspace(0,pi,400*pi); %指定一段频率值hf=freqs(B,A,W); %计算模拟滤波器的幅频响应subplot(2,1,1);plot(W/pi,abs(hf)/abs(hf(1))); %绘出巴特沃斯模拟滤波器的幅频特性曲线gridon;title('巴特沃斯模拟滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');[D,C]=bilinear(B,A,fs); %调用双线性变换法Hz=freqz(D,C,W); %返回频率响应subplot(2,1,2);plot(W/pi,abs(Hz)/abs(Hz(1))); %绘出巴特沃斯数字低通滤波器的幅频特性曲线gridon;title('巴特沃斯数字滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');进行试验输出图像：实验分析通过查workspace可以得出脉冲响应为：

H(滤波器幅度响应：H(ω)=|H滤波器相位响应：θ(ω)=利用巴特沃思模拟滤波器，通过双线性变换法设计数字带阻滤波器，数字滤波器的技术指标为0.90|H(0.90采样周期为T=1。程序代码T=1;%设置采样周期为1fs=1/T;%采样频率为周期倒数wp=[0.25*pi,0.75*pi];ws=[0.35*pi,0.65*pi];Wp=(2/T)*tan(wp/2);Ws=(2/T)*tan(ws/2);%设置归一化通带和阻带截止频率Ap=20*log10(1/0.9);As=20*log10(1/0.18);%设置通带最大和最小衰减[N,Wc]=buttord(Wp,Ws,Ap,As,'s');%调用butter函数确定巴特沃斯滤波器阶数[B,A]=butter(N,Wc,'stop','s');%调用butter函数设计巴特沃斯滤波器W=linspace(0,2*pi,400*pi);%指定一段频率值hf=freqs(B,A,W);%计算模拟滤波器的幅频响应subplot(2,1,1);plot(W/pi,abs(hf)); %绘出巴特沃斯模拟滤波器的幅频特性曲线gridon;title('巴特沃斯模拟滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');[D,C]=bilinear(B,A,fs);%调用双线性变换法Hz=freqz(D,C,W);%返回频率响应subplot(2,1,2);plot(W/pi,abs(Hz));%绘出巴特沃斯数字带阻滤波器的幅频特性曲线gridon;title('巴特沃斯数字滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');进行试验输出图像：实验分析通过查workspace可以得出脉冲响应为：

H(滤波器幅度响应：H(ω)=|H滤波器相位响应：θ(ω)=通过切比雪夫模拟滤波器，利用双线性变换法设计数字高通滤波器，数字滤波器的技术指标为0.90|H(采样周期为T=1，并观察T分别为0.01，0.1，0.25，0.5，0.8时数字滤波器的频率响应。程序代码T=1; %设置采样周期为1fs=1/T; %采样频率为周期倒数Wp=(2/T)*tan(0.75*pi/2);Ws=(2/T)*tan(0.65*pi/2);%设置归一化通带和阻带截止频率Ap=20*log10(1/0.9);As=20*log10(1/0.18);%设置通带最大和最小衰减[N,Wc]=cheb1ord(Wp,Ws,Ap,As,'s');%调用cheb1ord函数确定切比雪夫I型滤波器阶数[B,A]=cheby1(N,Ap,Wc,'high','s');%调用cheby1函数设计切比雪夫I型滤波器W=linspace(0,7*pi,400*pi);%指定一段频率值hf=freqs(B,A,W);%计算模拟滤波器的幅频响应subplot(2,1,1);plot(W/pi,abs(hf)); %绘出切比雪夫I型模拟滤波器的幅频特性曲线gridon;title('切比雪夫I型模拟滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');w=linspace(0,pi,400*pi); %指定一段频率值[D,C]=bilinear(B,A,fs); %调用双线性变换法Hz=freqz(D,C,w);%返回频率响应subplot(2,1,2);plot(w/pi,abs(Hz)); %绘出切比雪夫I型数字带阻滤波器的幅频特性曲线gridon;title('切比雪夫I型数字滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');进行试验输出图像如下：实验分析通过查workspace可以得出脉冲响应为：

H(滤波器幅度响应：H(ω)=|H滤波器相位响应：θ(ω)= 当T分别为0.01，0.1，0.25，0.5，0.8时的数字滤波器的频率响应如下图像所示：其中， T=0.01 “” T=0.1 “-------------” T=0.25 “—·—·—” T=0.5 “————” T=0.8 “” 由图像可以看出，当采样周期从0.01s增长到0.8s时，数字滤波器的幅度响应没有任何变化。故对于用双线性变换法设计的切比雪夫高通数字滤波器，由于事先做了预畸处理，消除了频率映射时引起的非线性关系，所得数字滤波器的幅度响应较为稳定。对比用双线性变换法和脉冲响应不变法得到的数字高通滤波器幅度响应形状可知：双线性变换法不会出现由于高频部分超过折叠频率而混叠到低频部分去的现象。这意味着，模拟滤波器的全部频率响应被压缩于等效的数字频率范围（0，π）之内。这是双线性变换法的主要优点，虽然这是以频率的严重非线性关系为代价的，但因为可以用预畸的方法进行补偿，所以仍可以得到比较理想的数字滤波器幅度响应。实验四用窗函数法设计FIR滤波器一、实验目的1.掌握窗函数法设计FIR滤波器的原理和方法，观察用几种常用窗函数设计的FIR数字滤波器技术指标；2.掌握FIR滤波器的线性相位特性；3.了解各种窗函数对滤波特性的影响。二、实验原理与方法如果所希望的滤波器的理想频率响应函数为Hd(ejω)，则其对应的单位脉冲响应为，用窗函数wN(n)将hd(n)截断，并进行加权处理，得到实际滤波器的单位脉冲响应h(n)=hd(n)wN(n)，其频率响应函数为。如果要求线性相位特性，则h(n)还必须满足。可根据具体情况选择h(n)的长度及对称性。三、实验步骤1.写出理想低通滤波器的传输函数和单位脉冲响应。2.写出用四种窗函数设计的滤波器的单位脉冲响应。3.用窗函数法设计一个线性相位FIR低通滤波器，用理想低通滤波器作为逼近滤波器，截止频率ωc=π/4rad，选择窗函数的长度N＝15，33两种情况。要求在两种窗口长度下，分别求出h(n)，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和阻带衰减；4用其它窗函数(汉宁窗(升余弦窗)、哈明窗(改进的升余弦窗)、布莱克曼窗)设计该滤波器，要求同1；比较四种窗函数对滤波器特性的影响。四、实验用MATLAB函数可以调用MATLAB工具箱函数fir1实现本实验所要求的线性相位FIR-DF的设计，调用一维快速傅立叶变换函数fft来计算滤波器的频率响应函数。fir1是用窗函数法设计线性相位FIRDF的工具箱函数，调用格式如下：hn=fir1(N,wc,‘ftype’,window)fir1实现线性相位FIR滤波器的标准窗函数法设计。hn=fir1(N,wc)可得到6dB截止频率为wc的N阶(单位脉冲响应h(n)长度为N+1)FIR低通滤波器，默认(缺省参数windows)选用hammiing窗。其单位脉冲响应h(n)满足线性相位条件：h(n)=h(N-1-n)其中wc为对π归一化的数字频率，0≤wc≤1。当wc=［wc1,wc2］时，得到的是带通滤波器。hn=fir1(N,wc,’ftype’)可设计高通和带阻滤波器。当ftype=high时，设计高通FIR滤波器；当ftype=stop时，设计带阻FIR滤波器。应当注意，在设计高通和带阻滤波器时，阶数N只能取偶数(h(n)长度N+1为奇数)。不过，当用户将N设置为奇数时，fir1会自动对N加1。hn=fir1(N,wc,window)可以指定窗函数向量window。如果缺省window参数，则fir1默认为hamming窗。可用的其他窗函数有Boxcar,Hanning,Bartlett,Blackman,Kaiser和Chebwin窗。例如： hn=fir1(N,wc,bartlett(N+1))使用Bartlett窗设计； hn=fir1(N,wc,chebwin(N+1,R))使用Chebyshev窗设计。 hn=fir1(N,wc,’ftype’,window)通过选择wc、ftype和window