数字信号实验报告 全.docx
《数字信号实验报告 全.docx》由会员分享，可在线阅读，更多相关《数字信号实验报告 全.docx（29页珍藏版）》请在冰豆网上搜索。
数字信号实验报告全
数字信号处理实验报告
实验一：
用FFT做谱分析
一、实验目的
1、进一步加深DFT算法原理和基本性质的理解。
2、熟悉FFT算法原理和FFT子程序的应用。
3、学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。
二、实验原理
用FFT对信号作频谱分析是学习数字信号处理的重要内容。
经常需要进行谱分析的信号是模拟信号和时域离散信号。
对信号进行谱分析的重要问题是频谱分辨率D和分析误差。
频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是2π/N≤D。
可以根据此时选择FFT的变换区间N。
误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。
周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。
如果不知道信号周期，可以尽量选择信号的观察时间长一些。
对模拟信号的频谱时，首先要按照采样定理将其变成时域离散信号。
如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。
三、实验内容和步骤
对以下典型信号进行谱分析：
对于以上信号，x1（n）~x5（n）选择FFT的变换区间N为8和16两种情况进行频谱分析。
分别打印其幅频特性曲线。
并进行对比、分析和讨论;；x6（t）为模拟周期信号，选择采样频率
，变换区间N=16,32,64三种情况进行谱分析。
分别打印其幅频特性，并进行分析和讨论。
四、实验代码如下所示：
clear;%清除变量
closeall;%关闭全部绘图窗口
b=menu（'请选择信号x1（n）--x6（n）','x1（n）=R4（n）','x2（n）=[12344321]',...
'x3（n）=[43211234]','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，则退出循环
ifb==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）;
ifb==1
x=[1,1,1,1,0,0,0,0,zeros（1,N-8）];
elseifb==2
x=[1,2,3,4,4,3,2,1,zeros（1,N-8）];
elseifb==3
x=[4,3,2,1,1,2,3,4,zeros（1,N-8）];
elseifb==4
x=cos（n*pi/4）;
elseifb==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）=[12344321]',...
'x3（n）=[43211234]','x4（n）=cos（npi/4）','x5（n）=sin（npi/8）',...
'x6（n）=cos（8pit）+cos（16pit）+cos（20pit）','Exit'）;
end
closeall;
五、实验结果图及分析
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=fs/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=fs/N=8Hz，因为截取的不是为周期序列的整数倍，而且频率分辨率不够，所得出的谱有较大的误差，所以FFT图包含一些频率分量，不能清楚看清原信号的频率f。
当N=16及32时，频率分辨率F0=fs/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=fs/N=4Hz，因为截取的不是x6里各周期序列的整数倍，所得出的谱有频谱泄漏，FFT图里可以看出信号cos（8pit）和cos（16pit）的频率f1=4,f2=8（在点N=1，2处有较大的幅值），而且频率分辨率不够高，不能分辨开第三个信号cos（20pit）的频率f3。
当N=32，频率分辨率F0=fs/N=2Hz，因此在FFT图里的点N=2有高幅值，在N=4有高幅值N=5也有高幅值。
因为截取的为周期序列的整数倍，所以所得出的谱正确。
当N=64，频率分辨率F0=fs/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：
x3（n）的波形N=8
图1-8：
x3（n）的波形N=16
图1-9：
x3（n）的波形N=32
图1-10：
x4（n）的波形N=8
图1-11：
x4（n）的波形N=16
图1-12：
x4（n）的波形N=32
图1-13：
x5（n）的波形N=8
图1-14：
x5（n）的波形N=16
图1-15：
x5（n）的波形N=32
图1-16：
x6（n）的波形N=16
图1-17：
x6（n）的波形N=32
图1-18：
x6（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、编写滤波器仿真程序，计算H（z）对心电信号采样序列x（n）的响应序列y（n）。
设yk（n）为第k级二阶滤波器Hk（z）的输出序列，yk－1（n）为输入序列，如图1所示。
由式2.2可以得到差分方程：
当k=1时，
。
所以H（z）对x（n）的总响应序列y（n）可以用顺序迭代算法得到。
即以此对k=1,2,3，求解差分方程式2.3，最后得到
。
仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。
也可以直接调用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）））;gridon;
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）））;gridon;
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））;
gridon;
ylabel（'|H|'）;xlabel（'w/π'）;
title（'幅度响应'）;
axis（[0,1,0,1.1]）;
subplot（3,2,2）,plot（w/pi,dbH）;
gridon;
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信号处理工具箱函数buttpbuttorbutter是巴特沃斯滤波器设计函数，其有5种调用格式，本实验中用到的是[N,wc]=butter（N,wc,Rp,As,’s’）,该格式用于计算巴特沃斯模拟滤波器的阶数N和3dB截止频率wc。
函数butter,cheby1和ellip设计IIR滤波器时都是默认的双线性变换法，所以在设计滤波器时只需要代入相应的实现函数即可。
1、设计流程图
语音信号滤波去噪—使用脉冲不变响应法设计的巴特沃斯滤波器的设计流程如图所示：
图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:
\USER\bin\Yuyin2.wav'）;%载入波形
N=length（y）%计算输入波形长度
Y=fft（y,N）;%进行快速傅里叶变换
%sound（y）;%播放声音
figure（5）;plot（abs（Y））;%滤波前频谱图
Fp=1200;%阻带截止频率
Fs=1100;%通带截止频率
Ft=8000;%采集频率
As=20;%通带波纹Ap=1
Ap=1;%阻带波纹As=20
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'）;%求S域的频率响应的参数
[num,den]=bilinear（b,a,1）;%双线性变换实现S域到Z域的变换
[h,w]=freqz（num,den）;%根据参数求出频率响应
figure（3）
plot（w*8000*0.5/pi,abs（h））;
z=filter（num,den,y）;
sound（z）;
m=z;%求滤波后的信号
figure
（1）
subplot（2,2,3）;
plot（abs（m）,'r'）;
title（'滤波前信号的频谱'）;
grid;
subplot（2,2,4）;
plot（z,'b'）;
title（'滤波前的信号波形'）;
grid;
subplot（2,2,2）;
plot（y,'b'）;
title（'滤波后信号的波形'）;
grid;
subplot（2,2,1）;
plot（abs（Y）,'r'）;
title（'滤波后信号的频谱'）;
grid;
figure
（2）;
p=angle（m）;
q=angle（Y）;
subplot（2,1,1）;plot（q,'b'）;title（'滤波前相位'）;
grid;
subplot（2,1,2）;plot（p,'b'）;title（'滤波后相位'）;
grid;
（2）实验结果:
所得波形如下图所示：
图3-7滤波前后的时域对比波形和滤波前后频域频谱图
图3-8滤波前后相位比较图
六、结果分析
数字滤波器的功能就是对输入离散信号的数字代码进行运算处理，以达到改变信号频谱的目的。
由图3-7中滤波前后波形比较可看出，经过滤波后的波形比原波形的振幅有所减小，去除了很多由于噪声所产生的干扰；从滤波前后的频谱比较可以看出经过滤波后除了原本的声音外，中间由于噪声产生的频谱波形已经滤除；由图3-8滤波前后相位比较图可看出由于经过滤波，相位变得稀疏且全都为正相位；经过MATLAB仿真，听滤波前后的声音，可以听出有明显的滤波效果。
因此利用脉冲响应不变法设计的巴特沃斯滤波器已经达到了设计的要求。
升级会员 