数字信号处理课设

第一篇：数字信号处理课设

信息科学与工程学院

数字信号处理课程设计实验报告

课题名称： 简单信号滤波演示系统 学生姓名： 学 号： 专业班级： 指导老师： 实验时间： 2014.10.8

目 录

第一章

概述.................................3

第二章

第三章

第四章第五章第六章

1.1 FIR、IIR概述.................................3 1.2题目要求......................................3 设计分析.............................5 2.1算法分析......................................5 2.2 在matlab中实现的分析........................6

程序实现.......................8 3.1 程序主体介绍..................................8 3.2 子程序........................................9 3.3 程序调试及运行结果............................9 3.4 结果分析及问题分析...........................16 心得体会............................17 参考文献............................18 源代码..............................19

MATLAB

第一章 概述

1.1 FIR、IIR概述

数字滤波器是指输入输出均为数字信号，通过一定的运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。数字滤波器与模拟滤波器相比数字滤波器具有精度高、稳定、体积小、重量小、灵活等特点。主要分为两种：有限脉冲响应FIR和无限脉冲响应IIR。设计滤波器的主要要求有两种，一是幅频特性，一是相频特性。一般的滤波器主要是对幅频特性作出要求，如果对输出相频特性也有要求，就需要用到线性相位滤波器。IIR滤波器的设计主要有两类，一是借助于模拟滤波器设计进行，二是直接在频域或时域中进行设计。FIR滤波器的设计不能借助于模拟滤波器，也有两类设计方法，一是窗函数法，二是频率采样法。还有一种比较有效的方法是切比雪夫等波纹逼近法，需通过计算机辅助进行。

1.2题目要求

设计一个工作流程如图所示的信号滤波演示系统：

图2 滤波演示图

⑴ 信号发生器—根据信号选择分为两大类：

① 静态型：直接输入(或从文件读取)测试信号序列；

② 动态型：输入由多个不同频率正弦信号叠加组合而成的模拟信号公式 100sin(2f1t)100sin(2f2t)...100sin(2fnt)采样频率（Hz）以及采样点数，动态生成该信号的采样序列，作为测试信号。

⑵ 频谱分析—使用FFT 对产生的测试信号进行频谱分析并展示其幅频、相频特

性，指定需要滤除或保留的频带，通过选择滤波器类型（IIR/FIR），确定对

应的滤波器（低通、高通、带通、带阻）技术指标。

⑶ 滤波器设计—根据IIR/FIR 数字滤波器技术指标设计滤波器，生成相应的滤

波器系数，并展示对应的滤波器幅频（衰减）、相频特性。

① IIR DF 设计：使用双线性变换法，可选择滤波器基型（巴特沃斯或切比雪夫型）；

② FIR DF 设计：使用窗口法，可选择窗口类型。

⑷ 数字滤波—根据设计的滤波器系数，对测试信号进行滤波。

① IIR DF：要求通过差分方程迭代实现滤波，未知初值置零处理； ② FIR DF：要求通过快速卷积实现滤波。可以选择使用重叠相加或重叠保留法进行卷积运算，并动态展示卷积运算的详细过程。

⑸ 输出信号分析—展示滤波后信号的幅频与相频特性，分析是否满足滤波要求。

对同一滤波要求，根据输出信号频谱，对比分析各类滤波器的差异。选做部分

将一段语音作为测试信号，通过频谱展示和语音播放，对比分析滤波前后语音信号的变化，进一步加深对数字信号处理的理解。

第二章 设计分析

2.1算法分析

此题目的实现可分为三个某块的设计实现：输入信号模块，设计滤波器模块，滤波模块。

首先明确各模块间的数据依赖关系：在输入信号模块得到信号后，对信号

进行频域分析，从而确定滤波器的相关技术指标，根据滤波器的技术指标与类型，在滤波器设计模块完成滤波器的设计，然后将滤波器的设计结果传递给滤波模块，滤波模块同时接收输入信号，从而通过运算，实现信号的滤波处理。从数据传递关系上分析，滤波模块的实现依赖于其他模块的数据输出，因此放在最后设计。先设计输入模块。

因为此设计相对复杂，分模块设计，通过参数传递和接口实现分模块设计即检验，提高程序的稳定性与健壮性。

输入的实现可以有两种方式：静态输入和动态输入。静态输入选择从文本输入数据，将信号取样值以矩阵的形式存放在文本中。采用文件的读取就可以实现，比较容易。动态输入，即输入由一系列频率的正弦信号相加的组成的信号，需要经过采样的，注意在设置采样频率时一定要符合奈奎斯特准则，提高采样点数，增加频谱分辨率。最后输出一采样信号向量，传递给其余两模块。

滤波器的设计，通过输入信号的频谱分析，设置滤波器的参数，然后才可以设计滤波器。第一步需要总结设计滤波器需要哪些参数，通过学习可以总结，所有滤波器的参数有四个：通带截止频率、阻带截止频率、通带最大衰减、阻带最小衰减。

对滤波器的设计分两类：FIR和IIR，二则所需的函数及设计方法不同。IIR采用借助于模拟滤波器的方式，包括巴特沃斯滤波器和切比雪夫滤波器两种类型。FIR采用窗函数方式，有矩形窗、三角窗、汉明窗、汉宁窗、布莱克曼和凯塞窗。通过调用不同的函数来实现滤波器的设计。特别在实现窗函数滤波器时，各个函数的主要区别是不同的频率采样，可以通过选择结构实现，简化程序。通过滤波器的设计最后可以得到滤波器的系统函数H(z)的系数。分析滤波器的幅频特性和相频特性，如果不符合要求重新设定滤波器参数或者换成其他滤波器类型。如果性能符合要求，则将系数传递给滤波模块。滤波模块不调用滤波函数，实现滤波功能根据滤波器类型的不同，有两种方式可以选择，一种是通过差分方程运算，一种是通过线性卷积运算。前者适合对IIR滤波器进行滤波，后者适合对FIR滤波器进行滤波。且线性卷积为实现对长序列的卷积运算，采用重叠相加法，采用动态变化展示运算过程。

2.2 在matlab中实现的分析

输入模块通过读取文件和直接输入数据运算可以很容易实现。在输入信号确定后，要对其进行频谱分析，从而确定滤波器的参数和类型（低通、高通、带通、带阻），此模块作用也就完成，将数据分别用全局变量传递给下面的模块。

设置模拟信号：

我采用的测试信号是两个正弦信号叠加而成的信号，信号为

y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)其中频率f1=30;频率f2=50;频率f3=200;采样频率fs=3000;采样间隔 Ts=1/fs;采样点数

N=1024;n=1:N 采集模拟信号的程序代码： f1=30;% 频率1

f2=50;% 频率2 f3=200;% 频率3

fs=3000;% 采样频率

Ts=1/fs;% 采样间隔

N=1024;% 采样点数

n=1:N;

