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

中国地质大学（武汉）中国地质大学（武汉）数字信号处理上机实习数字信号处理上机实习学生姓名：学生姓名：班班级：级：071132071132学学号：号：20131002013100指导老师：王晓莉指导老师：王晓莉题目一题目一离散卷积计算离散卷积计算一、实验题目一、实验题目设线性时不变（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 实现离散卷积计算。三、算法设计三、算法设计离散卷积定义为：y(n) kx(k)h(nk)n1、h(n) = (0.8)n，0 n 4，x(n) = u(n)-u(n-4)，y(n) x(n)h(n) mx(m)h(n m)(a) 当n 0时，y(n) 0；(b) 当0 n 3时,y(n) (0.8)n；m0n(c) 当4 n 7时，y(n) mn34(0.8)n；(d) 当n 7时，y(n) 0；2、h(n) (0.8) u(n)，x(n) u(n)-u(n4)，y(n) x(n)h(n) (a) 当n 0时，y(n) 0；nmx(m)h(n m)(b) 当0 n 3时,y(n) (0.8)n；m0n(c) 当4 n 20时，y(n) mn320n(0.8)n；(d) 当21n23时，y(n) (e) 当n 23时，y(n) 0；mn3(0.8)n；3、h(n) (0.8) u(n)，x(n) u(n)，y(n) x(n)h(n) (a) 当n 0时，y(n) 0；(b) 当0 n 70时,y(n) (0.8)n；m0nnmx(m)h(n m)(c) 当71 n140时，y(n) (d) 当n 140时，y(n) 0；mn7070(0.8)n；四、程序分析四、程序分析所用到的函数：（1）y=conv（x.,h）：卷积运算函数，计算y(n) x(n)*h(n);（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:4h=0.8.nx=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);六、运行结果六、运行结果图中从左至右三列依次对应x(n)、h(n)及卷积结果y(n)题目二题目二离散傅立叶变换离散傅立叶变换一、实验题目一、实验题目设有离散序列x(n) cos(0.48n)cos(0.52n)分析下列三种情况下的幅频特性。（1） 采集数据长度 N=16，分析16 点的频谱，并画出幅频特性。采集数据长度N=16，并补零到 64 点，分析其频谱，并画出幅频特性。（2） 采集数据长度 N=64，分析 64 点的频谱，并画出幅频特性。观察三幅不同的幅频特性图，分析和比较它们的特点及形成原因。二、实验目的二、实验目的1、理解掌握 DFT 及 FFT 算法；2、利用 FFT 算法计算信号的频谱。三、算法设计三、算法设计当抽样数 N=2 时，以下为蝶形算法图。x(0)x(1)0WNMX(0)W0Nx(2)10WNX(4)X(2)X(6)X(1)X(5)X(3)X(7)x(3)x(4)1W2N110WNWN011WN1WN21W3N1x(5)x(6)x(7)111WN0W2N0WN11、当 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.48 n)+cos(0.52 n)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 仅具有计算上的意义，并不是真正的、物理意义上的频谱。频谱分辨率的提高只能通过提高采样频率实现。第三幅图形取样点数较多，分辨率最高，而第一幅图形则显得较为稀疏，主要是因为取样点数太少。题目三题目三IIRIIR 滤波器的设计滤波器的设计一、实验题目一、实验题目1、设计一个切比雪夫数字低通滤波器，设计指标如下：通带截止频率：0.2 ，幅度衰减不大于 3 分贝阻带截止频率：0.3 ，幅度衰减大于 20 分贝2、分析不同滤波器的特点和结果。3、编程设计实现 IIR 滤波器。二、实验目的二、实验目的1.理解和掌握不同 IIR 滤波器的性质、特点。2.掌握用 Matlab 设计切比雪夫滤波器和椭圆滤波器的方法。三、算法设计三、算法设计1 1、滤波器类型、滤波器类型切比雪夫滤波器：切比雪夫滤波器，又名“车比雪夫滤波器”，是在通带或阻带上频率响应幅度等波纹波动的滤波器。2 2、变换方法、变换方法（a）冲激响应不变法冲激响应不变法的基本原理是从滤波器的冲激响应出发，对模拟滤波器冲激响应 h(t)进行取样，所得到的离散序列 h(nT)作为数字滤波器的单位取样响应。H(z)是由 H(s)通过下式的对应关系得到。11s sk1espkTz1（b） 双线性变换是在所得到满足性能指标要求的模拟滤波器的基础上,通过变换1 1 z1s T 1 z1,从而得到相应的数字滤波器。四、程序分析四、程序分析所用到的函数：（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 滤波器滤波后，由傅里叶变换后的频谱图可看出高频信号被滤除，低频信号被保留了下来。题目四题目四FIRFIR 滤波器的设计滤波器的设计一、实验内容一、实验内容选取合适窗函数设计一个线性相位 FIR 低通滤波器，使它满足如下性能指标：通带截止频率： p=0.66 ,通带截止频率处的衰减不大于 3 分贝；阻带截止频率： s=0.5 ，阻带衰减不小于 40 分贝二、实验目的二、实验目的1、掌握用窗函数法设计 FIR 滤波器的原理和方法。2、熟悉线性相位滤波器特性。3、了解各种窗函数对滤波器特性的影响。三、原理及算法概要三、原理及算法概要1 1、算法、算法通过其通带截止频率 p与阻带截止频率 s算出其过渡带的宽度与滤波器的长度，从而得到理想滤波器的截止频率，根据所要求的理想滤波器，得到 hd（n）。由于其通带截止频率处的衰减不大于3分贝与阻带衰减不小于40分贝，我选择最接近的汉宁窗， 最后调用函数h=hd.*win 及freqz(h,1,512)得到实际汉宁窗的响应和实际滤波器的幅度响应。2 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加窗后的频谱H(ej) Hd(ej)W(ej) 1jj()W(e)H(e)加窗后使2实际频响偏离理想频响，影响主要有两个方面：（1） 通带和阻带之间存在过渡带，过渡带宽度取决于窗函数频响的主瓣宽度。（2） 通带和阻带区间有纹波，这是由窗函数的旁瓣引起的，旁瓣越多，纹波越多。增加窗函数的宽度 N，其主瓣宽度减小，但不改变旁瓣的相对值。为了改善滤波器的性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；旁瓣衰减尽可能大，数量尽可能大，从而改善纹波状况，使实际频响H（ej）更好地逼近理想频响 Hd（ej）。除了矩形窗外，一般还可以采用以下几种窗函数22)WR()N 1N 122海明窗W(n) 0.54WR() 0.23WR()WR()N 1N 122布来克曼窗W(n) 0.42WR() 0.23WR()WR()N 1N 1汉宁窗：W(n) 0.5WR()0.25WR(四、程序设计四、程序设计 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;%语音信号采样频率为 8000 x1=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.滤波后信号的时域、频域波形五结果分析五结果分析所以从语音的频谱可以看出，通过低通滤波器之后，滤波后的时域波形对滤波前的时域波形看起来相对光滑一些，滤波后的频谱图，滤除了添加的高频噪声信号，保留了低频信号，语音效果也改善了不少。通过最后的频谱图可以看出，与之前语音信号的频谱图几乎一样。而通过比较双线性变换法的滤波器可以看出，它们的滤波效果相似。