基于Matlab的数字信号处理课程设计2.doc

沈阳理工大学信息学院基于Matlab的数字信号处理课程设计 基于Matlab的数字信号处理课程设计一、课程设计目的：1熟悉Matlab运行环境，熟练使用Matlab语言进行编程，进行数字信号处理。2全面复习数字信号课程所学理论知识，巩固所学知识重点和难点，将理论与实践很好地结合起来。3提高综合运用所学知识独立分析和解决问题的能力；二、 MATLAB简介MATLAB是功能强大的科学及计算软件，它不但具有以矩阵计算为基础的强大数学计算和分析功能，而且还具有丰富的可视化图形表现功能和方便的程序设计能力。MATLAB的应用领域极为广泛，除了数学计算和分析外，还被广泛地应用于自动控制、系统仿真、数字信号处理、图形图像分析、数理统计、人工智能、虚拟现实技术、通信工程、金融系统等领域，因此，MATLAB是面向21世纪的计算机程序设计及科学计算语言。三、MATLAB的主要组成部分MATLAB系统包括5个主要部分：（1）开发环境MATLAB开发环境由一组工具和组件组成，这些工具是图形化的用户界面，包括MATLAB桌面和命令窗口、命令历史窗口、帮助信息浏览器、文件和搜索路径浏览器。（2）MATLAB数学函数库MATLAB集成了丰富的数学函数库，其强大的计算能力覆盖了从基本函数（如求和、正弦、余弦和复数运算等）到高级函数（如矩阵求逆、矩阵特征值、贝塞尔函数和快速傅立叶变换等）的范围。（3）MATLAB语言MATLAB语言是一种以矩阵运算为基础的高级语言，包括控制流的描述、函数、数据结构、输入输出及面向对象的编程环境，既可以编制快速使用小程序，也可以编制大型复杂的应用程序。（4）图形功能MATLAB提供了功能强大的图形系统，既可以用高级命令完成二维和三维数据的可视化、图像处理、动画和图形表达等功能，也可以通过使用图形句柄完成复杂的图形功能，实现对所有图形对象的操作。（5）应用程序接口（API）MATLAB还提供了应用程序接口库函数，允许用户使用C或FORTRAN语言编写程序与MATLAB连接，功能包括与MATLAB的动态连接、调用MATLAB作为运算引擎、读写MAT文件等。MATLAB的广泛应用促进了其本身功能的迅速发展，以特定领域的应用为主要目的的应用程序工具箱的数量和内容在不断增多。四、利用Matlab熟悉画图（设计1）1、设计内容：能够利用Matlab熟悉地画图，内容包括：X、Y坐标轴上的label，每幅图上的title，绘画多条曲线时的legend，对图形进行适当的标注等。 （1）在一副图上画出多幅小图；（2）画出一组二维图形； （3）画出一组三维图形； （4）画出复数的实部与虚部。1. 实验结果：（1） 、在一副图上画出多幅小图；（2） 、画出一组二维图形；（3） 、画出一组三维图形；（4） 、画出复数的实部与虚部；（5）、完成对一个源程序进行详细注释。（见程序）% 能够利用Matlab熟悉地画图，内容包括：X、Y坐标轴上的label，% 每幅图上的title，绘画多条曲线时的legend，对图形进行适当的标注等。% （1）在一副图上画出多幅小图；（2）画出一组二维图形；% （3）画出一组三维图形； （4）画出复数的实部与虚部。% （5）完成对一个源程序进行详细注释。%clear;clc;x=linspace(-pi,pi,200); % 定义x轴坐标y1=sin(x); % 函数1y2=cos(x); % 函数2y3=exp(x); % 函数3y4=sin(x).*cos(x); % 函数4% 在一幅图上绘制多幅小图 %figure(1)subplot(221), % 图一plot(x,y1); % 绘制y1grid on; % 添加网格% axis(min(x) max(x) min(y1) max(y1); % 指定显示范围subplot(222), % 图二hist(y2); % 绘制y2grid on; % 添加网格% axis(min(x) max(x) min(y2) max(y2); % 指定显示范围subplot(223), % 图三stem(x,y3); % 绘制y3grid on; % 添加网格% axis(min(x) max(x) min(y3) max(y3); % 指定显示范围subplot(224), % 图四plot(x,y4,rp); % 绘制y4grid on; % 添加网格% axis(min(x) max(x) min(y4) max(y4); % 指定显示范围% 画出一组二维图形 %figure(2) a=linspace(0,2*pi,200);k=1,2,3;y=sin(a)*k;plot(a,y(:,1),a,y(:,2),a,y(:,3);legend(y= ,num2str(k(1),*sin(x),y= ,num2str(k(2),*sin(x). ,y= ,num2str(k(3),*sin(x);xlabel(x);ylabel(y);title(多条曲线绘制);grid on;% 画出一组三维图形 %figure(3) kk=30;x=linspace(-5,5,kk);y=linspace(-5,5,kk);for i=1:kk for j=1:kk z(i,j)=sin(x(i)+cos(y(j); endendsurf(x,y,z);xlabel(x);ylabel(y);zlabel(z);title(z=sin2(x)/(3*y-2);% 复数绘制 %figure(4) clearkk=1:1:10;for i=1:10 y(i)=sin(kk(i)+cos(kk(i)/kk(i)*j;endcompass(y);五、褶积的相关计算与编程（设计2）1. 设计内容：计算普通褶积与循环褶积，分别使用时间域与频率域两种方法进行正、反演计算，指出循环褶积计算时所存在的边界效应现象；编写一个做相关分析的源程序。2. 实验成果：（1）、构造如下两个向量：a=1,2,3,4; b=9,8,7,6; 计算得到：线性褶积9 26 50 80 65 46 24循环次数时间域频率域循环褶积555 50 50 80 6555 50 50 80 65633 26 50 80 65 4633 26 50 80 65 4679 26 50 80 65 46 249 26 50 80 65 46 2489 26 50 80 65 46 24 09 26 50 80 65 46 24 0从计算结果可以看出，当循环次数小于N1+N2-1时，线性褶积和循环褶积不等，当循环次数大于等于N1+N2-1时，线性褶积和循环褶积相等。（2）、构造两个向量：a=1 2 3; b=1 2 3 5;计算相关得到：频率域： 14 13 16 23时间域（6次）： 14 8 3 5 13 23程序：1、 褶 积 % 计算普通褶积与循环褶积，分别使用时间域与频率域两种方法% 进行正、反演计算，指出循环褶积计算时所存在的边界效应现象；% 编写一个做相关分析的源程序。 % 验证循环褶积和线性褶积在什么条件下相等%clc;clear;a=1,2,3,4; % a矩阵b=9,8,7,6; % b矩阵%kk1=length(a);kk2=length(b);% 循环褶积 时间域%xx=8; % 循环次数 手动更改aa=zeros(xx,1);aa(1:kk1,1)=a; % 扩展abb=zeros(xx,1);bb(1:kk2,1)=b; % 扩展bbb1=zeros(xx,xx); % 循环矩阵for i=1:xx bb1(i:xx,i)=bb(1:xx-i+1); bb1(1:i-1,i)=bb(xx-i+2:xx);endok1=bb1*aa; % 循环褶积结果% 线性褶积 时间域%bb2=zeros(kk1+kk2-1,kk1);for i=1:kk1 bb2(i:i+kk2-1,i)=b(1:kk2);endan=a;ok2=bb2*an; % 线性褶积结果% 循环褶积 频率域%af=fft(aa);bf=fft(bb);ok3=ifft(af.*bf);2、 相 关 % 计算普通褶积与循环褶积，分别使用时间域与频率域两种方法进行% 正、反演计算，指出循环褶积计算时所存在的边界效应现象；% 编写一个做相关分析的源程序。 % 求两个一维数组的相关系数%clear;clc;a=1 2 3; % 数组a 较短时可以手动输入b=1 2 3 5; % 数组b 较短时可以手动输入% a=csvread(*.csv,0,0); % 数组a 较长时可以读取csv文件% b=csvread(*.csv,0,0); % 数组a 较长时可以读取csv文件% 频率域求相关%a2=zeros(1,max(length(a),length(b);b2=zeros(1,max(length(a),length(b);a2(1,1:length(a)=a;b2(1,1:length(b)=b;ffta=fft(a2);fftb=fft(b2);fftb=conj(fftb);result=(ifft(ffta.*fftb) % 频率域求相关结果% 时间域求相关%long=length(a)+length(b)-1;aa=zeros(long,1);bb=zeros(long,1);aa(1:length(a),1)=a;bb(1:length(b),1)=b;temp=zeros(long);for i=1:long temp(1:long-i+1,i)=aa(i:long); temp(long-i+2:long,i)=aa(1:i-1);endresult2=temp*(bb) % 时间域求循环相关六、设计一维滤波处理程序（设计3）1. 实验内容：设计一个一维滤波处理程序（1、分别做低通、高通、带通、带阻等理想滤波器进行处理；2、窗函数）2. 实验结果：（1） 、一维信号理想滤波器（2）、窗函数程序：3、 一 维 理 想 滤 波 器 % 一维信号滤波器设计%function varargout = num4_1lvbo(varargin)% Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct(gui_Name, mfilename, . gui_Singleton, gui_Singleton, . gui_OpeningFcn, num4_1lvbo_OpeningFcn, . gui_OutputFcn, num4_1lvbo_OutputFcn, . gui_LayoutFcn, , . gui_Callback, );if nargin & ischar(varargin1) gui_State.gui_Callback = str2func(varargin1);endif nargout varargout1:nargout = gui_mainfcn(gui_State, varargin:);else gui_mainfcn(gui_State, varargin:);endfunction num4_1lvbo_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject;guidata(hObject, handles);function varargout = num4_1lvbo_OutputFcn(hObject, eventdata, handles) varargout1 = handles.output;% 全局变量说明% yuchu 用于存放信号FFT变换结果% long 用于存放信号长度% ys 用于存放原始信号% 读取原始信号并进行预处理%function bt_read_Callback(hObject, eventdata, handles)FileName,PathName = uigetfile(*.csv) global ysys=csvread(FileName);axes(handles.wave_ys);plot(ys(:,1),ys(:,2); grid on;axis(min(ys(:,1) max(ys(:,1) min(ys(:,2)*3/2 max(ys(:,2)*3/2);global yuchu yuchu=fft(ys(:,2);global long long=length(ys(:,2);% 低通滤波%function bt_dt_Callback(hObject, eventdata, handles)global yuchu long ysf1=str2double(get(handles.input_f1,string);for i=1:long if i=(fix(long/2)+1)-(fix(long/2)+1)*f1+. (fix(long/2)+1) dt(i)=yuchu(i); else dt(i)=0; endendok=ifft(dt);axes(handles.wave_dt);plot(ys(:,1),ok);grid on;% 高通滤波%function bt_gt_Callback(hObject, eventdata, handles)global yuchu long ysf1=str2double(get(handles.input_f1,string);for i=1:long if i=(fix(long/2)+1)-(fix(long/2)+1)*f1+. (fix(long/2)+1) gt(i)=0; else gt(i)=yuchu(i); endendok=ifft(gt);axes(handles.wave_gt);plot(ys(:,1),ok);grid on;% 带阻滤波%function bt_dz_Callback(hObject, eventdata, handles)global yuchu long ysf1=str2double(get(handles.input_f1,string);f2=str2double(get(handles.input_f2,string);for i=1:long if (i=(fix(long/2)+1)*f1 & i=(fix(long/2)+1)*f2) | . (i=(fix(long/2)+1)-(fix(long/2)+1)*f2+(fix(long/2)+1) dz(i)=0; else dz(i)=yuchu(i); endendok=ifft(dz);axes(handles.wave_dz);plot(ys(:,1),ok);grid on;% 带通滤波%function bt_dat_Callback(hObject, eventdata, handles)global yuchu long ysf1=str2double(get(handles.input_f1,string);f2=str2double(get(handles.input_f2,string);for i=1:long if (i=(fix(long/2)+1)*f1 & i=(fix(long/2)+1)*f2) | . (i=(fix(long/2)+1)-(fix(long/2)+1)*f2+(fix(long/2)+1) dz(i)=yuchu(i); else dz(i)=0; endendok=ifft(dz);axes(handles.wave_dat);plot(ys(:,1),ok);grid on;4、 窗 函 数 function varargout = num4_windows(varargin)% 本程序用于对一维信号进行窗函数滤波% 其中，f1、f2、f3、f4依次为% 带通滤波镶边后的四个频率% 为了使滤波效果较好，建议使 f1-f2=f3-f4%gui_Singleton = 1;gui_State = struct(gui_Name, mfilename, . gui_Singleton, gui_Singleton, . gui_OpeningFcn, num4_windows_OpeningFcn, . gui_OutputFcn, num4_windows_OutputFcn, . gui_LayoutFcn, , . gui_Callback, );if nargin & ischar(varargin1) gui_State.gui_Callback = str2func(varargin1);endif nargout varargout1:nargout = gui_mainfcn(gui_State, varargin:);else gui_mainfcn(gui_State, varargin:);endfunction num4_windows_OpeningFcn(hObject, eventdata, handles, varargin)handles.output = hObject;guidata(hObject, handles);function varargout = num4_windows_OutputFcn(hObject, eventdata, handles) varargout1 = handles.output;% 以下为信号读取以及原始信号的绘制%function read_Callback(hObject, eventdata, handles)FileName,PathName = uigetfile(*.csv) % read file global signsign=csvread(FileName);axes(handles.axes_sign)plot(1:length(sign),sign(:,2)grid onfunction plot_Callback(hObject, eventdata, handles)% 采样频率100Hz，Wpl=20Hz，Wph=30Hz，Wsl=15Hz，Wsh=35Hz，% 频率/采样频率*2*pi% Wpl=0.4*pi，Wph=0.6*pi，Wsl=0.2*pi，Wsh=0.8*pi，% Rp=1;As=20clc;T = 1; Fs = 1/T; % Set T=1%按照比值进行带通的四个频率成份的界定 %f1=str2double(get(handles.input_f1,string);f2=str2double(get(handles.input_f2,string);f3=str2double(get(handles.input_f3,string);f4=str2double(get(handles.input_f4,string);%-global signkk=length(sign); % 一维信号长度Wpl=tan(f2/2/Fs*pi);Wph=tan(f3/2/Fs*pi);Wsl=tan(f1/2/Fs*pi);Wsh=tan(f4/2/Fs*pi);Rp = 1;As= 40;% 计算模拟低通原型的参数OmigaP=Wph-Wpl; % 低通滤波器通带截止频率OmigaS=Wsh-Wsl; % 低通滤波器通带截止频率N,Wn=buttord(OmigaP,OmigaS,Rp,As,s);B,A=butter(N,Wn,s);BT,AT=lp2bp(B,A,sqrt(Wph*Wpl),Wph-Wpl);b,a=bilinear(BT,AT,Fs);H=freqz(b,a,kk,whole);axes(handles.axes_ifft) % 输出滤波器plot(1:kk,ifft(abs(H).*fft(sign(:,2),r);grid on;axes(handles.axes_fft) % 输出滤波后的信号plot(abs(H),r);grid on;六、验证褶积性质（设计4）1. 设计内容：验证时间域的循环褶积对应的是频率域的乘积；线性褶积则不然。2. 实验结果：a，b两向量取如下值：a=1,2,3,4; b=1,2,3;计算得到如下结果： 当循环次数依次为4，5，6，7时，在频率域和时间域计算得到的循环褶积均相等；当循环次数小于6（N1+N2-1）时，时间域和频率域计算所得的结果不同，但是当循环次数大于6时，两种方法计算得到的结果相等。因此可以得到如下结论:(1) 两个向量的循环褶积所对应的频谱是两信号各自频谱在频率域的乘积；(2) 当两个向量在频率域计算线性褶积时，要对两个信号进行补零扩展，使两个信号同大小。当扩展出的新信号长度小于临界长度6（N1+N2-1）时，频率域和时间域计算得到的线性褶积不相等；当扩展出的新信号长度大于等于临界长度6（N1+N2-1）时，频率域和时间域计算得到的线性褶积相等。 例如：设length(a)=n1 length(b)=n2 length(新信号)=n3 如果 n3=n1+n2-1 两种方法所得结果相等；