y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)+sin(2*pi*f1*n*Ts）;% 正弦波混合

频谱分析：

使用FFT对产生的测试信号进行频谱分析并展示其幅频特性与相频特性，指定需要滤除的频带，通过选择滤波器类型IIR，确定对应的滤波器（低通、高通）技术指标Fp、Fc、Rp、Rs。

滤波器的设计：

根据以上技术指标（通带截止频率、通带最大衰减、阻带截止频率、阻带最

小衰减），设计数字滤波器，生成相应的滤波器系数，并画出对应的滤波器

幅频特性与相频特性。分别设计巴特沃斯滤波器、切比雪夫I型滤波器、切

比雪夫II型滤波器、椭圆滤波器、yulewalk滤波器。

巴特沃斯和切比雪夫的滤波函数调用为：

[N,Wc]=buttord(wp,ws,rp,rs);[N,Wc]=cheb1ord(wp,ws,rp,rs);[B,A]=butter(N,Wc,’property’);[B,A]=cheby1(N,rp,Wc,’property’);

property对于低通和高通为’’，带通’high’,带阻’stop’;窗函数滤波器设计的调用函数：

求窗函数的阶数：

N=ceil((h*pi)/wdel);%wdel为窗函数的过渡带宽，h对应不同窗函数的值 wn=boxcar(N+1);wn=triang(N+1);wn=hanning(N+1);wn=hamming(N+1);wn=blackman(N+1);wn=kaiser(N+1,beta);%bata为kaiser的a参数 B=fir1(N,ws,'property',wn);property对于低通和高通为’’，带通’high’,带阻’stop’;

数字滤波：

根据设计的滤波器系数，对测试信号通过差分方程迭代实现滤波数字滤波，展示滤波后信号的幅频特性与相频特性，分析是否满足滤波要求（对同一滤波要求，对比分析各类滤波器的差异）。

需要注意的是，窗函数对滤波参数的使用，需要通过运算得到一窗函数阶数，对通带最大衰减无使用。需要注意FIR和IIR对参数的不同利用。

滤波模块：运用差分方程运算的滤波函数，可以用循环调用简单实现。动态展示线性卷积的重叠相加法可以用流程图来说明： 第三章 MATLAB程序实现

3.1 程序主体介绍

此程序界面也是采用GUI，不过是采用图形用户界面开发环境。和第一个对比，布局更加容易，很容易调整界面。但是在一些函数操作行会受一些限制，比如函数参数的传递，目前来说还是采用全局变量才能完成。

pushbutton1_Callback(hObject, eventdata, handles)%获得输入信号

uipanel11_SelectionChangeFcn(hObject, eventdata, handles)%获取滤波器的设计参数

uipanel15_SelectionChangeFcn(hObject, eventdata, handles)%设计滤波器

pushbutton4_Callback(hObject, eventdata, handles)%进行滤波

以上是一些功能控件，需要注意的是他们的view ballcak属性，对于groupbutton一组按钮，需要选用SelectionChangeFcn，其余选择callback属性即可。另外，handles函数作用为传递的句柄，是一结构体，调用时注意此处。其余和直接用函数创建GUI界面无太大区别。

3.2 子程序

3.2.1

IIR滤波器

FIR滤波器

巴特沃斯滤波器

3.2.2

矩形窗

三角窗

汉宁窗

哈明窗 3.3 程序调试及运行结果

原始信号幅频特性及相频特性

巴特沃斯的幅频特性及相频特性

信号滤波后的幅频特性及相频特性

矩形窗的幅频特性及相频特性

信号滤波后的幅频特性及相频特性

三角窗的幅频特性及相频特性

信号滤波后的幅频特性及相频特性

汉宁窗的幅频特性及相频特性

信号滤波后的幅频特性及相频特性

哈明窗的幅频特性及相频特性

信号滤波后的幅频特性及相频特性

3.4 结果分析及问题分析

用双线性变换法设计无限脉冲响应数字滤波器（IIF DF）时，先把数字滤波器指标转换成模拟滤波器的指标，然后根据模拟滤波器的指标设计模拟滤波器，再经过线性变换把模拟滤波器转换成数字滤波器。该系统要能够设计巴特沃兹型低通、带通、高通滤波器，并能够输入数字滤波器的性能指标，显示出滤波器的阶数和系数。该系统的关键部分是滤波器的设计部分，按照双线性变换法设计滤波器的步骤进行设计即可。

第四章 心得体会

本次课程设计可以称为matlab和数字信号处理的综合设计，因为这次有一半的时间在研究matlab，有一半的时间在思考数字信号问题的解决。几天下来收获是很多的，对数字信号处理有了一次全面的回顾，而且也看到了matlab的神秘面容。只要是自己真正的努力去做了，就绝对不会后悔在课程设计上花费的时间。现在发现自己确实会的太少，而这次又学了太多，对自己以后的学习还是有鞭策意义的。

这次课程设计也使我认识到了matlab的强大，或者进一步说工具软件有着你意想不到的功效，能节省你的时间，同时又能够让你从实践上更深层次的认识到所学知识的奇妙，而且可以同时明白了学习与实践的相辅相成。对课程设计也是一直保持很高兴趣的，虽然有时为它焦头烂额，但是也会因此有柳暗花明的喜悦，任何事情都折射着一个普通的道理——付出才有回报。浅显的道理，却是需要我们用毅力来坚持见证的。

具体来说，通过本次课程设计，我掌握了MATLAB的基本运用，尤其是其中GUI可视化图形用户界面的设计，包括函数调用与在图形用户界面开发环境下的调用。函数的调用与参数的传递是两个简单却很容易出错的地方，自由细心才可以保证程序的健壮性。

对数字信号本身内容的理解来说，全书的内容确实是可以融合在这两个课程设计题目中的，一个是DFT运用，另一个是滤波器的设计和利用。对全书的内容可以做出最好的概括。其中滤波器的设计中，调用了许多滤波函数，这是本次实验有点欠缺的地方，但是仍然从整体上把握了整个滤波器的设计过程。

对于课程设计中出现的问题，解决的方式有两种，一是自己解决。可以上网，查阅图书和matlab的本身的帮助文件，不断尝试，自己修改错误，可以更好的反省。二是与同学相商。一加一不是等于二那么简单的，相互交流才是进步的最好方式。其中解决问题最重要的应该是坚持不懈，在遇到问题时，不放弃才有可能称为最后的胜利者。

每次课程设计都收获颇多，而且每次都对自己的学习态度做一次调整，这个的作用确实是很大的，希望以后更加珍惜的课程设计的时间。

第五章 参考文献

[1] 丁玉美等.数字信号处理 [M].西安:西安电子科技大学出版社,2002

[2] 程佩青.数字信号处理教程,第二版[M].北京:清华大学出版社,2001

[3] 赵树杰等.数字信号处理[M].西安:西安电子科技大学出版社,1997

[4] 陈怀琛等.MATLAB 及在电子信息课程中的应用[M],北京:电子工业出版社出

版,2002

[5] 余成波.数字信号处理及MATLAB实现，第一版[M].北京：清华大学出版社，2008

[6] S.K.Mitra.Digital Signal Processing: A Computer Based Approach, 3rdEdition[M],New York, USA:McGraw-Hill,2000

[7] R.G.Lyons.Understanding Digital Signal Processing,2nd Edition[M].New Jersey, USA:Prentice Hall,2005

第六章 源代码

function varargout = df(varargin)

gui_Singleton = 1;

gui_State = struct('gui_Name', mfilename,...'gui_Singleton', gui_Singleton,...'gui_OpeningFcn', @df_OpeningFcn,...'gui_OutputFcn', @df_OutputFcn,...'gui_LayoutFcn', [] ,...'gui_Callback', []);if nargin && ischar(varargin{1})

gui_State.gui_Callback = str2func(varargin{1});end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});else

gui_mainfcn(gui_State, varargin{:});end

function df_OpeningFcn(hObject, eventdata, handles, varargin)

handles.output = hObject;

guidata(hObject, handles);

function varargout = df_OutputFcn(hObject, eventdata, handles)

varargout{1} = handles.output;

function edit1_Callback(hObject, eventdata, handles)

function edit1_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit2_Callback(hObject, eventdata, handles)

function edit2_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit3_Callback(hObject, eventdata, handles)

function edit3_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit4_Callback(hObject, eventdata, handles)

function edit4_CreateFcn(hObject, eventdata, handles)if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit5_Callback(hObject, eventdata, handles)

function edit5_CreateFcn(hObject, eventdata, handles)if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function pushbutton1_Callback(hObject, eventdata, handles)f=zeros(1,4);

f(1)=str2num(get(handles.edit2,'String'));f(2)=str2num(get(handles.edit1,'String'));f(3)=str2num(get(handles.edit3,'String'));f(4)=str2num(get(handles.edit4,'String'));

T=str2num(get(handles.edit12,'String'));%采样周期% fre=str2num(get(handles.edit5,'String'));%采样频率% t=0:1/fre:T;

xn=10*sin(2*pi*t*f(1))+10*sin(2*pi*t*f(2))+10*sin(2*pi*t*f(3))+10*sin(2*pi*t*f(4));

set(handles.axes15,'userdata',xn);%将Xn放在用户数据userdata% yn3=abs(fft(xn));%快速傅立叶变换（符频特性）% n1=[0:length(yn3)-1]/length(yn3)*fre;%横坐标% axes(handles.axes12);%坐标系编号% stem(n1,yn3,'.');

axis([0,fre/2,0,max(yn3)]);%坐标轴单位控制% title('信号的幅频特性');xlabel('频率（Hz）');ylabel('|X(ejw)|');

yn4=angle(fft(xn));%相频特性%

n4=[0:length(yn4)-1]/length(yn4)*fre;axes(handles.axes16);stem(n4,yn4,'.');

axis([0,fre/2,min(yn4),max(yn4)]);title('信号的相频特性');xlabel('频率（Hz）');ylabel('相位');

function edit6_Callback(hObject, eventdata, handles)

function edit6_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit7_Callback(hObject, eventdata, handles)

function edit7_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit8_Callback(hObject, eventdata, handles)

function edit8_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit9_Callback(hObject, eventdata, handles)

function edit9_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function pushbutton2_Callback(hObject, eventdata, handles)xn=get(handles.axes15,'userdata');fre=str2num(get(handles.edit5,'String'));

fp=str2num(get(handles.edit6,'String'));%通带最大频率 % fs=str2num(get(handles.edit7,'String'));%阻带最小频率% rp=str2num(get(handles.edit8,'String'));rs=str2num(get(handles.edit9,'String'));wp=2*fp/fre;ws=2*fs/fre;

[N,wc]=buttord(wp,ws,rp,rs,'s');%求阶数，截止频率%

[B,A]=butter(N,wc,'s');%巴特沃兹模拟低通滤波器系数%

[Bz,Az]=bilinear(B,A,1);

[H,w]=freqz(Bz,Az);%分析数字滤波器% axes(handles.axes14);plot(w/pi,abs(H));

axis([0,1,0,max(abs(H))]);title('巴特沃兹的幅频特性');xlabel('w/π');

ylabel('|X(ejw)|');axes(handles.axes17);plot(w/pi,angle(H));title('巴特沃兹的相频特性');xlabel('频率（Hz）');ylabel('相位');

yn2=filter(Bz,Az,xn);%迭代法求解滤波信号% yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信号滤波后的幅频特性');xlabel('频率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信号滤波后的相频特性');xlabel('频率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function edit10_Callback(hObject, eventdata, handles)function edit10_CreateFcn(hObject, eventdata, handles)

if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function edit11_Callback(hObject, eventdata, handles)

function edit11_CreateFcn(hObject, eventdata, handles)if ispc

set(hObject,'BackgroundColor','white');else

set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));end

function pushbutton5_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;

Bt=ws-wp;

N0=ceil(1.8*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,boxcar(N));%单位序列相应%

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('矩形窗的损耗函数');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);

