数字信号处理报告—实习报告

1、 数字信号处理报告 学 校： 中国地质大学（武汉） 姓 名： 班 级： 学 号： 指导老师： 目录1-1：离散卷积的运算1-2：离散傅里叶变换及其应用1-3：滤波器的设计1-4：低通滤波器的设计1-5：用设计一信号处理系统1-6:频率抽选1-7：实习心得体会1-1离算卷积的应用1：基本要求设线性时不变（lti）系统的冲激响应为h(n)，输入序列为x(n)1、h(n)=(0.8)n，0n4； x(n)=u(n)-u(n-4)2、h(n)=(0.8)nu(n), x(n)=u(n)-u(n-4)3、h(n)=(0.8)nu(n), x(n)=u(n)求以上三种情况下系统的输出y(n)，显示输入和输

2、出波形。2：实验目的及算法原理实验目的：（1）：掌握离散卷积计算机实现。（2）：进一步对离散信号卷积算法的理解。算法原理：把冲激响应h(n)与输入序列x(n)分别输入到程序中，然后调用离散卷积函数y=conv(x.,h)即可得到所要求的结果。离散卷积定义为当序列为有限长时则3：程序代码x1=1 1 1 1 ;nx1=0:3; h1=1 0.8 0.64 0.83 0.84;nh1=0:4; y1=conv(x1,h1); subplot(3,3,1);stem(nx1,x1);title(序列x1); xlabel(n);ylabel(x1(n); subplot(3,3,2);stem(nh

3、1,h1);title(序列h1); xlabel(n);ylabel(h1(n); subplot(3,3,3);stem(y1);title(序列y1); xlabel(n);ylabel(y1(n); x2=1 1 1 1;nx2=0:3; nh2=0:1:20; h2=(0.8).nh2; y2=conv(x2,h2); subplot(3,3,4);stem(nx2,x2);title(序列x2); xlabel(n);ylabel(x2(n); subplot(3,3,5);stem(nh2,h2);title(序列h2); xlabel(n);ylabel(h2(n); subp

4、lot(3,3,6);stem(y2);title(序列y2); xlabel(n);ylabel(y2(n) ; nx3=0:1:20;x3=1.nx3; nh3=0:1:20; h3=(0.8).nh3; y3=conv(x3,h3); subplot(3,3,7);stem(nx3,x3);title(序列x3); xlabel(n);ylabel(x3(n); subplot(3,3,8);stem(nh3,h3);title(序列h3); xlabel(n);ylabel(h3(n); subplot(3,3,9);stem(y3);title(序列y3); xlabel(n);yl

5、abel(y3(n);4：结果分析 1-2离散傅里叶变换及应用1：基本要求设有离散序列 x(n)=cos(0.48n)+cos(0.52n)分析下列三种情况下的幅频特性。（1） 采集数据长度n=16，分析16点的频谱，并画出幅频特性。采集数据长度n=16，并补零到64点，分析其频谱，并画出幅频特性。（2） 采集数据长度n=64，分析64点的频谱，并画出幅频特性。观察三幅不同的幅频特性图，分析和比较它们的特点及形成原因。2：实习目的及相关算法实习目的：（1）了解dft及fft的性质和特点。（2）利用fft算法计算信号的频谱。相关算法：读入离散序列x(n) 采集长度为n=16的数据，调用matla

6、b中的函数fft（x,16）与fft（x,64）对其作离散傅里叶变换得到16点、64点的频谱 采集数据长度为n=64，调用matlab中的函数fft（x,46）对其作离散傅里叶变换得到64点的频谱。3：相关实现代码n=0:1:15;n1=n;x1=cos(0.48*pi*n)+cos(0.52*pi*n); g1=abs(fft(x1,16); subplot(3,2,1);stem(x1);title(x1); subplot(3,2,2);stem(g1);title(g1); n2=n; x2=cos(0.48*pi*n2)+cos(0.52*pi*n2); x2=x2 zeros(1,

7、48);g2=abs(fft(x2,64);subplot(3,2,3);stem(x2);title(x2);subplot(3,2,4);stem(g2);title(g2); n3=n;x3=cos(0.48*pi*n3)+cos(0.52*pi*n3);g3=abs(fft(x3,64); subplot(3,2,5);stem(x3);title(x3); subplot(3,2,6);stem(g3);title(g3); 4：运行结果 1-3：iir滤波器的设计1：基本要求(1)设计一个butterworth数字低通滤波器，设计指标如下： 通带截止频率：0.2，幅度衰减不大于1分

8、贝 阻带截止频率：0.3，幅度衰减大于15分贝(2)分析不同滤波器的特点和结果。(3)编程设计实现iir滤波器。2：实习目的及相关算法实习目的：（1）掌握不同iir滤波器的性质、特点。（2）通过实验学习如何设计各种常用的iir滤波器，以便在实际工作中能根据具体情况使用iir滤波器。相关算法：输入通带截止频率wp，阻带截止频率ws，通带衰减rp，阻带衰减rs，通过这些数值调用n wn=buttord(wp,ws,rp,rs) 函数计算巴特沃斯数字滤波器的阶数n和截止频率wn，再根据阶数n通过函数b,a=butter(n,wn)，即可得到所要的巴特沃斯滤波器。设计一个正弦波信号，再调用函数a=fi

9、lter(b,a,i)让正弦波信号通过滤波器，得到滤波信号。3:实验代码wp=0.2; ws=0.3; rp=1; rs=15; n wn=buttord(wp,ws,rp,rs) %用于计算巴特沃斯数字滤波器的阶数n和截止频率wn b,a=butter(n,wn); %计算n阶巴特沃斯数字滤波器系统函数分子、分母多项式的系数向量b、a,设计所需的低通滤波器 h,omega=freqz(b,a,512);%返回量h包含了离散系统频响 ，调用中若n默认，默认值为512。 plot(omega/pi,20*log10(abs(h);grid; xlabel(omega/pi);ylabel(gai

10、n,db); title(iir butterworth lowpass filter); wp=0.2; ws=0.3; rp=1; rs=15; n1,wn1=buttord(wp,ws,rp,rs) %用于确定阶次 b,a=butter(n,wn);%用于直接设计巴特沃兹数字滤波器，即为iir滤波器 %freqz(b,a); t=1:300; i=sin(0.1*pi*t)+sin(0.4*pi*t);%设计正弦波 figure;plot(i); a=filter(b,a,i);%正弦波通过滤波器 figure;plot(a); 4:运行结果巴特沃斯滤波器单独运行结果输入为滤波的正弦波信

11、号以及滤波后的效果 1-4：fir低通滤波器的设计 1：基本要求选取合适窗函数设计一个线性相位fir低通滤波器，使它满足如下性能指标：通带截止频率：p=0.5,通带截止频率处的衰减不大于3分贝；阻带截止频率：s=0.66，阻带衰减不小于40分贝2实习目的及相关算法实习目的：（1） 掌握用窗函数法设计fir滤波器的原理和方法。（2） 熟悉线性相位滤波器特性。（3） 了解各种窗函数对滤波器特性的影响。相关算法：通过其通带截止频率p与阻带截止频率s算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率，根据所要求的理想滤波器，得到hd（n）。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不

12、小于40分贝，我选择最接近的汉宁窗，最后调用函数h=hd.*win 及freqz(h,1,512)得到实际汉宁窗的响应和实际滤波器的幅度响应。3：实现代码wp =0.5*pi; ws=0.66*pi; wdelta =ws-wp; %过渡带宽度 n=ceil(8*pi/wdelta) %滤波器长度 if rem(n,2)=0 n=n+1; end nw =n; wc =(wp+ws)/2; %理想低通滤波器的截止频率 n =0: n-1; alpha =(n-1)/2; m =n-alpha+0.00001; hd =sin(wc*m)./(pi*m); %一个响应 win =(hanning

13、(nw); %汉宁窗 h=hd.*win; %实际汉宁窗的响应 freqz(h,1,512); %实际滤波器的幅度响应 4：运行结果 1-5综合实习1：基本要求录制一段自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数法或双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的语音信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，用matlab设计一信号处理系统界面。 2：实习目的及相关算法实习目的：（1）能够对之前四个实习中的内容加以很好的应用（2）锻炼

14、自己的实际能力，使其不在局限简单的小模块相关算法：调用函数function pushbutton1_callback(hobject, eventdata, handles)实现一个信号处理系统界面。选择左键时，用双线性变换法设计滤波器来对信号进行处理，选择右键时，用窗函数法设计滤波器来对信号进行处理。读取语音信号，对语音信号进行f=8000的频率进行采样，调用函数y1=fft(x1,2048)对所采集的点做2048点fft变换。先设计butterworth模拟滤波器，再用双线性变换法实现模拟滤波器到数字滤波器的转换。最后调用函数f1=filter(bz,az,x2)对加了噪声的语音信号进行滤

15、波，得到滤波后的频谱图。3:实现代码function varargout = ts(varargin)% ts m-file for ts.fig% ts, by itself, creates a new ts or raises the existing% singleton*.% h = ts returns the handle to a new ts or the handle to% the existing singleton*.% ts(callback,hobject,eventdata,handles,.) calls the local% function named c

16、allback in ts.m with the given input arguments.% ts(property,value,.) creates a new ts or raises the% existing singleton*. starting from the left, property value pairs are% applied to the gui before ts_openingfunction gets called. an% unrecognized property name or invalid value makes property applic

17、ation% stop. all inputs are passed to ts_openingfcn via varargin.% *see gui options on guides tools menu. choose gui allows only one% instance to run (singleton).% see also: guide, guidata, guihandles% copyright 2002-2003 the mathworks, inc.% edit the above text to modify the response to help ts% la

18、st modified by guide v2.5 22-nov-2013 10:55:52% begin initialization code - do not editgui_singleton = 1;gui_state = struct(gui_name, mfilename, . gui_singleton, gui_singleton, . gui_openingfcn, ts_openingfcn, . gui_outputfcn, ts_outputfcn, . gui_layoutfcn, , . gui_callback, );if nargin & ischar(var

19、argin1) gui_state.gui_callback = str2func(varargin1);endif nargout varargout1:nargout = gui_mainfcn(gui_state, varargin:);else gui_mainfcn(gui_state, varargin:);end% end initialization code - do not edit% - executes just before ts is made visible.function ts_openingfcn(hobject, eventdata, handles, v

20、arargin)% this function has no output args, see outputfcn.% hobject handle to figure% eventdata reserved - to be defined in a future version of matlab% handles structure with handles and user data (see guidata)% varargin command line arguments to ts (see varargin)% choose default command line output

21、 for tshandles.output = hobject;% update handles structureguidata(hobject, handles);% uiwait makes ts wait for user response (see uiresume)% uiwait(handles.figure1);% - outputs from this function are returned to the command line.function varargout = ts_outputfcn(hobject, eventdata, handles) % vararg

22、out cell array for returning output args (see varargout);% hobject handle to figure% eventdata reserved - to be defined in a future version of matlab% handles structure with handles and user data (see guidata)% get default command line output from handles structurevarargout1 = handles.output;% - exe

23、cutes on button press in pushbutton1.function pushbutton1_callback(hobject, eventdata, handles)% hobject handle to pushbutton1 (see gcbo)% eventdata reserved - to be defined in a future version of matlab% handles structure with handles and user data (see guidata)x1,fs=wavread(语音.wav);wavplay(x1,fs);

24、% - executes on button press in pushbutton2.function pushbutton2_callback(hobject, eventdata, handles)% hobject handle to pushbutton2 (see gcbo)% eventdata reserved - to be defined in a future version of matlab% handles structure with handles and user data (see guidata)fs=8000; %语音信号采样频率为8000 x1,fs=

25、wavread(语音.wav); t=(0:length(x1)-1)/8000; y1=fft(x1,2048); %对信号做2048点fft变换 f=fs*(0:1023)/2048; figure(1);subplot(2,2,1);plot(t,x1); %做原始信号的时域波形 grid on;axis tight; title(原始语音信号); xlabel(time(s); ylabel(幅度); subplot(2,2,2); plot(f,abs(y1(1:1024); %做原始信号的fft频谱 grid on;axis tight; title(原始语音信号的fft频谱) ;

26、xlabel(hz); ylabel(幅度); %双线性变换法设计的巴特沃斯滤波器 a1=0.05;a2=0.10; d=a1*cos(2*pi*3800*t)+a2*sin(2*pi*3600*t); d(:,2)=d(:,1);x2=x1+d; wp=0.8*pi; ws=0.85*pi; rp=1; rs=15; fs=8000; ts=1/fs; wp1=2/ts*tan(wp/2); %将模拟指标转换为数字指标 ws1=2/ts*tan(ws/2); n,wn=buttord(wp1,ws1,rp,rs,s); %选择滤波器最小阶数 z,p,k=buttap(n); %创建butte

27、rworth模拟滤波器 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,2,3) ;plot(w*fs/(2*pi),abs(h) ;grid on;axis tight; xlabel(频率（hz）) ;ylabel(频率响应) ;title(butterworth) ;f1=filter(bz,az,x2); figure(2) ;subplot(2,2,1); plot(t,x2) ; %画

28、出滤波前的时域图 grid on;axis tight; title(滤波前的时域波形); subplot(2,2,2); plot(t,f1); %画出滤波后的时域图 grid on;axis tight; title(滤波后的时域波形); y3=fft(f1,2048); y2=fft(x2,2048); subplot(2,2,3); plot(f,abs(y2(1:1024); %画出滤波前的频谱图 grid on;axis tight; title(滤波前的频谱) ;xlabel(hz); ylabel(幅度); subplot(2,2,4) plot(f,abs(y3(1:1024

29、); %画出滤波后的频谱图 grid on;axis tight; title(滤波后的频谱) ;xlabel(hz); ylabel(幅度); wavplay(f1,fs);% - executes on button press in pushbutton3.function pushbutton3_callback(hobject, eventdata, handles)% hobject handle to pushbutton3 (see gcbo)% eventdata reserved - to be defined in a future version of matlab%

30、handles structure with handles and user data (see guidata)fs=8000; %语音信号采样频率为8000 x1,fs=wavread(语音.wav); t=(0:length(x1)-1)/8000; y1=fft(x1,2048); %对信号做2048点fft变换 f=fs*(0:1023)/2048; figure(1) ;subplot(2,1,1); plot(t,x1) ; %做原始信号的时域波形 grid on;axis tight; title(原始语音信号); xlabel(time(s); ylabel(幅度); su

31、bplot(2,1,2) ;plot(f,abs(y1(1:1024) ; %做原始信号的fft频谱 grid on;axis tight; title(原始语音信号的fft频谱) ;xlabel(hz); ylabel(幅度); %窗函数设计滤波器 t=(0:length(x1)-1)/8000; f=fs*(0:2047)/4096; a1=0.05;a2=0.10; d=a1*cos(2*pi*3600*t)+a2*sin(2*pi*3800*t); d(:,2)=d(:,1);x2=x1+d; wp=0.8*pi; ws=0.85*pi; wdelta=ws-wp; n=ceil(6.

32、6*pi/wdelta); %取整 wn=(0.8+0.85)*pi/2; bz,az=fir1(n,wn/pi,hamming(n+1); %选择窗函数，并归一化截止频率 figure(2) ;freqz(bz,az); grid on;axis tight; f2=filter(bz,az,x2); figure(3) ;subplot(2,2,1); plot(t,x2); grid on;axis tight; title(滤波前的时域波形); subplot(2,2,2) ;plot(t,f2); grid on;axis tight; title(滤波后的时域波形); y3=fft

33、(f2,4096); f=fs*(0:2047)/4096; y2=fft(x2,4096); subplot(2,2,3); plot(f,abs(y2(1:2048); grid on;axis tight; title(滤波前的频谱) ;xlabel(hz); ylabel(幅度); subplot(2,2,4) ;plot(f,abs(y3(1:2048); grid on;axis tight; title(滤波后的频谱); xlabel(hz); ylabel(幅度); wavplay(f2,fs);4：运行结果按键显示结果选择左边的双线性按键结果显示如下：选择右边的按键结果显示如

34、下：汉明窗显示如下：1-6 频率抽选1. 流程图2. 程序代码 主函数#include #include #include#include dsp.hvoid main() complex *x,*xtemp; int i,nn,mm; double y,q; system(cls); printf(n*); printf(n* filename :fftl.c *); printf(n* function :fft analysis for signals *); printf(n* signal :x(n) = exp(-n) *); printf(n*); printf(n input

35、m = ); scanf(%d,&mm); nn = 1mm; x = (complex*)calloc(nn,sizeof(complex); xtemp = (complex*)calloc(nn,sizeof(complex); x0.real = 1.0; x0.imag = 0.0; for(i = 1;inn;i+) xi.real = (float)exp(-i);xi.imag = 0.0; for(i = 0; i nn ; i+)xtempi = xi; /* fft analysis for xn*/ fft(x,mm); printf(finished! nn); /*

36、 output result in a format */ printf( k xr(k) xi(k) |x(k)| q(k) n); printf(-n);for(i = 0;i0) q = 0.0;if(xi.real0) q = 4.0*(atan(1.0); q = atan(xi.imag/xi.real);printf(%4d %9f %9f %9f %9fn,i,xi.real,xi.imag,sqrt(y),q); printf(nn);子函数#include /#include #include /#include/* define struct complex */type

37、def structfloat real;float imag;complex;extern void fft (complex *x,int m);extern void ifft (complex *x,int m);void fft(complex *x,int m)static complex *w;static int mstore = 0;static int n = 1;complex u,temp,tm;complex *xi,*xip,*xj,*wptr;int i,j,k,l,le,windex;double arg,w_real,w_imag,wrecur_real,wr

38、ecur_imag,wtemp_real;if (m != mstore)if (mstore != 0) free(w);mstore = m;if (m = 0)return;n = 1m;le = n/2;w = (complex *)calloc(le-1,sizeof(complex);if (!w)printf(nunable to allocate complex w arrayn);exit(1);arg = 4.0*atan(1.0)/le;wrecur_real = w_real = cos(arg);wrecur_imag = w_imag = -sin(arg);xj

39、= w;for (j = 1;j real = (float)wrecur_real;xj-imag = (float)wrecur_imag;xj+;wtemp_real = wrecur_real * w_real - wrecur_imag *w_imag;wrecur_imag = wrecur_real * w_imag + wrecur_imag *w_real;wrecur_real = wtemp_real;/* 开始fft变换 */le = n;windex = 1;for (l = 0;l m;l+)le =le/2;/* first interation with no multip