13年地大数字信号处理要点

1、数字信号处理报告中国地质大学（武汉）学 校:姓 名:班 级:学 号:指导老师:目录1-1 :离散卷积的运算1-2 :离散傅里叶变换及其应用1-3 : IIR滤波器的设计1-4 : FIR低通滤波器的设计1-5 :用MATLAE设计一信号处理系统1-6:频率抽选1-7 :实习心得体会1-1离算卷积的应用1基本要求设线性时不变(LTI)系统的冲激响应为 h(n),输入序列为x(n)1、 h(n)=(0.8)n, Ow nW 4； 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)求以上

2、三种情况下系统的输出y(n),显示输入和输出波形。2：实验目的及算法原理实验目的：(1 )：掌握离散卷积计算机实现。(2)：进一步对离散信号卷积算法的理解。算法原理：把冲激响应h(n)与输入序列x(n)分别输入到程序中，然后调用离散卷积函数y=co nv(x.,h)即可得到所要求的结果。离散卷积定义为 y(n)二x(k)h(n-k)k =.：n当序列为有限长时则y( n) =x(k)h (n-k)k=03 :程序代码x1=1 1 1 1 ;n x1=0:3;h1=1 0.8 0.64 0.8A3 0.8A4; nh1=0:4;y1=co nv(x1,h1);subplot(3,3,1);ste

3、m(nx1,x1);title('序列 x1');xlabel(' n');ylabel('x1( n)');subplot(3,3,2);stem(nh1,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;n h2=0:1:20

4、;h2=(0.8).A nh2;y2=co nv(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)');subplot(3,3,6);stem(y2);title('序列 y2');xlabel(' n');ylabel

5、('y2( n)');nx3=0:1:20;x3=1.a nx3;nh3=0:1:20;h3=(0.8)A nh3;y3=co nv(x3,h3);subplot(3,3,7);stem( nx3,x3);title(' xlabel(' n');ylabel('x3( n)');subplot(3,3,8);stem( nh3,h3);title(' xlabel(' n');ylabel('h3( n)'); subplot(3,3,9);stem(y3);title(' xlabel(

6、' n');ylabel('y3( n)');序列x3');序列h3');序列y3');序列hi序列灯0.5JZ4瘁列cQ -24n序列x2f 224序列h2°0序列yd序列幻°0序列h30 60n4:结果分析1-2离散傅里叶变换及应用1：基本要求设有离散序列x(n)=cos(0.48 n n)+cos(0.52 n n)分析下列三种情况下的幅频特性。(1) 采集数据长度N=16，分析16点的频谱，并画出幅频特性。采集数据长度N=16，并补零到64点，分析其频谱，并画出幅频特性。(2) 采集数据长度N=64，分析64点

7、的频谱，并画出幅频特性。观察三幅不同的幅频特性图，分析和比较它们的特点及形成原因。2:实习目的及相关算法实习目的：(1) 了解DFT及FFT的性质和特点。(2) 利用FFT算法计算信号的频谱。相关算法：读入离散序列x(n) 采集长度为N=16的数据，调用matlab中的函数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

8、(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,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*

9、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 n，幅度衰减不大于1分贝 阻带截止频率：0.3 n ,幅度衰减大于15分贝(2) 分析不同滤波器的特点和结果。(3) 编程设计实现IIR滤波器。2:实习目的及相关算法实习目的：(1 )掌握不同IIR滤波器的性质、特点。(2)通

10、过实验学习如何设计各种常用的IIR滤波器，以便在实际工作中能根据具体情况使用IIR滤波器。相关算法：输入通带截止频率 Wp阻带截止频率 Ws通带衰减Rp,阻带衰减Rs,通过这些数值调用N Wn=buttord(Wp,Ws,Rp,Rs)函数计算巴特沃斯数字滤波器的阶数N和截止频率wn,再根据阶数N通过函数b,a=butter(N,Wn),即可得到所要的巴特沃斯滤波器。设计一个正弦波信号，再调用函数 A=filter(b,a,I)让正弦波信号通过滤波器，得到滤波信号。3:实验代码Wp=0.2;Ws=0.3;Rp=1;Rs=15;N Wn=buttord(Wp,Ws,Rp,Rs) %用于计算巴特沃斯

11、数字滤波器的阶数 N和截止 频率wnb,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 n,dB');title('IIR Butterworth Lowpass Filter');Wp=0.2;Ws=0.3;Rp=1;Rs=15;N1,

12、Wn1=buttord(Wp,Ws,Rp,Rs) % 用于确定阶次 b,a=butter(N,Wn);%用于直接设计巴特沃兹数字滤波器，即为IIR滤波器%freqz(b,a);t=1:300;I=s in (0.1*pi*t)+si n(0.4*pi*t);%设计正弦波figure;plot(I);A=filter(b,a,l);%正弦波通过滤波器figure;plot(A);4:运行结果巴特沃斯滤波器单独运行结果输入为滤波的正弦波信号以及滤波后的效果o o o 2 2mp) apn_c6e 乏-400 11111111111Norma

13、lized Frequency (xz rad/sample)0-200-400-60001Normalized Frequency (xx rad/sample)1-4 : FIR低通滤波器的设计1：基本要求选取合适窗函数设计一个线性相位FIR低通滤波器，使它满足如下性能指标：通带截止频率：3 p=0.5 n，通带截止频率处的衰减不大于3分贝；阻带截止频率：3 s=0.66n，阻带衰减不小于 40分贝2实习目的及相关算法实习目的：(1) 掌握用窗函数法设计 FIR滤波器的原理和方法。(2) 熟悉线性相位滤波器特性。(3) 了解各种窗函数对

14、滤波器特性的影响。相关算法：通过其通带截止频率 3p与阻带截止频率3 s算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率，根据所要求的理想滤波器，得到hd ( n)。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不小于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)=0N=N+1;endNw =N;wc =(wp+

15、ws)/2; %理想低通滤波器的截止频率n =0: N-1;alpha =(N-1)/2;m =n-alpha+0.00001;hd =si n(wc*m)./(pi*m); %一个响应wi n =(ha nnin g(Nw)'% 汉宁窗h=hd.*wi n;%实际汉宁窗的响应freqz(h,1,512);%实际滤波器的幅度响应4 :运行结果mp)-snBEWO0 10 20 3040,50 60 70 00 91Normalized Frequency (x rad/sample)O5o ooo-2O1Normalized Fre

16、quency rad sample)o ooo1-5综合实习1：基本要求录制一段自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频 谱图；给定滤波器的性能指标，采用窗函数法或双线性变换设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的语音信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；最后，用MATLAB设计一信号处理系统界面。2:实习目的及相关算法实习目的：(1) 能够对之前四个实习中的内容加以很好的应用(2) 锻炼自己的实际能力，使其不在局限简单的小模块相关算法：调用函数 fun ctio n p