n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('矩形窗的相频特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);%重叠相加法求滤波后序列% yn=abs(fft(yn2));n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信号滤波后的幅频特性');xlabel('频率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信号滤波后的相频特性');xlabel('频率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function pushbutton6_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;

Bt=ws-wp;

N0=ceil(6.1*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,bartlett(N));

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('三角窗的损耗函数');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('三角窗的相频特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信号滤波后的幅频特性');xlabel('频率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信号滤波后的相频特性');xlabel('频率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function pushbutton7_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;if(fp

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,hanning(N));else

Bt=wp-ws;

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,'high',hanning(N));end

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('汉宁窗的损耗函数');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);

n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('汉宁窗的相频特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信号滤波后的幅频特性');xlabel('频率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信号滤波后的相频特性');xlabel('频率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

function pushbutton8_Callback(hObject, eventdata, handles)

xn=get(handles.axes15,'userdata');fp=str2num(get(handles.edit10,'String'));fs=str2num(get(handles.edit11,'String'));fre=str2num(get(handles.edit5,'String'));wp=2*fp/fre*pi;ws=2*fs/fre*pi;if(fp

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,hamming(N));else

Bt=wp-ws;

N0=ceil(6.2*pi/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2/pi;

hn=fir1(N-1,wc,'high',hamming(N));end

yn=abs(fft(hn));yn=20*log10(yn);

n1=[0:length(yn)-1]/length(yn)*2;axes(handles.axes14);plot(n1,yn);

title('哈明窗的损耗函数');xlabel('w/π');

ylabel('20log|Hg(w)|');axis([0,1,min(yn),max(yn)]);axes(handles.axes17);yn=angle(hn);

n1=[0:length(yn)-1]/length(yn)*2;plot(n1,yn);

title('哈明窗的相频特性');xlabel('w/π');

ylabel('相位');

axis([0,1,min(yn),max(yn)]);

yn2=fftfilt(hn,xn,512);yn=abs(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes15);stem(n1,yn,'.');

title('信号滤波后的幅频特性');xlabel('频率（Hz）');ylabel('|X(ejw)|');axis([0,fre/2,0,max(yn)]);yn=angle(fft(yn2));

n1=[0:length(yn)-1]/length(yn)*fre;axes(handles.axes18);stem(n1,yn,'.');

title('信号滤波后的相频特性');xlabel('频率（Hz）');

ylabel('相位');

axis([0,fre/2,min(yn),max(yn)]);

pushbutton1_Callback(hObject, eventdata, handles);

%---Executes during object creation, after setting all properties.function axes12_CreateFcn(hObject, eventdata, handles)% hObject handle to axes12(see GCBO)

% eventdata reservedhandles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes12

%---Executes on mouse press over axes background.function axes12_ButtonDownFcn(hObject, eventdata, handles)% hObject handle to axes12(see GCBO)

% eventdata reservedto be defined in a future version of MATLAB % handles structure with handles and user data(see GUIDATA)

% Hints: get(hObject,'String')returns contents of edit12 as text % str2double(get(hObject,'String'))returns contents of edit12 as a double

%---Executes during object creation, after setting all properties.function edit12_CreateFcn(hObject, eventdata, handles)% hObject handle to edit12(see GCBO)

% eventdata reservedhandles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.% See ISPC and COMPUTER.if ispc

&&

isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');end

%---If Enable == 'on', executes on mouse press in 5 pixel border.%---Otherwise, executes on mouse press in 5 pixel border or over edit5.function edit5_ButtonDownFcn(hObject, eventdata, handles)% hObject handle to edit5(see GCBO)

% eventdata reservedto be defined in a future version of MATLAB % handles structure with handles and user data(see GUIDATA)

%---Executes during object deletion, before destroying properties.function edit12_DeleteFcn(hObject, eventdata, handles)% hObject handle to edit12(see GCBO)

% eventdata reserved-to be defined in a future version of MATLAB % handles structure with handles and user data(see GUIDATA)

第二篇：数字信号处理知识总结课案（模版）

1.傅里叶变换

x(n)的一个周期； 有限长序列

x(n)

可看成周期序列

~

x(n)看成x把 ~(n)的以N为周期的周期延拓。有限长序列的离散傅里叶变换（DFT）：

N1nk X(k)DFT[x(n)]x(n)WNx(n)x((n))N~n0 ~1N1nkx(n)x(n)RN(n)x((n))NRN(n)x(n)IDFT[X(k)]X(k)W NNk0

① 长度为N的有限长序列 x(n)，其离散傅里叶变换 X(k)仍是一个长度为N 的有限长序列；

② x(n)与X(k)是一个有限长序列离散傅里叶变换对，已知x(n)就能唯一地确定 X(k)；同样已知X(k)也就唯一地确定x(n)。实际上x(n)与 X(k)都是长度为 N 的序列（复序列）都有N个独立值，因而具有等量的信息； ③ 有限长序列隐含着周期性。

2.循环卷积（有可能会让画出卷积过程或结果）

循环卷积过程为：

最后结果为：

3.（见课本）

课本

3、线性卷积（有可能会让画出卷积过程或结果）

以下为PPT上的相关题目：

4.计算分段卷积：重叠相加法和重叠保留法（一定会考一种）

重叠相加法解题基本步骤:

将长序列均匀分段，每段长度为M；

基于DFT快速卷积法，通过循环卷积求每一段的线性卷积；

依次将相邻两段的卷积的N-1个重叠点相加，得到最终的卷积结果。

4.级联、并联、直接形（画图）以下为课后作业相关题目：

1.已知系统用下面差分方程描述: 311y(n)＝y(n1)－y(n2)＋x(n)x(n1)483

试分别画出系统的直接型、级联型和并联型结构。式中x(n)和y(n)分别表示系统的输入和输出信号。

解： 将原式移项得

311y(n)y(n1)y(n2)x(n)x(n1)483

将上式进行Z变换，得到

311121Y(z)Y(z)zY(z)zX(z)X(z)z483

11z13H(z)311z1z248(1)按照系统函数

H(z)，根据Masson公式，画出直接型结构如题1解图

（一）所示。1z11(2)将H(z)的分母进行因式分解：

H(z)  3111z11z1

24111z按照上式可以有两种级联型结构: 13H(z)  111z11z1 241画出级联型结构如题1解图

（二）（a)所示 1z11 H(z)  3111z11z1

24画出级联型结构如题1解图

（二）（b)所示

111z11z133H(z)311211111zz(1z)(1z)4824

(3)

将H(z)进行部分分式展开：

1zH(z)AB3H(z)11111111z(z)(z)zz(1z)(1z)2424 2411z***zzH(z)33333H(z)311111111z(z)z1z1zzz24242410

根据上式画出并联型结构如题1解图

（三）所示。

3.设系统的差分方程为

y(n)=(a+b)y(n－1)－aby(n－2)+x(n－2)+(a+b)x(n－1)+ab 式中, |a|<1，|b|<1，x(n)和y(n)分别表示系统的输入和输出信号, 试画出系统的直接型和级联型结构。

解：

(1)直接型结构。将差分方程进行Z变换，得到

Y(z)=(a+b)Y(z)z－1－abY(z)z－2+X(z)z－2－(a+b)X(z)z－1+ab Y(z)ab(ab)zzH(z)X(z)1(ab)z1abz2按照Masson公式画出直接型结构如题3解图

（一）所示。

12

(2)级联型结构。将H(z)的分子和分母进行因式分解，得到

(az1)(bz1)H(z)H1(z)H2(z)11(1az)(1bz)

按照上式可以有两种级联型结构:①

z1aH1(z)1az1画出级联型结构如题3解图

（二）(a)所示

z1bH2(z)1bz11

zaH1(z)11bz画出级联型结构如题3解图

（二）(b)所示 1zbH2(z)11az

四．设计模拟滤波器（考试时不能编代码）一般步骤：

根据Ap、As、Ωs、Ωp，确定滤波器阶次N和截止频率Ωc。

P161 【例6.2.2】

设计一个模拟低通巴特沃斯滤波器，指标如下：

(1)通带截止频率：Ωp=0.2π；通带最大衰减：Ap=7 dB。

(2)阻带截止频率：Ωs=0.3π；阻带最小衰减：As=16dB。解： lg[(100.71)/(101.61)] N2.7932lg(0.2/0.3)

由Ωp，得:

0.2

Qc0.49850.76 101由Ωs，得：

0.3Q0.5122 c1.66101

在上面两个Ωc之间选Ωc=0.5。

最后可得（级联型）：

0.125Ha(s)

(s0.5)(s20.5s0.25)

五、脉冲响应不变法（P177 第6.3节）156-158页

脉冲响应不变法的优点：

 时域逼近。使数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，即时域逼近良好。 线性频率关系。

模拟频率Ω和数字频率ω之间呈线性关系ω=ΩT。脉冲响应不变法的缺点：  混叠失真效应

因此，只适用于限带的模拟滤波器（例如衰减特性很好的低通或带通滤波器），而且高频衰减越快，混叠效应越小；而对于高通和带阻滤波器，由于它们在高频部分不衰减，因此会产生混叠现象。

六、双线性变换法

七，与实验相关

本题中老师会给出类似于下列表达式的信号：

S(t)23cos(250t30/180)1.5cos(275t90/180)

要求用脉冲相应不变法或双线性法编写主要的代码（如下面代码）来达到滤除其中的部分信号，并画出你所设计的滤波器的频响曲线，并标明Ωs、Ωp，以及滤波后信号的时域波形（波形中要体现相位特征）。1）脉冲响应不变法滤除第三个信号： Fs=256;

% 采样频率 fp=60;

% 通带截止频率 fs=70;

% 阻带截止频率 Rp=1;

Rs=25;Wp=(fp/Fs)*2*pi;

