数字信号处理实习报告.doc

中国地质大学（武汉） 数字信号处理上机实习 学生姓名： 班 级：071132 学 号：2013100 指导老师：王晓莉 题目一 离散卷积计算一、 实验题目设线性时不变（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)，显示输入和输出波形。二、 实验目的1.理解和掌握离散卷积计算；2.学习如何用Mtalab实现离散卷积计算。三、算法设计离散卷积定义为：1、 ， (a) 当 时，；(b) 当时,(0.8)n ；(c) 当时，(0.8)n；(d) 当时，；2、， (a) 当 时，；(b) 当时,(0.8)n；(c) 当时，(0.8)n；(d) 当时，(0.8)n；(e) 当时，；3、， (a) 当 时，；(b) 当时,(0.8)n；(c) 当时，(0.8)n；(d) 当时，；四、 程序分析所用到的函数：（1）y=conv（x.,h）：卷积运算函数，计算;（2）n1=0:4：n1取04；（3）subplot(m,n,p):subplot()函数是将多个图画到一个平面上的工具。其中，m表示是图排成m行，n表示图排成n列，也就是整个figure中有n个图是排成一行的，一共m行，如果m=2就是表示2行图。p表示图所在的位置，p=1表示从左到右从上到下的第一个位置。（4）title(content ):title()函数的功能是为当前坐标系添加标题“content ”。五、程序设计n=0:4;h=0.8.n;x=1 1 1 1subplot(331);stem(x); title(x(n);y=conv(x,h);subplot(332);stem(h); title(h(n);subplot(333);stem(y); title(y(n); n=0:40;h=0.8.n;x=1 1 1 1subplot(334);stem(x); title(x(n);y=conv(x,h);subplot(335);stem(h); title(h(n);subplot(336);stem(y); title(y(n); n=0:40;h=0.8.n;x=zeros(1,0),ones(1,40);subplot(337);stem(x); title(x(n);y=conv(x,h);subplot(338);stem(h); title(h(n);subplot(339);stem(y); title(y(n);六、运行结果图中从左至右三列依次对应、及卷积结果题目二 离散傅立叶变换一、 实验题目设有离散序列 分析下列三种情况下的幅频特性。（1） 采集数据长度N=16，分析16点的频谱，并画出幅频特性。采集数据长度N=16，并补零到64点，分析其频谱，并画出幅频特性。（2） 采集数据长度N=64，分析64点的频谱，并画出幅频特性。观察三幅不同的幅频特性图，分析和比较它们的特点及形成原因。二、 实验目的1、理解掌握DFT及FFT算法；2、利用FFT算法计算信号的频谱。三、算法设计当抽样数N=2M时，以下为蝶形算法图。1、 当N=2M时，则要进行M次分解，即进行M级蝶形单元的计算；2、按自然顺序输入，输出是码位倒置；3、每一级包含N/2个基本蝶形运算；4、第L级有2L-1个蝶群，蝶群间隔为N/2L-1；四、程序分析所用到的函数：（1） fft（）：求x的一维傅里叶变换，计算X（k）；使用方法Xk=fft(xn)，例如：N=8;n=0:N-1;xn=4 3 2 6 7 8 9 0;Xk=fft(xn)（2）stem（）:显示函数图像。（3）xlabel():横坐标轴的名称；（4）ylabel():纵坐标轴的名称；五、程序设计%x(n)=cos(0.48n)+cos(0.52n)n=0:15; %产生序列x(n) 取16点x=cos(0.48*pi*n)+cos(0.52*pi*n);subplot(231);stem(n,x);title(采集数据长度N=16); %显示x(n)x1k=fft(x,16);%进行16点傅里叶变换subplot(234);stem(0:15,abs(x1k),.); %显示X（k）xlabel(n);ylabel(X1(k); title(16点傅立叶变换频谱); x=x zeros(1,48); %取16点，补零到64点x2k=fft(x,64);%进行64点傅里叶变换subplot(232);stem(x);title(采集N=16,补零到64);subplot(235);stem(0:63,abs(x2k),.);xlabel(n);ylabel(X2(k);title(64点傅立叶变换频谱); n=0:63; %产生序列x(n)，取64点x=cos(0.48*pi*n)+cos(0.52*pi*n);subplot(233);stem(n,x);title(采集数据长度N=64);x3k=fft(x);%做64点傅里叶变换subplot(236);stem(0:63,abs(x3k),.);xlabel(n);ylabel(X3(k);title(64点傅立叶变换频谱);六、运行结果七结果分析N点DFT的频谱分辨率是2 /N。一节指出可以通过补零观察到更多的频点，但是这并不意味着补零能够提高真正的频谱分辨率。这是因为xn 实际上是x(t) 采样的主值序列，而将xn补零得到的xn 周期延拓之后与原来的序列并不相同，也不是x(t) 的采样。因此是不同离散信号的频谱。对于补零至M点的x的DFT，只能说它的分辨率2 /M仅具有计算上的意义，并不是真正的、物理意义上的频谱。频谱分辨率的提高只能通过提高采样频率实现。第三幅图形取样点数较多，分辨率最高，而第一幅图形则显得较为稀疏，主要是因为取样点数太少。题目三 IIR滤波器的设计一、 实验题目1、设计一个切比雪夫数字低通滤波器，设计指标如下： 通带截止频率：0.2，幅度衰减不大于3分贝 阻带截止频率：0.3，幅度衰减大于20分贝2、分析不同滤波器的特点和结果。3、编程设计实现IIR滤波器。二、实验目的1.理解和掌握不同IIR滤波器的性质、特点。2.掌握用Matlab设计切比雪夫滤波器和椭圆滤波器的方法。三、算法设计1、滤波器类型切比雪夫滤波器：切比雪夫滤波器，又名“车比雪夫滤波器”，是在通带或阻带上频率响应幅度等波纹波动的滤波器。2、变换方法（a）冲激响应不变法冲激响应不变法的基本原理是从滤波器的冲激响应出发，对模拟滤波器冲激响应h(t)进行取样，所得到的离散序列h(nT)作为数字滤波器的单位取样响应。H(z)是由H(s)通过下式的对应关系得到。（b）双线性变换是在所得到满足性能指标要求的模拟滤波器的基础上,通过变换,从而得到相应的数字滤波器。四、 程序分析所用到的函数：（1）N,wn=buttord(Wp,Ws,ap,as,s)：阶数，截止频率=buttord（通带截止频率，阻带截止频率，通带衰减，阻带衰减，模拟）；（2）num,den=bilinear(b,a,fs)：双线性变换可选频率预畸变，在IIR数字滤波器设计中有把模拟的原型滤波器转变成相应的数字滤波器，即通过已知S变换中模拟滤波器的系数求数字滤波器的系数。其中有两种方法来进行这种变换，一种是激冲不变法，另一种是双线性Z变换法。bilinear 函数便是双线性Z变换法。函数中的b和a是原模拟滤波器S变换中的分子和分母系数，num和den是数字滤波器的系数。（3）grid on：增加主要网格线为当前轴；（4）h,w=freqz(num,den)：求频率响应。函数的输出：a.滤波器的频率响应H（N点）b.频率向量W(N点，且单位为弧度)，频率向量W是均匀分布在滤波器的上半区，即：0：pi，这些点上的频率响应都将通过此函数计算出来。；滤波器的系数：分子为num，分母为den。 五、程序设计clc;Rs=20; Rp=3;Wp1=0.2*pi;Ws1=0.3*pi;N,Wn=cheb1ord(Wp1,Ws1,Rp,Rs,s);num1,den1=cheby1(N,Rp,Wn,s);num,den=bilinear(num1,den1,0.5); subplot(2,2,1);zplane(num,den);title(零极点图)w=0:pi/256:pi;h=freqz(num,den,w);g=abs(h);g1=angle(h);subplot(223);plot(w/pi,g);gridaxis(0 1 0 1);xlabel(频率); ylabel(幅度); title(幅频响应);subplot(224);plot(w/pi,g1);grid onaxis(0 1 -5 5);xlabel(频率); ylabel(相位); title(相频响应);%测试t=1:300;I=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t);%设计正弦波figure;subplot(311);plot(I);title(原信号x(n)=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t);A=filter(num,den,I);%正弦波通过滤波器subplot(312);plot(A);title(滤波后y(n);六、 运行结果（1） 零极点图、幅频响应、相频响应（2）输入为I=sin(0.1*pi*t)+sin(0.8*pi*t)+sin(0.95*pi*t)时的滤波结果七、 结果分析由幅频特性可知，椭圆滤波器在通带内的频率特性是平坦的,并且随着频率的增加而衰减。通带截止频率：0.2，幅度衰减不大于3分贝，阻带截止频率：0.3，幅度衰减大于20分贝。正弦信号在经过IIR滤波器滤波后，由傅里叶变换后的频谱图可看出高频信号被滤除，低频信号被保留了下来。题目四 FIR滤波器的设计一、实验内容选取合适窗函数设计一个线性相位FIR低通滤波器，使它满足如下性能指标：通带截止频率：p=0.66,通带截止频率处的衰减不大于3分贝；阻带截止频率：s=0.5，阻带衰减不小于40分贝二、实验目的1、掌握用窗函数法设计FIR滤波器的原理和方法。2、熟悉线性相位滤波器特性。3、了解各种窗函数对滤波器特性的影响。三、原理及算法概要1、算法通过其通带截止频率p与阻带截止频率s算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率，根据所要求的理想滤波器，得到hd（n）。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不小于40分贝，我选择最接近的汉宁窗，最后调用函数h=hd.*win 及freqz(h,1,512)得到实际汉宁窗的响应和实际滤波器的幅度响应。2、原理概要利用窗函数法设计FIR 滤波器，FIR滤波器的最大特点是其相位特性可以设计为严格的线性，而其幅值可以任意设置，这样输出波形就不会相位失真。理想低通滤波器的单位取样响应hd（n）是无限长的，所以要用一个有限长的因果序列h（n）进行逼近，最有效的方法是截断hd（n），即用有限长的窗函数w（n）来截取hd（n），表示为h（n）=hd（n）w（n）。为获得线性相位的FIR滤波器，h（n）必须满足中心对称条件，序列h（n）应有一定的延迟，且=（N-1）/2频率响应逼近hd（ejw）的FIR滤波器，最简单的窗函数为矩形窗： 1 n(N-1)/2加窗后的频谱 加窗后使实际频响偏离理想频响，影响主要有两个方面：（1） 通带和阻带之间存在过渡带，过渡带宽度取决于窗函数频响的主瓣宽度。（2） 通带和阻带区间有纹波，这是由窗函数的旁瓣引起的，旁瓣越多，纹波越多。增加窗函数的宽度N，其主瓣宽度减小，但不改变旁瓣的相对值。为了改善滤波器的性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；旁瓣衰减尽可能大，数量尽可能大，从而改善纹波状况，使实际频响H（ej）更好地逼近理想频响Hd（ej）。除了矩形窗外，一般还可以采用以下几种窗函数 汉宁窗：海明窗布来克曼窗四、 程序设计 clc; wp=0.66*pi;ws=0.5*pi;wd=wp-ws;%hanning窗N_hann=ceil(8*pi/wd); wc=(wp+ws)/2;if mod(N_hann,2)=0 N_hann=N_hann+1;endhann_win=hanning(N_hann);b_hann=fir1(N_hann-1,wc/pi,high,hann_win);freq_axis=0:pi/512:pi-pi/512;freqn=0:511/512;H_hann=freqz(b_hann,1,512);figure(1);subplot(311);plot(freqn,20*log10(abs(H_hann);xlabel(归一化频率w/pi);ylabel(幅度谱（dB）);title(采用hanning窗设计fir滤波器幅度谱);subplot(312);plot(freqn,angle(H_hann);xlabel(归一化频率w/pi);ylabel(相位谱);title(采用hanning窗设计fir滤波器相位谱);subplot(313);stem(b_hann);xlabel(n);ylabel(h(n);title(单位脉冲响应序列);%hamming窗N_hamm=ceil(8*pi/wd);wc=(wp+ws)/2;if mod(N_hamm,2)=0 N_hamm=N_hamm+1;endhamm_win=hamming(N_hamm);b_hamm=fir1(N_hamm-1,wc/pi,high,hamm_win);H_hamm=freqz(b_hamm,1,512);figure(2); subplot(311);plot(freqn,20*log10(abs(H_hamm);xlabel(归一化频率w/pi);ylabel(幅度谱（dB）);title(采用hamming窗设计fir滤波器幅度谱);subplot(312);plot(freqn,angle(H_hamm);xlabel(归一化频率w/pi);ylabel(相位谱);title(采用hamming窗设计fir滤波器相位谱);subplot(313);stem(b_hann);xlabel(n);ylabel(h(n);title(单位脉冲响应序列);五 程序运行结果（1）采用hanning窗设计FIR滤波器（2）采用hamming窗设计FIR滤波器六 结果分析上图为选用汉宁窗和汉明窗的FIR 滤波器的相频特性、幅频特性及单位取样响应hd（n）。由运行得出的幅频特性可知，满足设计指标通带截止频率：p=0.66；阻带截止频率：s=0.5。得出的窗主瓣宽度较宽，旁瓣宽度较窄，也合乎要求。题目五一、实验内容录制一段电话的拨号音，并对录制的信号进行采样；画出采样后信号的时域波形和频谱图；必要的话对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，根据信号的频谱解析出拨号的数字。最后，用MATLAB设计一信号处理系统界面。二、实验目的1、掌握用滤波器处理语音信号的原理和方法。2、掌握用matlab设计界面。三、原理及算法概要调用函数function pushbutton1_Callback(hObject, eventdata, handles)实现一个信号处理系统界面。选择button1时，用双线性变换法设计滤波器来对信号进行处理。读取语音信号，对语音信号进行f=8000的频率进行采样，调用函数y1=fft(x1,2048)对所采集的点做2048点FFT变换。先设计butterworth模拟滤波器，再用双线性变换法实现模拟滤波器到数字滤波器的转换。最后调用函数f1=filter(bz,az,x2)对加了噪声的语音信号进行滤波，得到滤波后的频谱图。四 程序% - Executes 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)clc;fs=8000; %语音信号采样频率为8000x1=wavread(D:matlab2011workshua.wav);t=(0:length(x1)-1)/8000;y1=fft(x1,2048); %对信号做2048点FFT变换f=fs*(0:1023)/2048;plot(handles.axes1,t,x1); %做原始信号的时域波形grid on;axis tight;title(handles.axes1,原始语音信号);xlabel(handles.axes1,time(s);ylabel(handles.axes1,幅度);plot(handles.axes2,f,abs(y1(1:1024); %做原始信号的FFT频谱grid on;axis tight;title(handles.axes2,原始语音信号的FFT频谱);xlabel(handles.axes2,Hz);ylabel(handles.axes2,幅度);%双线性变换法设计的巴特沃斯滤波器A1=0.05;A2=0.10;d=A1*cos(2*pi*3800*t)+A2*sin(2*pi*3600*t);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); %创建butterworth模拟滤波器Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn); bz,az=bilinear(b,a,Fs); %用双线性法实现模拟到数字的转换H,W=freqz(bz,az); %绘制频率响应曲线plot(handles.axes3,W*Fs/(2*pi),abs(H);grid on;axis tight;xlabel(handles.axes3,频率（Hz）);ylabel(handles.axes3,频率响应);title(handles.axes3,Butterworth);f1=filter(bz,az,x2);plot(handles.axes4,t,x2); %画出滤波前的时域图grid on;axis tight;title(handles.axes4,加噪声信号滤波前的时域波形);plot(handles.axes5,t,f1); %画出滤波后的时域图grid on;axis tight;title(handles.axes5,滤波后的时域波形);y3=fft(f1,2048);y2=fft(x2,2048);plot(handles.axes6,f,abs(y2(1:1024); %画出滤波前的频谱图grid on;axis tight;title(handles.axes6,加噪声信号滤波前的频谱);xlabel(handles.axes6,Hz);ylabel(handles.axes6,幅度);plot(handles.axes7,f,abs(y3(1:1024); %画出滤波后的频谱图grid on;axis tight;title(handles.axes7,滤波后的频谱);xlabel(handles.axes7,Hz);ylabel(handles.axes7,幅度);四 实验结果 1. 实验界面2. 原始信号和加噪声信号的时域、频域波形3. 巴特沃斯滤波器4. 滤波后信号的时域、频域波形 五结果分析 所以从语音的频谱可以看出，通过低通滤波器之后，滤波后的时域波形对滤波前的时域波形看起来相对光滑一些，滤波后的频谱图，滤除了添加的高频噪声信号，保留了低频信号，语音效果也改善了不少。通过最后的频谱图可以看出，与之前语音信号的频谱图几乎一样。而通过比较双线性变换法的滤波器可以看出，它们的滤波效果相似。