17、ushbutto n1_Callback(hObject, even tdata, ha ndles)实现一个信号处理系统界面。选择左键时，用双线性变换法设计滤波器来对信号进行处理，选择右键时，用窗函数法设计滤波器来对信号进行处理。读取语音信号，对语音信号进行f=8000的频率进行采样，调用函数y仁fft(x1,2048)对所采集的点做 2048点FFT变换。先设计butterworth模拟滤波器，再用双线性变换法实现模拟滤波器到数字滤波器的转换。最后调用函数 f1=filter(bz,az,x2)对加了噪声的语音信号进行滤波，得到滤波后的频谱图。3:实现代码 function varargo

18、ut = ts(vararg in)% TS M-file for ts.fig% TS, by itself, creates a new TS or raises the exist ing%sin glet on*.% H = TS retur ns the han dle to a new TS or the han dle to%the existi ng sin glet on*.%TS('CALLBACK',hObject,eve ntData,ha ndles,.) calls the local%fun ctio n named CALLBACK in TS.

19、M with the give n in put argume nts.% TS('Property','Value',.) creates a new TS or raises the%existi ng sin glet on*. Starting from the left, property value pairs are% applied to the GUI before ts_Ope ningFunction gets called. An%un recog ni zed property n ame or in valid value makes

20、 property applicati on% stop. All in puts are passed to ts_Ope nin gFc n via vararg in.% *See GUI Optio ns on GUIDE's Tools me nu. Choose "GUI allows only one%in sta nee to run (sin glet on)".% See also: GUIDE, GUIDATA, GUIHANDLES % Copyright 2002-2003 The MathWorks, I nc.% Edit the ab

21、ove text to modify the resp onse to help ts% Last Modified by GUIDE v2.5 22-Nov-2013 10:55:52% Begi n ini tialization code - DO NOT EDIT gui_S in glet on = 1;gui_State = struct('gui_Name', mfile name, .'gui_S in glet on', gui_S in glet on,'gui_Ope nin gFc n', ts_Ope nin gFc n

22、,'gui_OutputFc n', ts_OutputFc n, 'gui_LayoutFc n',.'gui_Callback', );if n argin && ischar(varargi n1) gui_State.gui_Callback = str2fu nc(vararg in 1);endif n argoutvarargout1: nargout = gui_ma infcn( gui_State, vararg in :); elsegui_ma infcn( gui_State, vararg in :);

23、end% End ini tialization code - DO NOT EDIT% - Executes just before ts is made visible.function ts_Ope nin gFc n( hObject, even tdata, han dles, vararg in)% This fun cti on has no output args, see OutputFc n.% hObject han dle to figure% even tdata reserved - to be defi ned in a future version of MAT