%临界频率采用角频率表示 Ws=(fs/Fs)*2*pi;

%临界频率采用角频率表示 OmegaP=Wp*Fs;OmegaS=Ws*Fs;[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');[Bz,Az]=impinvar(b,a,Fs);

2）双线性法滤除第三个信号： Fs=256;

% 采样频率 fp=60;

% 通带截止频率 fs=70;

% 阻带截止频率 Rp=1;

Rs=25;Wp=(fp/Fs)*2*pi;

% 临界频率采用角频率表示 Ws=(fs/Fs)*2*pi;

% 临界频率采用角频率表示 OmegaP=2*Fs*tan(Wp/2);

% 频率预畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');

[Bz,Az]=bilinear(b,a,Fs);注：要好好看实验中关于低通，高通，带通，带阻的设计代码。

带通：

fp1=40;

% 通带截止频率 fs1=30;

% 阻带截止频率 fp2=60;

% 通带截止频率 fs2=70;

% 阻带截止频率 Rp=1;

Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=[Wp1,Wp2];

% 向量 Ws=[Ws1,Ws2];

% 向量

带阻：

fp1=30;

% 通带截止频率 fs1=40;

% 阻带截止频率 fp2=70;

% 通带截止频率 fs2=60;

% 阻带截止频率 Rp=1;

Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];若S(5t)信号表达

i(n2

式为

t2703ctosts则相关代码为：

1)低通滤波器代码

fp=110;% 通带截止频率 fs=130;% 阻带截止频率

Rp=1;Rs=25;Wp=(fp/Fs)*2*pi;Ws=(fs/Fs)*2*pi;%临界频率采用角频率表示

(1):脉冲响应不变法

OmegaP=Wp*Fs;OmegaS=Ws*Fs;[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');% 指明为高通滤波器 [Bz,Az]=impinvar(b,a,Fs);

（2）双线性变换法

OmegaP=2*Fs*tan(Wp/2);OmegaS=2*Fs*tan(Ws/2);% 频率预畸 [n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(n,Wc,'s');[Bz,Az]=bilinear(b,a,Fs);

2)高通滤波器

fp=280;% 通带截止频率 fs=260;% 阻带截止频率 Rp=1;Rs=25;Wp=(fp/Fs)*2*pi;%临界频率采用角频率表示 Ws=(fs/Fs)*2*pi;%临界频率采用角频率表示(2):双线性变换法

OmegaP=2*Fs*tan(Wp/2);% 频率预畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(2*n,Wc,'high','s');[Bz, Az]=bilinear(b,a,Fs);3)带通滤波器

fp1=130;% 通带截止频率 fs1=110;% 阻带截止频率 fp2=255;% 通带截止频率 fs2=265;% 阻带截止频率 Rp=1;Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];(2):双线性变换法

OmegaP=2*Fs*tan(Wp/2);% 频率预畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(2*n,Wc,'s');[Bz,Az]=bilinear(b,a,Fs);4)带阻滤波器的代码如下：

fp1=110;% 通带截止频率 fs1=240;% 阻带截止频率 fp2=265;% 通带截止频率 fs2=255;% 阻带截止频率 Rp=1;Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;

Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];(2):双线性变换法

OmegaP=2*Fs*tan(Wp/2);% 频率预畸 OmegaS=2*Fs*tan(Ws/2);[n,Wc]=buttord(OmegaP,OmegaS,Rp,Rs,'s');[b,a]=butter(2*n,Wc,'stop','s');[Bz,Az]=bilinear(b,a,Fs);

第三篇：基于微课的数字信号处理教学方法的探讨

基于微课的数字信号处理教学方法的探讨

【摘 要】基于目前数字信号处理课程在教学方法及教学手段上的不足，本文提出了基于“微课”的教学方法，通过对该课程重点难点内容做成微课视频的形式，供学生课前或课后学习，以提高学生的学习兴趣及自主学习能力，再结合双向互动式研究型教学方法及教学实践力度的加大，从而改善教学效果，使独立院校学生能更好的掌握这门课程。

【关键词】数字信号处理；微课；教学改革；实践教学

0 引言

数字信号处理作为通信、电子信息类专业的一门专业基础课程，要求学生有微积分、信号与系统分析等课程的基础，另外数字信号处理是一门理论性非常强的课程，抽象概念较多，涉及的数学知识也较多，而独立学院学生数学基础、信号与系统分析基础都相对比较薄弱，所以学习本课程有一定的难度，易出现学生学习兴趣低、怕学以及厌学等现象。针对该问题，本文提出了一种基于微课的教学方式，将课程内容制作成微课，分享给学生，以提高学生的学习兴趣。

“微课”采用了精细化管理的教学模式，把“慕课”平台上的“高级软件工程”课程加以优化及改编之后得到“微课”平台，加州大学伯克利分校的 Armando Fox 教授将“微课”平台引入了校园，并取得了意想不到的好效果。2013 年，哈佛大学也开始尝试“微课”平台的教学，效果很好，且学术委员会主席Robert Lue 教授对“微课”给出了极高的评价，认为相对于“慕课”教学方式而言，“微课”体现了“几乎不可避免的进化”。科罗拉多大学全球分校最初采用“慕课”大班教学，经过调研，对学生的需求进行分析，最终也将所有课程都改成采用“微课”的在线授课模式。

从人的认知规律出发，人在学习的时候前十分钟注意力最集中，所以“微课”长度一般设置在5-15分钟的。“微课”在阐述问题的方式上除了教师讲解外，还结合动画、视频等各种信息技术，使得对知识的讲解过程更为直观。且由于一个微课视频时间较短，可以突破传统课堂对时间和地点的限制，学生可以随时随地学习，也可以根据个人学习情况调整学习进度，使得学生真正成为学习的主角。数字信号处理教学中出现的问题

近几年来，数字信号处理课程多媒体与黑板板书相结合的教学方式，再结合Matlab软件实验，教学效果有一定的改善，但仍存在不足之处，主要体现为一下几个方面：

1.1 教学内容方面

公式推导比较多，概念比较抽象，内容枯燥，学生无法将所学的内容与工程实践相结合，再者独立学院学生基础较差，所以易出现“太难，听不懂”等不学、厌学现象。

1.2 教学方法方面

传统教学多采用“灌输式”的方法进行授课，教师作为课堂主体，单向地向学生灌输知识，教学的方式过于单调，不足以吸引学生的注意力，和学习兴趣的提高。

1.3 实践教学方面

实验目前仅局限于理论仿真，而没有硬件的实践操作，学生很难把所学的理论知识和理论仿真同实际应用联系起来，而独立学院学生更应该把学习的重心往实践方向倾斜。数字信号处理教学改革措施

针对数字信号处理教学过程中出现的3个问题，我们可以从以下3个方面进行改进。

2.1 教学内容的改革

数字信号处理课程的知识结构可分为变换域和数字滤波器两部分。对于变换域部分的知识，如离散时间信号的时域频域分析部分与信号与系统分析Z变换部分是重复的，所以在实际教学中，可以对这部分内容进行精简和提炼，并且删除一些复杂的公式推导。关于时域与变换域之间的关系、波形变化以及滤波器的原理等抽象的内容，设计原理的讲解可以精简，而应侧重如何通过Matlab实现具体滤波器的设计，另外，在课堂中增加数字信号应用实例的MALTAB程序讲解，演示仿真波形，通过多波形的分析，使学生理解时域信号好频域信号的特性，从而进一步加深学生对理论知识的理解。

2.2 教学方法和教学手段的改革

2.2.1 引入微课形式的教学方法

在传统“灌输式”、PPT讲解的基础上，借助视频学习、慕课、微课丰富教学手段，将课程重要内容制作成微课，在课堂上播放，吸引学生注意力，提高学生学习兴趣；也可放到网络平台上，供学生学习，充分利用学生的碎片化时间提供在线学习的机会。甚至能够翻转课堂，学生可在课前通过微课自主完成课堂学习，而在课堂上对重点难点进行讨论，使学生在讨论的过程中深化对知识点的理解，同时可以活跃课堂气氛，让学生成为课堂的主体。

2.2.2 双向互动式研究型学习

针对课程部分重点难点教学内容，提取出合适的讨论课题，进行定期或不定期的专题讨论。首先，通过课堂或是其他网络平台向学生发布题目，让学生提前做好讨论或者报告准备，采用多媒体上讲台进行专题讨论或报告。学生通过对课题的准备及讨论过程可对该课程内容有更深入的理解，且能培养学生的辩证思考的能力，另外教师对课题讨论的主要参与人做相应的奖励，比如加平时分，来诱导学生自主学习。

2.3 实践教学的改革

部分独立学院目前已经开展了数字信号处理的软件实验课程，即基于Matlab软件，通过编程实现对信号的仿真，可一定程度上加深学生对信号时域频域的形象理解，可提高学生的编程能力，但毕竟Matlab软件实验只是做信号建模及仿真，所以很多学生仍然存在“学了这门课，却不知道怎么在实际中应用”的疑问，因此，增加DSP硬件实验是必要的。然而增加DSP硬件实验需要校方新建DSP实验室，该问题不是朝夕之间能解决的问题。所以针对此问题，可先进行探索性实验，即通过收集整理其他院校的DSP硬件实现演示视频，使学生对理论的仿真到实际的硬件实现有一个形象全面的了解，比如较抽象的知识点像卷积运算、FIR/IIR数字滤波器以及快速傅里叶变换FFT的DSP实现，通过实验课使学生把抽象的知识与实际应用相结合。总结

