数字信号实验报告 (全).doc

_数字信号处理实验报告实验一：用 FFT 做谱分析一、 实验目的 1、进一步加深 DFT 算法原理和基本性质的理解。 2、熟悉 FFT 算法原理和 FFT 子程序的应用。 3、学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用 FFT。二、实验原理用FFT对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是2/ND。可以根据此时选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。对模拟信号的频谱时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。三、实验内容和步骤对以下典型信号进行谱分析： 对于以上信号，x1(n)x5(n) 选择FFT的变换区间N为8和16 两种情况进行频谱分析。分别打印其幅频特性曲线。 并进行对比、分析和讨论;；x6(t)为模拟周期信号，选择 采样频率，变换区间N=16,32,64 三种情况进行谱分析。分别打印其幅频特性，并进行分析和讨论。四、实验代码如下所示：clear; %清除变量 close all; %关闭全部绘图窗口 b=menu(请选择信号x1(n)-x6(n),x1(n)=R4(n),x2(n)=1 2 3 4 4 3 2 1,. x3(n)=4 3 2 1 1 2 3 4,x4(n)=cos(npi/4),x5(n)=sin(npi/8),. x6(n)=cos(8pit)+cos(16pit)+cos(20pit),Exit); i=0; A=8,16,32,64; while(b=7) %当选择EXIT时，返回值7，则退出循环 if b=6 m=menu(请选择FFT变换区间长度N,N=16,N=32,N=64); N=A(m+1); fs=64; n=0:(N-1); x=cos(8*pi*n/fs)+cos(16*pi*n/fs)+cos(20*pi*n/fs); else m=menu(请选择FFT变换区间长度N,N=8,N=16,N=32); N=A(m); n=0:(N-1); if b=1 x=1,1,1,1,0,0,0,0,zeros(1,N-8); elseif b=2 x=1,2,3,4,4,3,2,1,zeros(1,N-8); elseif b=3 x=4,3,2,1,1,2,3,4,zeros(1,N-8); elseif b=4 x=cos(n*pi/4); elseif b=5 x=sin(n*pi/8); end end %先画出信号源图 i=i+1; figure(i); %创建绘图窗口 subplot(2,2,1); %指定1号子图 xlabel(n); %标记X坐标 stem(n,x,.r); ylabel(x(n); title(x,num2str(b),(n)的波形); %进行FFT f=fft(x,N); %再画出FFT波形 subplot(2,2,3); stem(n,abs(f),.b); xlabel(k); ylabel(|X(k)|); title(x,num2str(b),(n)的N=,num2str(N),点FFT); b=menu(请选择信号x1(n)-x6(n),x1(n)=R4(n),x2(n)=1 2 3 4 4 3 2 1,. x3(n)=4 3 2 1 1 2 3 4,x4(n)=cos(npi/4),x5(n)=sin(npi/8),. x6(n)=cos(8pit)+cos(16pit)+cos(20pit),Exit); end close all;五、实验结果图及分析1、实验结果图分析：（1）x1(n)的波形如图1-1、图1-2和图1-3所示，由3张图可知道,N值越大，频率分辨率越高。（2）x2(n)的波形如图1-4、图1-5和图1-6所示，由3张图可知道,N值越大，频率分辨率越高。（3）x3(n)的波形如图1-7、图1-8和图1-9所示，由3张图可知道,N值越大，频率分辨率越高。（4）x4(n)的波形如图1-10、图1-11和图1-12所示。根据参数可得出X4(t)的频率 f=8Hz，当N=8、16、32时，频率分辨率为F0=f s/N=8Hz、4Hz、2Hz，因此在FFT图里分别在N=1、2、4有高幅值,因为截取的为周期序列的整数倍，所以所得出的谱正确。（5）x5(n)的波形如图1-13、图1-14和图1-15所示。根据参数可得出X5(t)的频率 f=4Hz。当N=8时，频率分辨率F0=f s/N=8Hz，因为截取的不是为周期序列的整数倍，而且频率分辨率不够，所得出的谱有较大的误差，所以FFT图包含一些频率分量，不能清楚看清原信号的频率f。当N=16及32时，频率分辨率F0=f s/N=4Hz、2Hz，因此在FFT图里在N=1、2有高幅值,因为截取的为周期序列的整数倍，所以所得出的谱正确。（6） x6(n)的波形如图1-16、图1-17和图1-18所示。根据参数可得出X6(t)里包含3个频率，分别为f1=4，f2=8，f3=10。当N=16,频率分辨率F0=f s/N=4Hz，因为截取的不是x6里各周期序列的整数倍，所得出的谱有频谱泄漏，FFT图里可以看出信号cos(8pit)和cos(16pit)的频率f1=4,f2=8（在点N=1，2处有较大的幅值），而且频率分辨率不够高，不能分辨开第三个信号cos(20pit)的频率f3。当N=32，频率分辨率F0=f s/N=2Hz，因此在FFT图里的点N=2有高幅值，在N=4有高幅值 N=5也有高幅值。因为截取的为周期序列的整数倍，所以所得出的谱正确。当N=64，频率分辨率F0=f s/N=1Hz，因此在FFT图里的点N=4有高幅值，在N=8有高幅值 N=10也有高幅值。因为截取的为周期序列的整数倍，所以所得出的谱正确。变换区间N=64 时频谱幅度是变换区间N=32 时2倍，这种结果正好验证了用FFT对中期序列谱分析的理论。2、 误差分析误差产生的原因：（1）对周期序列的截取不当，造成频谱泄漏（2）抽样点数N太少，频率分辨率不够用FFT做谱分析时参数的选择:（1） 抽样频率要满足奈奎斯特准则，不小于信号最高频率的2倍（2） 在抽样频率一定的情况下，抽样点数N要适当。太小会造成频率分辨力不够，太大会造成数据冗余。对周期序列，最好截取周期的整数倍进行谱分析图1-1：x1(n)的波形:N=8图1-2：x1(n)的波形N=16图1-3：x1(n)的波形N=32图1-4：x2(n)的波形N=8图1-5：x2(n)的波形:N=16图1-6：x2(n)的波形N=32图1-7：x 3(n)的波形N=8图1-8：x 3(n)的波形N=16图1-9：x 3(n)的波形N=32图1-10：x 4(n)的波形N=8图1-11：x 4(n)的波形N=16图1-12：x 4(n)的波形N=32图1-13：x 5(n)的波形N=8图1-14：x 5(n)的波形N=16图1-15：x 5(n)的波形N=32图1-16： x 6(n)的波形N=16图1-17： x 6(n)的波形N=32图1-18： x 6(n)的波形N=64实验二：用双线性变换法设计 IIR 数字滤波器一、实验目的 1、熟悉用双线性变换法设计 IIR 数字滤波器的原理与方法。 2、掌握数字滤波器的计算机仿真方法。 二、实验内容及步骤1、用双线性变换法设计一个butterworth低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2时，最大衰减小于1dB，在阻带内0.3,频率区间上，最小衰减大于15dB。2、打印出数字滤波器在频率区间0,上的幅频响应特性曲线。 3、运用MATLAB产生两个正弦信号，信号频率为50Hz和400Hz，采样频率为1000Hz。两个正弦信号相叠加为输入信号y(t)。设计一滤波器，保留源信号中50Hz的低频信号，对y(t)信号进行滤波。观察滤波前后信号的频谱特性，评价滤波器效果。三、实验步骤1、双线性变换法设计 butterworth 低通 IIR 数字滤波器复习有关butterworth模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数。其中满足本实验要求的数字滤波器系统函数为： （2.1）式中： （2.2）根据设计指标，调用MATLAB信号处理工具箱buttord和butter，也可以得到。由公式（2.1）和（2.2）可见，滤波器由三个二阶滤波器、和级联而成，如图2-1所示。图2-1滤波器的组成2、编写滤波器仿真程序，计算 H(z)对心电信号采样序列 x(n)的响应序列 y(n)。 设 yk(n)为第k 级二阶滤波器 Hk(z)的输出序列，yk1(n)为输入序列，如图 1所示。由式2.2 可以得到差分方程：当 k=1 时，yk-1n=x(n)。所以 H(z)对 x(n)的总响应序列 y(n)可以用顺序迭代算法得到。即以此对 k=1,2,3，求解差分方程式 2.3，最后得到y3n=y(n) 。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。也可以直接调用Matlab 的 filter 库函数实现仿真。 四、实验代码。1、双线性变换法设计 butterworth 低通 IIR 数字滤波器的代码如下所示：T=1;Fs=1/T;wpz=0.2;wsz=0.3;wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=15;N,wc=buttord(wp,ws,rp,rs,s);B,A=butter(N,wc,s);fk=0:1/512:1;wk=2*pi*fk;Hk=freqs(B,A,wk);subplot(2,2,1);plot(fk,20*log10(abs(Hk);grid on;xlabel(f/Hz);ylabel(幅度(dB);axis(0,1,-100,5);title(模拟滤波器的幅频响应曲线);N,wc=buttord(wpz,wsz,rp,rs);Bz,Az=butter(N,wc);wk=0:pi/512:pi;Hz=freqz(Bz,Az,wk);subplot(2,2,4);plot(wk/pi,20*log10(abs(Hz);grid on;xlabel(omega/pi);ylabel(幅度(dB);axis(0,1,-100,5);title(数字滤波器的幅频响应曲线);2、根据叠加输入信号y(t)设计的滤波器代码如下：wpd=0.2*pi; %滤波器的通带截止频率 wsd=0.3*pi; %滤波器的阻带截止频率 Rp=1;As=15; %输入滤波器的通阻带衰减指标 Fs=1;T=1/Fs;wp=(2/T)*tan(wpd/2); ws=(2/T)*tan(wsd/2); n,wc=buttord(wp,ws,Rp,As,s) z0,p0,k0=buttap(n); ba=k0*poly(z0); aa=poly(p0); ba1,aa1=lp2lp(ba,aa,wc);%变换为模拟低通滤波器%用双线性变换法计算数字滤波器系数 bd,ad=bilinear(ba1,aa1,Fs) %双线性变换 %求数字系统的频率特性 H,w=freqz(bd,ad); dbH=20*log10(abs(H)/max(abs(H);%化为分贝值 subplot(3,2,1),plot(w/pi,abs(H);grid on;ylabel(|H|);xlabel(w/);title(幅度响应);axis(0,1,0,1.1);subplot(3,2,2),plot(w/pi,dbH);grid on;title(幅度响应(dB);ylabel(dB);xlabel(w/);axis(0,1,-40,5);fs=1000;t=0:1/fs:0.1;x1=sin(2*pi*50*t);x2=sin(2*pi*400*t);y=x1+x2;subplot(3,2,3),plot(t,y);title(滤波前信号);N=1024;T=1/fs;k=0:N-1;fw=fft(y,N);subplot(3,2,4);plot(k,abs(fw);title(滤波前频谱);sf=filter(bd,ad,y);%用所设计的滤波器进行滤波subplot(3,2,5),plot(t,sf);title(滤波后信号);fw1=fft(sf,N);subplot(3,2,6),plot(k,abs(fw1);title(滤波后频谱);五、 实验结果图及分析1、 双线性变换法设计 butterworth 低通 IIR 数字滤波器的实验图像如图2-2所示：图2-2：模拟和数字滤波器的幅频响应图像2、 根据叠加输入信号y(t)设计的滤波器实验图像如2.3所示：图2-3：滤波器的幅频响应图像及y(t)滤波前后的图像双线性变换法的特点：（1）模拟滤波器经过双线性变换后，不存在频率特性的混叠失真，因而对模拟滤波器的频率响应函数Ha(s) 无限带要求，而且能够直接用于设计低通、高通、带通、带阻等各种类型的数字滤波器。（2）与冲激响应不变法中模拟频率与数字频率之间的线性关系=Ts不同的是，双线性变换法中模拟滤波器的频率与所转换成的数字滤波器的频率之间是非线性关系，但是，如果事先进行频率预畸变，这种非线性关系不会使所设计的数字滤波器的幅频特性受到影响。由于频率的非线性关系会产生相频特性失真，所以若对数字滤波器的相位特性要求较严，则不宜采用这种变换方法。（3）双线性变换方法比较容易，不需要将模拟系统函数进行部分分式分解。实验三：语音信号谱分析及去噪处理一、实验目的 1、通过对实际采集的语音信号进行分析和处理，获得数字信号处理实际应用的感性认识。 2、掌握数字信号谱分析的知识。 3、掌握数字滤波器设计的知识，并通过对语音信号的去噪处理，获得数字滤波器实际应用的感性知识。二、实验内容 1、用麦克风自行采集两段语音信号高频噪声、人声+高频噪声（.wav格式）。 2、通过 Matlab 读入采集信号，观察其采样频率，并绘图采样信号。 3、通过 Matlab 对语音信号进行谱分析，分析出噪声的频带。 3、 设计一滤波器，对叠加入噪声的语音信号进行去噪处理。绘图并发声去噪后的信号。三、实验原理1、FFT原理：运用快速傅里叶变换得信号的频谱，快速傅里叶变换是可以用计算机编程实现的一种快速求得信号频谱的方法。在MATLAB编程中可以直接调用FFT函数。2、FIR数字滤波器的设计步骤：（1）将给定的数字滤波器性能指标转化成相应的模拟滤波器性能指标。（2）将模拟滤波器的性能指标变换成模拟低通滤波器的性能指标。（3）用所得的模拟低通滤波器的性能指标，利用某种模拟滤波器逼近方法，设计得出该滤波器。低通滤波器的频率响应：设计滤波器时主要考虑四个指标：通带截止频率、阻带起始频率、通带衰减、阻带衰减。四、实验步骤1、利用麦克风采集一段 5s 以内的语音信号。利用格式工厂软件对语音信号进行预处理。通常语音信号为单声道，采样频率为 8000Hz，语音信号为.wav 格式。 2、通过 Matlab 读入语音信号及其采样频率（使用 Matlab库函数 wavread）， 在Matlab 软件的 workspace 工作平台上观察读入的语音信号，在 Matlab 中，对入的语音信号为一维矩阵。应注意，库函数 wavread 自动将语音信号幅度归一化-1,1区间范围。使用 Matlab 库函数 plot 绘图语音信号，并使用库函数 sound 发音语音信号。 3、分析噪声的频谱。在这里进行谱分析的目的，是了解噪声信号的频谱特性，为去噪滤波器的技术指标提供依据。叠加入噪声的语音信号，如图 3.1 所示： （4）通过 Matlab 对语音信号进行谱分析。应注意，对信号进行谱分析，在实验一中已经详细介绍过。在这里进行谱分析的目的，是了解本段语音信号的频谱特性，为去噪滤波器的技术指标提供依据。 （5）根据语音信号及噪声信号的频谱特性，自行设计一滤波器，对叠加入噪声的语音信号进行去噪处理。最后绘图并发声去噪后的信号。应注意，数字滤波器的实际应考虑实际需求，合理制定滤波器的技术指标。5、 实验内容MATLAB信号处理工具箱函数buttp buttor butter是巴特沃斯滤波器设计函数，其有5种调用格式，本实验中用到的是N,wc=butter(N,wc,Rp,As,s),该格式用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率wc。函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。 1、设计流程图语音信号滤波去噪使用脉冲不变响应法设计的巴特沃斯滤波器的设计流程如图所示：开始语音信号的采集（wavread函数），画时域图快速傅里叶变换，并且画频谱图设定滤波器性能指标，通带截止频率fb=1100，阻带截止频率fc=1200，通带波纹Ap=1，阻带波纹As=20脉冲响应不变法法设计巴特沃斯滤波器验证并进行频谱分析设计好的滤波器进行滤波处理比较滤波前后语音信号的波形及频谱回放语音信号结束图3-2脉冲响应不变法巴特沃斯滤波器对语音信号去噪流程图2、语言信号的采集 下图是基于PC机的语音信号采集过程:通过计算机录音系统录制5秒： 图3-3语音采集过程3、语音信号的频谱分析如图所示为原始语音信号图3-3：原始语音信号图像图3-4：含有噪声的原始语音图像图3-5含有噪声的语音频谱图4、滤波器设计将数字滤波器的设计指标设为通带截止频率fb=1100HZ,阻带频率fc=1200HZ,通带波纹Ap=1dB，阻带波纹As=20dB,要求确定H(z)。设计程序如下：Fp=1200;Ft=8000;As=20;Ap=1;wp=2*pi*Fp/Ft;ws=2*pi*Fs/Ft;fp=2*Ft*tan(wp/2);fs=2*Ft*tan(ws/2);n,wn=buttord(wp,ws,Ap,As,s); b,a=butter(n,wn,s); num,den=bilinear(b,a,1); h,w=freqz(num,den); Buttor波形如下图3.4：图3-6 ：IIR低通滤波器图像5、 信号滤波过程（1）原始代码如下：y,fs,nbits=wavread(E:USERbinYuyin2.wav);%载入波形N=length(y)%计算输入波形长度Y=fft(y,N);%进行快速傅里叶变换%sound(y);%播放声音figure(5);plot(abs(Y);%滤波前频谱图Fp=1200;%阻带截止频率Fs=1100;%通带截止频率Ft=8000;%采集频率As=20;%通带波纹Ap=1Ap=1;%阻带波纹As=20wp=2*pi*Fp/Ft;ws=2*pi*Fs/Ft;fp=2*Ft*tan(wp/2);fs=2*Ft*tan(w