24、LAB % han dles structure with han dles and user data (see GUIDATA) % varargin comma nd line argume nts to ts (see VARARGIN)% Choose default comma nd line output for ts han dles.output = hObject;% Update han dles structure guidata(hObject, han dles);% UIWAIT makes ts wait for user respo nse (see UIRE

25、SUME)% uiwait(ha ndles.figure1);% - Outputs from this fun ctio n are returned to the comma nd line, function varargout = ts_OutputFc n( hObject, eve ntdata, han dies)% varargout cell array for retur ning output args (see VARARGOUT); % hObject han die to figure% even tdata reserved - to be defi ned i

26、n a future version of MATLAB % han dies structure with han dies and user data (see GUIDATA)% Get default comma nd line output from han dies structure varargout1 = han dles.output;% - Executes on butt on press in pushbutt onl.function pushbutt on 1_Callback(hObject, eve ntdata, han dies)% hObject han

27、dle to pushbutt onl (see GCBO)% even tdata reserved - to be defi ned in a future version of MATLAB % han dies structure with han dies and user data (see GUIDATA) x1,fs=wavread(' 语音.wav');wavplay(x1,fs);% - Executes on butt on press in pushbutt on2.function pushbutt on 2_Callback(hObject, eve

28、 ntdata, han dies)% hObject handle to pushbutt on2 (see GCBO)% even tdata reserved - to be defi ned in a future version of MATLAB % han dies structure with han dies and user data (see GUIDATA)做原始信号的时域波形fs=8000;%x1,fs=wavread(' 语音.wav'); t=(0:le ngth(x1)-1)/8000;y1=fft(x1,2048);%f=fs*(0:1023)

29、/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);%语音信号采样频率为 8000对信号做2048点FFT变换%做原始信号的FFT频谱grid on; axis tight;title（'原始语音信号的FFT频谱'）;xlabel（'Hz'）;ylabel（'幅度&#

30、39;）;%双线性变换法设计的巴特沃斯滤波器A1=0.05;A2=0.10;d=A1*cos(2*pi*3800*t)+A2*si n( 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*ta n(wp/2);%ws1=2/Ts*ta n(ws/2);N,Wn=buttord(wp1,ws1,Rp,Rs,'s'); % Z,P,K=buttap(N);%Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn);