实践证明，基于微课及双向互动式的教学方法能充分调动学生的学习积极性，提高学生的自学能力。另外，引入DSP实践教学，使学生更能将理论知识与工程实践相结合。

【参考文献】

[1]曾伟梁，刘颖，鲍曼，徐宏艳，赵沛丰.用 MATLAB 实践数字信号处理课程的辅助教学[J].黑龙江生态工程职业学院学报，2014（1）.[2]高静，王凤文，舒冬梅.《数字信号处理》课程教学实践与探索[J].科技创新导报，2011（4）：153-154.[3]曾伟梁，刘颖，鲍曼，徐宏艳，赵沛丰.用MATLAB 实践数字信号处理课程的辅助教学[J].黑龙江生态工程职业学院学报，2014（1）.[4]刘伯红，阎英，方义秋.高校计算机实践教学质量保障体系改革探索与实践[J].实验室研究与探索，2012，31（12）：121-123，150.[责任编辑：杨玉洁]

第四篇：数字信号处理课程设计..

课程设计报告

课程名称： 数字信号处理 课题名称： 语音信号的处理与滤波

姓 名： 学 号： 院 系： 专业班级： 指导教师： 完成日期： 2013年7月2日

目录

第1部分 课程设计报告………………………………………3 一．设计目的……………………………………………3 二．设计内容……………………………………………3 三．设计原理……………………………………………3 四．具体实现……………………………………………5 1.录制一段声音…………………………………5 2.巴特沃斯滤波器的设计………………………8 3.将声音信号送入滤波器滤波…………………13 4.语音信号的回放………………………………19 5.男女语音信号的频谱分析……………………19 6.噪声的叠加和滤除……………………………22 五． 结果分析……………………………………………27 第2部分 课程设计总结………………………………28 一． 参考文献……………………………………………28

第1部分 课程设计报告

一．设计目的

综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

二．设计内容

录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除；

三．设计原理

1.在Matlab软件平台下，利用函数wavrecord()，wavwrite(),wavread(),wavplay()对语音信号进行录制，存储，读取，回放。

2.用y=fft(x)对采集的信号做快速傅立叶变换，并用[h1,w]=freqz(h)进行DTFT变换。

3.掌握FIR DF线性相位的概念，即线性相位对h(n)、H()及零点的约束，了解四种FIR DF的频响特点。

4.在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波。

5.抽样定理

连续信号经理想抽样后时域、频域发生的变化（理想抽样信号与连续信号频谱之间的关系）

理想抽样信号能否代表原始信号、如何不失真地还原信号即由离散信号恢复连续信号的条件（抽样定理）

理想采样过程描述： 时域描述：

ˆa(t)xa(t)T(t)xa(t)(tnT)xa(nT)(tnT)xnnT(t)频域描述：利用傅氏变换的性质，时域相乘频域卷积，若

n(tnT)ˆa(t)Xa(j)xXa(j)xa(t)T(j)T(t)

则有

ˆ(j)1X(j)(j)XaaT2121ˆXa(j)Xa(jjk)Xa(jjks)TkTTkˆ(j)与X(j)的关系：理想抽样信号的频谱是连续信号频谱的Xaa

周期延拓,重复周期为s(采样角频率)。如果：

X(j)Xa(j)a0s/2s/2即连续信号是带限的，且信号最高频率不超过抽样频率的二分之一，则可不失真恢复。

奈奎斯特采样定理：要使实信号采样后能够不失真还原，采样频率必须大于信号最高频率的两倍：s2h 或 fs2fh

四．具体实现

1.录制一段声音

1.1录制并分析

在MATLAB中用wavrecord、wavread、wavplay、wavwrite对声音进行录制、读取、回放、存储。

程序如下：

Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号 fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令

wavplay(x,Fs);%按任意键播放语音信号

wavwrite(x,Fs,'C:UsersacerDesktop数字信号sound.wav');%存储语音信号

N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(2);subplot(2,1,1);plot(x);%录制信号的时域波形 title('原始信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格

y0=fft(x);%快速傅立叶变换 figure(2);subplot(2,1,2);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 title('原始信号的频谱图');%加标题

grid;%加网格

图1.1 原始信号的时域与频谱图

1.2滤除无效点

针对实际发出声音落后录制动作半拍的现象，如何拔除对无效点的采样的问题: 出现这种现象的原因主要是录音开始时，人的反应慢了半拍，导致出现了一些无效点，而后而出现的无效的点，主要是已经没有声音的动作，先读取声音出来，将原始语音信号时域波形图画出来，根据己得到的信号，可以在第二次读取声音的后面设定采样点，取好有效点，画出滤除无效点后的语音信号时域波形图，对比可以看出。这样就可以解决这个问题。

x=wavread('C:UsersacerDesktop数字信号sound.wav', 7

[4000,24000]);%从4000点截取到24000结束 plot(x);%画出截取后的时域图形 title('截取后的声音时域图形');%标题 xlabel('频率');ylabel('振幅');grid;%画网格

图1.2 去除无效点

2.巴特沃斯滤波器的设计

2.1设计巴特沃思低通滤波器

MATLAB程序如下。滤波器图如图3.3所示。

%低通滤波

fp=1000;fs=1200;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯低通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');gridon;9

图2.1 巴特沃思低通滤波器

2.2设计巴特沃思高通滤波器

MATLAB程序如下。滤波器图如图3.5所示。%高通滤波

fp=4800;fs=5000;Fs=22050;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);10

[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'high','s');[Bz,Az]=bilinear(B,A,Fs1);figure(1);[h,w]=freqz(Bz,Az,512,Fs1*22050);plot(w,abs(h));title('巴特沃斯高通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');grid on;

图2.2巴特沃思高通滤波器

2.3设计巴特沃思带通滤波器

MATLAB程序如下。滤波器图如图3.7所示。%带通滤波

fp=[1200,3000];fs=[1000,3200];Fs=8000;rp=1;rs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;T=1;Fs1=1;wap=2*tan(wp/2);was=2*tan(ws/2);[N,wc]=buttord(wap,was,rp,rs,'s');[B,A]=butter(N,wc,'s');[Bz,Az]=bilinear(B,A,Fs1);figure(4);[h,w]=freqz(Bz,Az,512,Fs1*1000);plot(w,abs(h));title('巴特沃斯带通滤波器');xlabel('频率（HZ）');ylabel('耗损（dB）');grid on;12

图2.3巴特沃思带通滤波器

3．将声音信号送入滤波器滤波

x=wavread('C:UsersacerDesktop数字信号sound.wav');%播放原始信号

wavplay(x,fs);%播放原始信号 N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(4);subplot(4,2,1);plot(x);%录制信号的时域波形

title('原始信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格

y0=fft(x);%快速傅立叶变换 subplot(4,2,3);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 title('原始信号的频谱图');%加标题 grid;%加网格

3．1低通滤波器滤波 fs=8000;beta=10.056;wc=2*pi*1000/fs;ws=2*pi*1200/fs;width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi /width);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);

ys=fftfilt(h,x);%信号送入滤波器滤波，ys为输出 fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 figure(4);subplot(4,2,2);%在四行两列的第二个窗口显示图形 plot(ys);%信号的时域波形

title('低通滤波后信号的时域波形');%加标题 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 grid;%网格

subplot(4,2,4);%在四行两列的第四个窗口显示图形 plot(f, abs(fftwave(n1)));%绘制模值 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义

title('低通滤波器滤波后信号的频谱图');%标题 grid;%加网格

wavplay(ys,8000);%播放滤波后信号

3．2高通滤波器滤波 fs=8000;beta=10.056;ws=2*5000/fs;wc=2*4800/fs;

width=ws-wc;wn=(ws+wc)/2;n=ceil(12.8*pi/width);h=fir1(n,wn/pi, 'high',kaiser(n+2,beta));[h1,w]=freqz(h);ys=fftfilt(h,x);%将信号送入高通滤波器滤波 subplot(4,2,5);%在四行两列的第五个窗口显示图形 plot(ys);%信号的时域波形 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('高通滤波后信号的时域波形');%标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%网格

fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 subplot(4,2,7);%在四行两列的第七个窗口显示图形 plot(f,abs(fftwave(n1)));%绘制模值 axis([0 1 0 50]);xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义

title('高通滤波器滤波后信号的频谱图');%标题 grid;%加网格

wavplay(ys,8000);%播放滤波后信号

3．3带通滤波器 fs=8000;beta=10.056;wc1=2*pi*1000/fs;wc2=2*pi*3200/fs;ws1=2*pi*1200/fs;ws2=2*pi*3000/fs;width=ws1-wc1;wn1=(ws1+wc1)/2;wn2=(ws2+wc2)/2;wn=[wn1 wn2];n=ceil(12.8/width*pi);h=fir1(n,wn/pi,'band',kaiser(n+1,beta));[h1,w]=freqz(h);ys1= fftfilt(h,x);%将信号送入高通滤波器滤波 figure(4);subplot(4,2,6);%在四行两列的第六个窗口显示图形 plot(ys1);%绘制后信号的时域的图形 title('带通滤波后信号的时域波形');%加标题 xlabel('频率w/pi');ylabel('幅值/A');%显示纵坐标表示的意义 grid;%网格

fftwave=fft(ys1);%对滤波后的信号进行快速傅立叶变换 subplot(4,2,8);%在四行两列的第八个窗口显示图形

plot(f, abs(fftwave(n1)));%绘制模值 axis([0 1 0 50]);xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 title('带通滤波器滤波后信号的频谱图');%加标题 grid;%网格

wavplay(ys1,8000);%播放滤波后信号 图形如下：

原始信号的时域波形幅值/A0-1012x 10原始信号的频谱图34幅值/A1低通滤波后信号的时域波形0.50-0.5012频率w/pi3400.51频率w/pi高通滤波后信号的时域波形幅值/A0幅值/A0幅值/Ax 10高通滤波器滤波后信号的频谱图5012频率w/pi34幅值/A0.20-0.2幅值/A2001000x 10低通滤波器滤波后信号的频谱图200100000.51频率w/pi带通滤波后信号的时域波形0.50-0.501234频率w/pix 10带通滤波器滤波后信号的频谱图50幅值 00.5频率w/pi1000.5频率w/pi1

分析：三个滤波器滤波后的声音与原来的声音都发生了变化。其中低

通的滤波后与原来声音没有很大的变化，其它两个都又明显的变化

4.语音信号的回放

sound(xlow,Fs,bits);%在Matlab中，函数sound可以对声音进行回放，其调用格式: sound(xhigh, Fs,bits);%sound(x, Fs, bits);sound(xdaitong, Fs,bits);5.男女语音信号的频谱分析

5.1 录制一段异性的声音进行频谱分析

Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号 fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令 wavplay(x,Fs);%按任意键播放语音信号

wavwrite(x,Fs,'C:UsersacerDesktop数字信号sound2.wav');%存储语音信号

5.2 分析男女声音的频谱

x=wavread(' C:UsersacerDesktop数字信号sound2.wav ');%播放原始信号,解决落后半拍

wavplay(x,fs);%播放原始信号 N=length(x);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;

f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(1);subplot(2,2,1);plot(x);%录制信号的时域波形

title('原始女生信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格

y0=fft(x);%快速傅立叶变换 subplot(2,2,2);plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始女生信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 grid;%加网格

[y,fs,bits]=wavread(' C:UsersacerDesktop数字信号sound.wav ');% 对语音信号进行采样

wavplay(y,fs);%播放原始信号 N=length(y);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 subplot(2,2,3);plot(y);%录制信号的时域波形

title('原始男生信号的时域波形');%加标题 ylabel('幅值/A');%显示纵坐标的表示意义 grid;%加网格

y0=fft(y);%快速傅立叶变换

subplot(2,2,4);%在四行两列的第三个窗口显示图形 plot(f,abs(y0(n1)));%原始信号的频谱图 title('原始男生信号的频谱图');%加标题 xlabel('频率w/pi');%显示横坐标表示的意义 ylabel('幅值 ');%显示纵坐标表示的意义 grid;%加网格

5.3男女声音的频谱图

原始女生信号的时域波形0.50-0.5-1150100原始女生信号的频谱图幅值/A幅值 012345000x 10原始男生信号的时域波形0.50.5频率w/pi原始男生信号的频谱图1300200幅值/A0幅值 012x 1034100-0.5000.5频率w/pi1

图5.3男女声音信号波形与频谱对比

分析：就时域图看，男生的时域图中振幅比女生的高，对于频谱图女生的高频成分比较多

6.噪声的叠加和滤除

6.1录制一段背景噪声

Fs=8000;%抽样频率 time=3;%录音时间 fprintf('按Enter键录音%ds',time);%文字提示 pause;%暂停命令 fprintf('录音中......');x=wavrecord(time*Fs,Fs,'double');%录制语音信号

fprintf('录音结束');%文字提示 fprintf('按Enter键回放录音');pause;%暂停命令 wavplay(x,Fs);%按任意键播放语音信号 wavwrite(x,Fs,'C:UsersacerDesktop数字信号噪音.wav');%存储语音信号

6.2 对噪声进行频谱的分析

[x1,fs,bits]=wavread(' C:UsersacerDesktop数字信号噪音.wav ');%对语音信号进行采样

wavplay(x1,fs);%播放噪声信号 N=length(x1);%返回采样点数 df=fs/N;%采样间隔

n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(5);subplot(3,2,1);plot(x1);%信号的时域波形 title('噪声信号的时域波形');grid;ylabel('幅值/A');y0=fft(x1);%快速傅立叶变换

subplot(3,2,2);plot(f,abs(y0(n1)));%噪声信号的频谱图 ylabel('幅值');title('噪声信号的频谱图');

6.3原始信号与噪音的叠加

fs=8000;[x,fs,bits]=wavread(' C:UsersacerDesktop数字信号sound.wav ');%对录入信号进行采样

[x1,fs,bits]=wavread(' C:UsersacerDesktop数字信号噪音.wav ');%对噪声信号进行采样

yy=x+x1;%将两个声音叠加

6.4叠加信号的频谱分析：

wavplay(yy,fs);%播放叠加后信号 N=length(yy);%返回采样点数 df=fs/N;%采样间隔 n1=1:N/2;f=[(n1-1)*(2*pi/N)]/pi;%频带宽度 figure(5);subplot(3,2,3);plot(yy,'LineWidth',2);%信号的时域波形

title('叠加信号的时域波形');xlabel('时间/t');ylabel('幅值/A');grid;y0=fft(yy);%快速傅立叶变换 subplot(3,2,4);plot(f,abs(y0(n1)));%叠加信号的频谱图 title('叠加信号的频谱图');xlabel('频率w/pi');ylabel('幅值/db');grid;

6.5 设计一个合适的滤波器将噪声滤除 fs=18000;%采样频率 Wp=2*1000/fs;%通带截至频率 Ws=2*2000/fs;%阻带截至频率 Rp=1;%最大衰减 Rs=100;%最小衰减

[N,Wn]=buttord(Wp,Ws,Rp,Rs);%buttord函数（n为阶数，Wn为截至频率）

[num,den]=butter(N,Wn);%butter函数（num为分子系数den为分母系数）

[h,w]=freqz(num,den);%DTFT变换

ys=filter(num,den,yy);%信号送入滤波器滤波，ys为输出 fftwave=fft(ys);%将滤波后的语音信号进行快速傅立叶变换 figure(5);subplot(3,2,5);plot(ys);%信号的时域波形

title('低通滤波后信号的时域波形');%加标题 ylabel('幅值/A');%显示标表示的意义 grid;%网格 subplot(3,2,6);plot(f, abs(fftwave(n1)));%绘制模值 title('低通滤波器滤波后信号的频谱图');%标题 xlabel('频率w/pi');ylabel('幅值/A');%显示标表示的意义 grid;%加网格

wavplay(ys,8000);%播放滤波后信号 grid;图形如下：

噪声信号的时域波形1100噪声信号的频谱图幅值/A0-1幅值0123450000.5叠加信号的频谱图1x 10叠加信号的时域波形10-101时间/t2200幅值/db34幅值/A100000.5频率w/pi1x 10低通滤波后信号的时域波形0.5低通滤波器滤波后信号的频谱图200幅值/A0-0.5幅值/A012x 1034100000.5频率w/pi1

图6.1噪音的叠加与滤除前后频谱对比

7.结果分析

1.录制刚开始时，常会出现实际发出声音落后录制动作半拍，可在[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav')加 窗[x,fs,bits]=wavread('d:matlavworkwomamaaiwo.wav'，[100 10000])，窗的长度可根据需要定义。

2.语音信号通过低通滤波器后，把高频滤除，声音变得比较低沉。当通过高通滤波器后，把低频滤除，声音变得比较就尖锐。通过带通滤波器后，声音比较适中。

3.通过观察男生和女生图像知：时域图的振幅大小与性别无关，只与说话人音量大小有关，音量越大，振幅越大。频率图中，女生高 27

频成分较多。

4.叠加噪声后，噪声与原信号明显区分,但通过低通滤波器后，噪声没有滤除，信号产生失真。原因可能为噪声与信号频率相近无法滤除。

第2部分 课程设计总结

通过本次课程设计，使我们对数字信号处理相关知识有了更深刻的理解，尤其是对各种滤波器的设计。在设计的过程中遇到了很多问题，刚刚开始时曾天真的认为只要把以前的程序改了参数就可以用了，可是问题没有我想象中的那么简单，单纯的搬程序是不能解决问题的。通过查阅资料和请教同学收获了很多以前不懂的理论知识。再利用所学的操作，发现所写的程序还是没有能够运行，通过不断地调试，运行，最终得出了需要的结果。整个过程中学到了很多新的知识，特别是对Matlab的使用终于有些了解。在以后的学习中还需要深入了解这方面的内容。在这次的课程设计中让我体会最深的是：知识来不得半点的马虎。也认识到自己的不足，以后要进一步学习。

八．参考文献

[1]数字信号处理教程(第三版)程佩青 清华大学出版社 [2]MATLAB信号处理 刘波 文忠 电子工业出版社 [3]MATLAB7.1及其在信号处理中的应用 王宏 清华大学出版社

[4]MATLAB基础与编程入门 张威 西安电子科技大学出版社

[5] 数字信号处理及其MATLAB实验 赵红怡 张常 化学工业出版社

[6]MATLAB信号处理详解 陈亚勇等 人民邮电出版社 [7] 数字信号处理

钱同惠 机械工业出版社 29

第五篇：数字信号处理实验报告

JIANGSU

UNIVERSITY OF TECHNOLOGY

数字信号处理实验报告

学院名称： 电气信息工程学院

专 业：

班 级： 姓 名： 学 号： 指导老师： 张维玺（教授）

2013年12月20日

实验一 离散时间信号的产生

一、实验目的

数字信号处理系统中的信号都是以离散时间形态存在的，所以对离散时间信号的研究是数字信号的基本所在。而要研究离散时间信号，首先需要产生出各种离散时间信号。使用MATLAB软件可以很方便地产生各种常见的离散时间信号，而且它还具有强大绘图功能，便于用户直观地处理输出结果。

通过本实验，学生将学习如何用MATLAB产生一些常见的离散时间信号，实现信号的卷积运算，并通过MATLAB中的绘图工具对产生的信号进行观察，加深对常用离散信号和信号卷积和运算的理解。

二、实验原理

离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用x（n）来表示，自变量必须是整数。常见的离散信号如下：（1）单位冲激序列δ（n）

如果δ（n）在时间轴上延迟了k个单位，得到δ（n-k），即长度为N的单位冲激序列δ（n）可以通过下面的MATLAB命令获得。

n=-（N-1）：N-1 x=[zeros(1,N-1)1 zeros(1,N-1)]； stem(n,x)延迟K个采样点的长度为N的单位冲激序列δ（n-k）（k

n=0：N-1 y=[zeros(1,M)1 zeros(1,N-M-1)]； stem（n，y）

（2）单位阶跃序列u（n）

如果u（n）在时间轴上延迟了k个单位，得到u（n-k），即长度为N的单位阶跃序列u（n）可以通过下面的MATLAB命令获得。

n=-（N-1）：N-1 x=[zeros(1,N-1)ones(1,N)]； stem(n,x)延迟的单位阶跃序列可以使用类似于单位冲激序列的方法获得。（3）矩形序列

矩形序列有一个重要的参数，就是序列的宽度N。矩形序列与u（n）之间的关系为矩形序列等= u（n）— u（n-N）。

因此，用MATLAB表示矩形序列可利用上面的单位阶跃序列组合而成。（4）正弦序列x（n）

这里，正弦序列的参数都是实数。与连续的正弦信号不同，正弦序列的自变量n必须为整数。可以证明，只有当2π/w为有理数时，正弦序列具有周期性。

长度为N的正弦序列x（n）可以通过下面的MATLAB命令获得。n=0：N-1 x=A*cos（2*pi*f*n/Fs+phase）（5）单边实指数序列x（n）

长度为N的实指数序列x（n）可以通过下面的MATLAB命令实现。n=0：N-1 x=a.^n stem(n,x)单边指数序列n的取值范围为n>=0。当|a|>1时，单边指数序列发散；当|a|<1时，单边指数序列收敛。当a>0时，该序列均取正值；当a<0时，序列在正负摆动。

（6）负指数序列x（n）

当a=0时，得到虚指数序列x（n）。

与连续负指数信号一样，我们将负指数序列实部和虚部的波形分开讨论，得到如下结论：

1）当a>0时，负指数序列x（n）的实部和虚部分别是按指数规律增长的正弦振荡序列；

2)当a<0时，负指数序列x（n）的实部和虚部分别是按指数规律衰减的正弦振荡序列；

3)当a=0时，负指数序列x（n）即为虚指数序列，其实部和虚部分别是等幅的正弦振荡序列；