31、bz,az=bili near(b,a,Fs);%H,W=freqz(b z, az);%subplot(2,2,3);plot(W*Fs/(2*pi),abs(H);grid on; axis tight; xlabel('频率(Hz)'); ylabel('频率响应');title('Butterworth'); f1=filter(bz,a z,x2);figure(2); subplot(2,2,1);plot(t,x2) ;%grid on; axis tight;title('滤波前的时域波形');subplot(2,

32、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);%grid on; axis tight;title(' 滤波后的频谱&#

33、39;);将模拟指标转换为数字指标选择滤波器最小阶数创建butterworth 模拟滤波器用双线性法实现模拟到数字的转换 绘制频率响应曲线画出滤波前的时域图画出滤波后的时域图画出滤波前的频谱图画出滤波后的频谱图xlabel('Hz'); ylabel('幅度'); wavplay(f1,fs);% - Executes on butt on press in pushbutt on3.function pushbutt on 3_Callback(hObject, eve ntdata, han dles)% hObject handle to pushbutt

34、 on3 (see GCBO)% even tdata reserved - to be defi ned in a future version of MATLAB % han dles structure with han dles and user data (see GUIDATA)fs=8000;%x1,fs=wavread(' 语音.wav'); t=(0:le ngth(x1)-1)/8000;y1=fft(x1,2048);%f=fs*(0:1023)/2048;figure(1);subplot(2,1,1);plot(t,x1) ;%grid on; axi

35、s tight;title('原始语音信号');xlabel('time(s)');ylabel('幅度'); subplot(2,1,2); plot(f,abs(y1(1:1024) ;%语音信号采样频率为 8000对信号做2048点FFT变换做原始信号的时域波形做原始信号的FFT频谱grid on; axis tight;title('原始语音信号的FFT频谱');xlabel('Hz');ylabel('幅度');%窗函数设计滤波器t=(0:le ngth(x1)-1)/8000; f=fs

36、*(0:2047)/4096;A1=0.05;A2=0.10;d=A1*cos(2*pi*3600*t)+A2*si n( 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.6*pi/wdelta);%wn=(0.8+0.85)*pi/2;b z, az=fir1(N,w n/pi,hammi ng(N+1);%figure(2);freqz(bz,az);grid on; axis tight;f2=filter(bz,a z,x2);figur

37、e(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(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('

38、;Hz');ylabel('幅度'); subplot(2,2,4);plot(f,abs(y3(1:2048);grid on; axis tight;title(' 滤波后的频谱'); xlabel('Hz');ylabel('幅度'); wavplay(f2,fs);4:运行结果 按键显示结果选择左边的双线性按键结果显示如下:O54204060time(s) Butterworth8000.5115Hz2x卅原始语音信号的F盯频谱15010050100020003000频率(Hz选择右边的按键结果显示如下:111 F

39、 a1>1 r wii|1111ii耳P i i i g111iiiklbUkG1 i i i11滤波前的频谱滤波后的频谱2502004ililii . !I*l4L _4iiiiIiL1i>1>111illid T亠 4 开 Hi .114*1>111 .Ji =1 一。I DU!i li!iI111ilogIII14I uu.I>l il!Il11uif r "' T =ii i id11LI0 0511.52HzX104汉明窗显示如下:Elp)upBcmcow50050-001-15011111111040 50 60

40、.70 8091Normalised Frequency i'xz rad-'sample'：1-6频率抽选w(1www1.流程图X(6)XfO)2.程序代码主函数#in elude <dos.h>#in elude <stdlib.h>#in elude<malloe.h>#in elude "dsp.h"void mai n()COMPLEX *x,*xtemp;in t i, nn ,mm;double y,q;system("els");printf("n*");pr

41、in tf("n* Function :fft an alysis for sig nals *");prin tf("n* Sig nal :x( n) = exp(-n)*");printf("n*");prin tf("n In put m =");prin tf("n* File name :fftl.e工sea nf("%d",&m m);nn = 1<<mm;x = (COMPLEX*)calloc( nn ,sizeof(COMPLEX);xtemp

42、= (COMPLEX*)calloc( nn ,sizeof(COMPLEX);xO.real = 1.0;xO.imag = 0.0;for(i = 1;i< nn ;i+)xi.real = (float)exp(-i);xi.imag = 0.0;for(i = 0; i < nn ; i+)xtempi = xi;/* FFT an alysis for x n */fft(x,mm);prin tf("Fi nished! nn");/* output result in a format */printf(" k XR(K) XI(k) |

43、X(k)| Q(k) n");printf("n");for(i = 0;i< nn ;i+)y = xi.real*xi.real+xi.imag*xi.imag;if(xi.real != 0)if(xi.imag = 0.0)if(xi.real>0) q = 0.0;if(xi.real<0) q = 4.0*(ata n(1.0);q = ata n(xi.imag/xi.real);prin tf("%4d %9f %9f %9f %9fn",i,xi.real,xi.imag,sqrt(y),q); prin t

44、f("nn");子函数#in elude <stdio.h>#in elude <iostream>#in elude <math.h>#in elude<alloe.h>/* define struet COMPLEX */typedef struetfloat real;float imag;COMPLEX;extern void fft (COMPLEX *x,i nt m); extern void ifft (COMPLEX *x,i nt m);void fft(COMPLEX *x,i nt m)statie C

45、OMPLEX *w;statie int mstore = 0;statie int n = 1;COMPLEX u,temp,tm;COMPLEX *xi,*xip,*xj,*wptr;int i,j,k,l,le,wi ndex;double arg,w_real,wm ag,wrecur_real,wrecur_imag,wtemp_real;if (m != mstore)if (mstore != 0) free(w);mstore = m;if (m = 0) return;n = 1<<m;le = n/2;w = (COMPLEX *)calloc(le-1,siz

46、eof(COMPLEX); if (!w)prin tf("nUn able to allocate complex w array' n"); exit(1);arg = 4.0*ata n(1.0)/le; wrecur_real = w_real = cos(arg);wrecur_imag = w_imag = -si n(arg);xj = w;for (j = 1;j < le; j+) xj->real = (float)wrecur_real; xj->imag = (float)wrecur_imag; xj+;wtemp_rea

47、l = wrecur_real * w_real - wrecur_imag *wmagwrecur_imag = wrecur_real * w_imag + wrecur_imag *w_real;wrecur_real = wtemp_real;*开始fft*le = n;win dex = 1;for (l = 0;l < m;l+)le =le/2;*first in terati onwith no multiplies*for ( i = 0; i<n; i = i+ 2*le)xi = x + i;xip = xi + le;temp.real =(xi->r

48、eal + xip->real); temp.imag = (xi->imag + xip->imag); xip->real = (xi->real - xip->real);xip->imag = (xi->imag - xip->imag);*xi = temp;stored/*remai niteratio nsuse*/wptr = w + win dex - 1;for ( j = 1; j < le; j+)u = *wptr;for ( i = j; i< n; i = i + 2 *le )xi = x+1;xip = xi+le;temp.real = (xi->real+xip->real);temp.imag = (xi->imag+x