长度为N的实指数序列x（n）可以通过下面的MATLAB命令实现。n=0：N-1 x=exp（（a.+j*w）*n）stem(n,real(x))或

stem(n,imag(x))

三、实验内容及分析

1n01、编制程序产生单位冲激序列n“并绘出其图及n”学号后两位0n0形。程序:(1)N=4;

n=-(N-1):N-1;

x=[zeros(1,N-1)1 zeros(1,N-1)];stem(n,x);

title('单位冲激序列');

grid on;

(2)N=6;

M=1;%学号01 n=-(N-1):N-1;

y=[zeros(1,N-M+1)1 zeros(1,N-M-1)];stem(n,y);

title('单位冲激序列');grid on;

分析:在上图的基础上向右平移了1个单位。

1n02、编制程序产生单位阶跃序列un、un“学号后两位”及

0n0unun“学号后两位”，并绘出其图形。程序: 4

(1)N=5;

n=-(N-1):N-1;

x=[zeros(1,N-1)ones(1,N)];stem(n,x);

title('单位阶跃序列');grid on;

(2)N=6;

M=1;%学号01 n=-(N-1):N-1;

x=[zeros(1,N-M+1)ones(1,N-M)];stem(n,x);

title('单位阶跃序列');grid on;

分析:在上图的基础上平移了1个单位.(3)N=6;

M=1;%学号01 n=-(N-1):N-1;

x=[zeros(1,N-1)ones(1,N)];y=[zeros(1,N-M+1)ones(1,N-M)];z=x-y;stem(n,z);

title('单位阶跃序列');grid on;

2

3、编制程序产生正弦序列xncos2n、xncosn及

学号后两位xnsin2n并绘出其图形。

程序:(1)N=5;

A=1;

w=2*pi;phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信号');grid on;

分析:该序列具有周期性,且输出为余弦信号.(2)N=5;

A=1;

w=2*pi/1;%学号01 phi=0;n=0:0.05:N-1;x=A*cos(w*n+phi);stem(n,x);title('余弦信号');grid on;

;

分析:该序列具有周期性,且输出为余弦信号.(3)N=5;

A=1;

w=2*pi;phi=0;

n=0:0.05:N-1;x=A*sin(w*n+phi);stem(n,x);title('正弦信号');grid on;

分析:该序列具有周期性,且输出为正弦信号.4、编制程序产生复正弦序列xne(2j学号后两位）n，并绘出其图形。N=3;

n=0:0.2:N-1;

w=1;%学号01 x=exp((2+j*w)*n);subplot(2,1,1)

stem(n,real(x)),title('实部');grid on;subplot(2,1,2)

stem(n,imag(x)),title('虚部');grid on;

5、编制程序产生指数序列xnan，并绘出其图形。其中a=学号后两位、a=1/“学号后两位”。

(1)N=10;

n=0:N-1;

a=1;%学号01 x=a.^n;stem(n,x);title('指数序列');grid on;

(2)N=10;

n=0:N-1;

a=1;%学号01 x=a.^(-n);stem(n,x);title('指数序列');grid on;

实验三 离散时间信号的频域分析

一、实验目的

信号的频域分析是信号处理中一种有效的工具。在离散信号的频域分析中，通常将信号表示成单位采样序列的线性组合，而在频域中，将信号表示成复变量或的线性组合。通过这样的表示，可以将时域的离散序列映射到频域以便于进一步的处理。

在本实验中，将学习利用MATLAB计算离散时间信号的DTFT和DFT，并加深对其相互关系的理解。

二、实验原理

（１）DTFT和DFT的定义及其相互关系。

（２）使用到的MATLAB命令有基于DTFT离散时间信号分析函数以及求解序列的DFT函数。

三、实验内容及分析

（１）编程计算并画出下面DTFT的实部、虚部、幅度和相位谱。

X(e)jw0.05180.1553e11.2828ex(n)cosjwjw0.1553ej2w1.0388ej2w0.0518ej3w0.3418ej3w

（２）计算32点序列

5n16，0≦n≦31的32点和64点DFT，分别绘出幅度谱图形，并绘出该序列的DTFT图形。

３－１

clear;

x=[0.0518,-0.1553,0.1553,0.0518];y=[1,1.2828,1.0388,0.3418];w=[0:500]*pi/500 H=freqz(x,y,w);

magX=abs(H);angX=angle(H);realX=real(H);imagX=imag(H);subplot(221);plot(w/pi,magX);grid;

xlabel('frequency in pi unit');ylabel('magnitude');title('幅度 part');axis([0 0.9 0 1.1]);

subplot(223);plot(w/pi,angX);grid;

xlabel('frequency in pi unit');ylabel('radians');title('相位 part');axis([0 1-3.2 3.2]);

subplot(222);plot(w/pi,realX);grid;

xlabel('frequency in pi unit');ylabel('real part');title('实部 part');axis([0 1-1 1]);

subplot(224);plot(w/pi,imagX);grid;

xlabel('frequency in pi unit');ylabel('imaginary');title('虚部 part');axis([0 1-1 1.1]);

３－２

N=32;n=0:N-1;

xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('32点');figure N=64;n=0:N-1;

xn=cos(5*pi*n/16);k=0:1:N-1;Xk=fft(xn,N);subplot(2,1,1);stem(n,xn);subplot(2,1,2);stem(k,abs(Xk));title('64点');

（１）

（２）

实验四 离散时间LTI系统的Z域分析

一、实验目的

本实验通过使用MATLAB函数对离散时间系统的一些特性进行仿真分析，以加深对离散时间系统的零极点、稳定性，频率响应等概念的理解。学会运用MATLAB分析离散时间系统的系统函数的零极点；学会运用MATLAB分析系统函数的零极点分布与其时域特性的关系；学会运用MATLAB进行离散时间系统的频率特性分析。

二、实验原理

离散时间系统的系统函数定义为系统零状态响应的Z变化与激励的Z变化之比。

在MATLAB中系统函数的零极点可通过函数roots得到，也可借助函数tf2zp得到，tf2zp的语句格式为

[Z,P,K]=tf2zp(B,A)其中，B与A分别表示H（z）的分子与分母多项式的系数向量。它的作用是将H（z）的有理分式表示式转换为零极点增益形式。

若要获得系统函数H（z）的零极点分布图，可直接应用zplane函数，其语句格式为

Zplane（B，A）

其中，B与A分别表示H（z）的分子和分母多项式的系数向量。它的作用是在z平面上画出单位圆、零点与极点。

离散系统中z变化建立了时域函数h（n）与z域函数H（z）之间的对应关系。因此，z变化的函数H（z）从形式可以反映h（n）的部分内在性质。可根据系统的传递函数H（z）求单位冲激响应h（n）的函数impz、filter等。

利用系统的频率响应，可以分析系统对各种频率成分的响应特性，并推出系统的特性（高通、低通、带通、带阻等）。

MATLAB提供了求离散时间系统频响特性的函数freqz，调用freqz的格式主要有两种。一种形式为

[H,w]= reqz(B,A,N)其中，B与A分别表示H（z）分子和分母多项式的系数向量；N为正整数，默认值为512；返回值w包含[0，π]范围内的N个频率等分点；返回值H则是离散时间系统频率响应在0~π范围内N个频率处的值。另一种形式为

[H，w]= freqz（B，A，N，‘whole’）

与第一种方式不同之处在于角频率的范围由[0，π]扩展到[0，2π]。

三、实验内容与结果分析

已知LTI离散时间系统，要求由键盘实现系统参数输入，并绘出幅频和相频响应曲线和零极点分布图，进而分析系统的滤波特性和稳定性。

（一）程序

b=[0.0528,0.797,0.1295,0.1295,0.797,0.0528];

a=[1,-1.8107,2.4947,-1.8801,0.9537,-0.2336];w=[0:20:500]*pi/500;

x1=0.0528+0.797*exp(-1*j*w)+0.1295*exp(-2*j*w)+0.1295*exp(-3*j*w)+0.797*exp(-4*j*w)+0.0528*exp(-5*j*w);

x2=1-1.8107*exp(-1*j*w)+2.4947*exp(-2*j*w)+1.8801*exp(-3*j*w)+0.9537*exp(-4*j*w)+0.2336*exp(-5*j*w);x22=x2+(x2==0)*eps;x=x1./x22;magx=abs(x);

angx=angle(x).*180/pi;

subplot(2,2,3);zplane(b,a);title('零极点图');subplot(2,2,2);stem(w/pi,magx);title('幅度部分');ylabel('振幅');subplot(2,2,4);stem(w/pi,angx);

xlabel('以pi为单位的频率');title('相位部分');ylabel('相位');

（二）波形图

图4-1 幅频、相频响应曲线、零极点分布图

实验六 IIR数字滤波器的设计

一、实验目的

从理论上讲，任何的线性是不变（LTI）离散时间系统都可以看做一个数字滤波器，因此设计数字滤波器实际就是设计离散时间系统。数字滤波器你包括IIR（无限冲激响应）和FIR（有限冲激响应）型，在设计时通常采用不同的方法。

本实验通过使用MATLAB函数对数字滤波器进行设计和和实现，要求掌握IIR数字巴特沃斯滤波器、数字切比雪夫滤波器的设计原理、设计方法和设计步骤；能根据给定的滤波器指标进行滤波器设计；同时也加深学生对数字滤波器的常用指标和设计过程的理解。

二、实验原理

在IIR滤波器的设计中，常用的方法是：先根据设计要求寻找一个合适的模拟原型滤波器，然后根据一定的准则将此模拟原型滤波器转换为数字滤波器。

IIR滤波器的阶数就等于所选的模拟原型滤波器的阶数，所以其阶数确定主要是在模拟原型滤波器中进行的。

IIR数字滤波器的设计方法如下：（1）冲激响应不变法。（2）双线性变化法。

一般来说，在要求时域冲激响应能模仿模拟滤波器的场合，一般使用冲激响应不变法。冲激响应不变法一个重要特点是频率坐标的变化是线性的，因此如果模拟滤波器的频率响应带限于折叠频率的话，则通过变换后滤波器的频率响应可不失真地反映原响应与频率的关系。

与冲激响应不变法比较，双线性变化的主要优点是靠频率的非线性关系得到s平面与z平面的单值一一对应关系，整个值对应于单位圆一周。所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。

MATLAB提供了一组标准的数字滤波器设计函数，大大简化了滤波器的设计工程。

（1）butter。

（2）cheby1、cheby2。

三、实验内容及分析

利用MATLAB编程方法或利用MATLAB中fdatool工具设计不同功能的IIR数字滤波器。

1、基于chebyshev I型模拟滤波器原型使用冲激不变转换方法设计数字滤波器，要求参数为通带截止频率p0.4；通带最大衰减Ap1dB；阻带截止频率s0.4；阻带最小衰减As35dB。

程序:

wp=0.2*pi;

%通带边界频率

ws=0.4*pi;

%阻带截止频率 rp=1;

%通带最大衰减 rs=35;

%阻带最小衰减

Fs=1000;

%¼ÙÉè³éÑùÂö³å1000hz

[N,Wn]=cheb1ord(wp,ws,rp,rs,'s');

[Z,P,K]=cheby1(N,rp,Wn,'s');[H,W]=zp2tf(Z,P,K);

figure(1);freqs(H,W);[P,Q]=freqs(H,W);figure(2);plot(Q*Fs/(2*pi),abs(P));grid on;

xlabel('频率/Hz');ylabel('幅度');

2、基于Butterworth型模拟滤波器原型使用双线性变换方法设计数字滤波器的，要求参数为截止频率p0.4；通带最大衰减Ap1dB；阻带截止频率s0.25；阻带最小衰减AS40dB。程序: wp=0.4*pi;ws=0.25*pi;rp=1;rs=40;fs=500;ts=1/fs;wp1=wp*ts;ws1=ws*ts;

wp2=2*fs*tan(wp1/2);ws2=2*fs*tan(ws1/2);

[N,Wn]=buttord(wp2,ws2,rp,rs,'s');[Z,P,K]=buttap(N);[Bap,Aap]=zp2tf(Z,P,K);[b,a]=lp2lp(Bap,Aap,Wn);[bz,az]=bilinear(b,a,fs);[H,W]=freqz(bz,az);subplot(2,1,1);plot(W/pi,abs(H));grid on;xlabel('频率')ylabel('幅度')subplot(2,1,2);

plot(W/pi,20*log10(abs(H)));grid on;xlabel('频率');ylabel('幅度(dB)');

实验七 FIR数字滤波器的设计

一、实验目的

掌握用窗函数设计FIR数字滤波的原理及其设计步骤；熟悉线性相位数字滤波器的特性。学习编写数字滤波器的设计程序的方法，并能进行正确编程；根据给定的滤波器指标，给出设计步骤。

二、实验原理

如果系统的冲激响应h（n）为已知，则系统的输入输出关系为

y（n）=x（n）*h（n）

对于低通滤波器，只要设计出低通滤波器的冲激响应函数，就可以由式得到系统的输出了。

但是将h（n）作为滤波器的脉冲响应有两个问题：一是它是无限长的；二是它是非因果的。对此，采取两项措施：一是将h（n）截短；二是将其右移。

设计时，要根据阻带的最小衰减和过渡带宽度来选择恰当的窗函数类型和窗口长度N。常用的窗函数有矩形窗、海明窗和布莱克曼窗等。

窗函数设计FIR滤波器步骤如下：

（1）给定理想频率响应的幅频特性和相频特性；

（2）求理想单位脉冲响应，在实际计算中，可对理想频率响应采样。（3）根据过渡带宽度和阻带最小衰减，确定窗函数类型和窗口长度N；（4）求FIR滤波器单位脉冲响应；

（5）分析幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。

三、实验内容及分析

1、分别用海明窗和布莱克曼窗设计一个48阶的FIR带通滤波器，通带为Wn0.450.55。程序1:海明窗设计

N=48;

Window=hamming(N+1);w1=0.45;w2=0.55;ws=[w1,w2];

b=fir1(N,ws/pi,Window);freqz(b,1,512);title('海明窗');grid on;

程序2:莱克曼窗设计

N=48;

Window=blackman(N+1);w1=0.45;w2=0.55;ws=[w1,w2];

b=fir1(N,ws/pi,Window);freqz(b,1,512);title('布莱克曼窗');grid on;

2、用矩形窗设计一个线性相位高通滤波器。其中Hejwej00.3

00.3程序: N=9;

alpha=(N-1)/2;Wc=0.7*pi;n=(0:8);i=n-alpha;i=i+(i==0)*eps;

h=(-1).^n.*sin((i).*Wc)./((i).*pi);%矩形窗函数设计的系统脉冲响应 w=(0:1:500)*2*pi/500;

H=h*exp(-j*n'*w);%矩形窗函数设计的频响 magH=abs(H);% 矩形窗函数设计的振幅 subplot(211);stem(n,h);

axis([0,8,-0.4,0.4]);title('矩形窗设计h(n)');line([0,10],[0,0]);xlabel('n');ylabel('h');subplot(212);plot(w/pi,magH);

xlabel('以pi为单位的频率');ylabel('H振幅');axis([0,2,0,1.7]);title('矩形窗设计振幅谱');

实验心得体会：

这次实验使我进一步加深了对MATLAB软件的使用。从上次的信号系统实验的初步使用到这一次的深入了解，有了更深刻的认识。对这种语言环境也有了新的了解。

在实验的过程中，我对数字滤波器的整个过程有了很好的理解和掌握。IIR数字滤波器的设计让我知道了巴特沃思滤波器和切比雪夫滤波器的频率特性，还有双线性变换及脉冲响应不变法设计的滤波器的频率特性。做这两个实验的时候程序有点困难，但经过细心的改写图形最终出来了。FIR数字滤波器的设计出来的是两种窗的图形，通过两种窗的比较，我了解了他们各自的特点，幅频和相频特性。

最后，感谢张老师对我的谆谆教导！