数字信号处理报告.docx
《数字信号处理报告.docx》由会员分享，可在线阅读，更多相关《数字信号处理报告.docx（26页珍藏版）》请在冰点文库上搜索。
数字信号处理报告
《数字信号处理》
课
程
设
计
报
告
设计一信号的基本运算
【一】设计目的
熟悉信号的基本运算，通过仿真结果加深理解。
【二】设计内容
1.现有一段以采样频率8192Hz采样得到的语音数据x（k），为了模仿出回音的效果，可以对数据进行如下处理：
x（k）=x（k）+a*x（k-d），其中d为时延，a为时延信号的衰减幅度。
根据上述原理，进行仿真实验，并演示回声的效果。
（参数：
时延d=0.4秒，衰减幅度a=0.5）
2.对上述语音信号进行时间反转x（-k）、上采样x（2k）、下采样x（k/2）操作，并演示这些运算的效果。
【三】设计过程，调试，分析。
数字信号处理的主要研究对象是语音信号和图像信号,语音信号的研究可以从时域和频域两个方面来进行。
其中时域的分析处理有两种方法:
一种是进行语音信号分析,这属于线性处理的范畴,主要是通过信号的加减、时移、倍乘、卷积、求相关函数等来实现;另一种是生成和变换成各种调制信号,这属于非线性的范畴,主要是对信号平均累加器的动态范围进行压缩扩张,用门限方法对噪声的抑制。
对频域分析处理,即对信号的频率特性在频谱中加以分析研究,这拓展了信号分析的范围,是对不确定信号分析的主要方法。
在实际应用中,信号的时域频分析经常同时进行。
1语音信号的采集与时频分析
1.1语音信号的采集
利用PC机上的声卡和WINDOWS操作系统可以进行数字信号的采集。
将话筒输入计算机的语音输入插口上,启动录音机。
按下录音按钮,接着对话筒说话/大家好0,说完后停止录音,屏幕左侧将显示所录声音的长度。
点击放音按钮,可以实现所录音的重现。
以文件名/dajiahao0保存入g:
\MATLAB\work中。
可以看到,文件存储器的后缀默认wav。
在Matlab软件平台下可以利用函数wavread对语音信号进行采样,得到了声音数据变量x,同时把x的采样频率Fs=8000Hz和数据位Nbits=8Bit放进了MATALB的工作空间。
其程序如下
Fs=8000Hz;%给出抽样频率
Nbits=8Bit;%放入数据位
[x,fs,Nbits]=wavread（-g:
\Matlab\work\chaotian.wav）;%把语音信号进行加载入Matlab仿真软件平台中
sound（x,fs）;%对加载的语音信号进行回放
stem（x）;title（‘语音信号的时域波形’）;
1.2语音信号的时域分析
因录音采样频率无法准确。
故用余弦序列表示信号的迟延，反转，上采样下采样。
迟延：
题目要求时延d为0.4s.取点数k为0.4*Fs。
产生一组新序列将原序列依次推后k点。
【四】程序清单与运行结果
程序清单：
y=load（'E:
\erf1s1t0'）;
fs=8192;
n=length（y）;
m=n+round（0.4*fs）;
q=（m-1）/fs;
l=round（0.4*fs）;
z=zeros（1,m）;
x=z;
fori=1:
nx（i）=y（i）;end%原信号
fori=1:
nz（l+i）=y（i）;end%信号的延拖
o=x+0.4*z;%回声信号
sound（o）%声音
t=0:
1/fs:
q;
subplot（3,1,1）;
plot（t,x）
title（'原信号'）;
subplot（3,1,2）;
plot（t,z）
title（'时延信号'）;
subplot（3,1,3）;
plot（t,o）
title（'回声信号'）;
运行结果：
翻转：
K为序列长度，将序列第i个点值依次与第k-i点值进行调换产生新序列，即为翻转序列。
采样：
下采样，也就是对信号的抽取。
其实，上采样和下采样都是对数字信号进行重采，重采的采样率与原来获得该数字信号的采样率比较，大于原信号的称为上采样，小于的则称为下采样。
上采样的实质也就是内插或插值。
程序清单：
%对上述语音信号进行时间反转x（-k）、上采样x（2k）、下采样x（k/2）操作
data=load（'E:
\erf1s1t0'）;
x=data;%原声信号
d=0.4;%d为时延
Fs=8192;%Fs为采样频率
k=0.4*8192;
a=0.5;%a为时延信号的衰减幅度
z=linspace（0,0,k）;%linspace为线性间隔的向量
y=[zx'];
x=[x'z];
f=x+a*y;%回声信号
f1=fliplr（x）;%翻转信号
f2=downsample（x,2）;%下采样信号
f3=upsample（x,2）;%上采样信号
subplot（4,1,1）;
plot（x）;title（'原声信号'）
subplot（4,1,2）;
plot（f1）;title（'翻转信号'）
subplot（4,1,3）;
plot（f2）;title（'下采样'）
subplot（4,1,4）;
plot（f3）;title（'上采样'）
sound（f1）
运行结果：
设计二正余弦信号的谱分析
【一】设计目的
用DFT实现对正余弦信号的谱分析；
观察DFT长度和窗函数长度对频谱的影响；
对DFT进行谱分析中的误差现象获得感性认识。
【二】设计内容
1.对一个频率为10Hz，采样频率为64Hz的32点余弦序列进行谱分析，画出其频谱图；若将频率改为11Hz，其他参数不变，重新画出该序列的频谱图，观察频率泄漏现象，分析原因。
2.考察DFT的长度对双频率信号频谱分析的影响。
设待分析的信号为
（1.2）
令两个长度为16的正余弦序列的数字频率为
及
。
取N为四个不同值16，32，64，128。
画出四个DFT幅频图，分析DFT长度对频谱分辨率的影响。
3.在上题中若把两个正弦波的频率取得较近，令
，
，试问怎样选择FFT参数才能在频谱分析中分辨出这两个分量？
【三】设计过程，调试，分析。
连续信号的频谱分析在工程上有着广泛的应用，计算机的第一步就是把连续信号离散化，一是采样，二是截断。
由此会产生频率混叠误差、栅栏效应和截断误差。
（一）对一个频率为10Hz，和11Hz，采样频率为64Hz的32点余弦序列用Matlab计算它的DFT程序为
N=32;Fs=64;f=10;
n=0:
N-1;
x=cos（2*pi*10*n/64）
y=cos（2*pi*11*n/64）
subplot（2,1,1）
X=abs（fft（x,N））
stem（n,X,'.'）;
xlabel（'n'）;title（'10HZ32点变化幅频曲线'）
subplot（2,1,2）
Y=abs（fft（y,N））
stem（n,Y,'.'）;
xlabel（'n'）;title（'11HZ32点变化幅频曲线'）
运行结果：
通过图可看出：
频率为10Hz的余弦曲线DFT只有两个点不等于零，位于k=5和k=27处，k=5对应于频率10Hz，k=27对应于频率54Hz（也就是-10Hz）。
这样DFT确实正确的分辨了余弦信号的频率。
但是这样的理想结果是恰巧得到的，此时我们借去了五个完整的余弦周期（f*N/Fs=5）.
将频率改为11Hz，采样频率和窗长度依然为32点，计算图像可看出：
频谱图上k=5和k=27处都有较大的峰值，而其它的点上幅度不再为零。
这两个峰值对应的频率为10Hz和12Hz，所以，信号的峰值位于两者之间，本来是单一的11Hz频率的能量会分不到许多DFT频率上，这种现象叫频率泄露，来源于截断效应。
（二）对双信号的频谱分析程序为：
16点，32点，64点，128点的频谱分析程序：
N=16;n=0:
15;
f1=0.22;f2=0.34;
x=0.5*sin（2*pi*f1*n）+sin（2*pi*f2*n）;
subplot（4,2,1）,stem（n,x,'.'）;
%xlable（'n'）;ylable（'x1（n）'）;
title（'余弦序列'）;
X=abs（fft（x,N））;
subplot（4,2,2）
k=0:
N-1;
stem（k,X,'.'）
%xlable（'k'）;ylable（'X（k）'）
string=[num2str（N）,'点FFT幅频曲线'];
title（string）;
N=32;n=0:
15;
f1=0.22;f2=0.34;
x=0.5*sin（2*pi*f1*n）+sin（2*pi*f2*n）;
subplot（4,2,3）,stem（n,x,'.'）
%xlable（'n'）;ylable（'x1（n）'）;
title（'余弦序列'）;
X=abs（fft（x,N））;
subplot（4,2,4）
k=0:
N-1;
stem（k,X,'.'）
%xlable（'k'）;ylable（'X（k）'）
string=[num2str（N）,'点FFT幅频曲线'];
title（string）;
N=64;n=0:
15;
f1=0.22;f2=0.34;
x=0.5*sin（2*pi*f1*n）+sin（2*pi*f2*n）;
subplot（4,2,5）,stem（n,x,'.'）
%xlable（'n'）;ylable（'x1（n）'）;
title（'余弦序列'）;
X=abs（fft（x,N））;
subplot（4,2,6）
k=0:
N-1;
stem（k,X,'.'）
%xlable（'k'）;ylable（'X（k）'）
string=[num2str（64）,'点FFT幅频曲线'];
title（string）;
N=128;n=0:
15;
f1=0.22;f2=0.34;
x=0.5*sin（2*pi*f1*n）+sin（2*pi*f2*n）;
subplot（4,2,7）,stem（n,x,'.'）
%xlable（'n'）;ylable（'x1（n）'）;
title（'余弦序列'）;
X=abs（fft（x,N））;
subplot（4,2,8）
k=0:
N-1;
stem（k,X,'.'）
%xlable（'k'）;ylable（'X（k）'）
string=[num2str（128）,'点FFT幅频曲线'];
title（string）;
运行结果:
分析：
DFT样本值就是其DTFT在相应位置的采样。
在图一中很难看出两个峰值，因此要提高它的分辨率，故把N增大，逐渐可以看出它有两个峰值，将k换算成数字频率f=w/2*pi=k/N.这样可确定峰值的位置大体在f=0.21和0.35之附近，与信号的给定频率有一定的误差，这也是截断和泄露带来的问题，在这图上还可以看到一些较小的峰，这是很难判断是输入信号固有的，还是由泄露引起的。
说明了增加DFT长度N减小了相邻样本间的频率间距，提高频谱的视在分辨率，因而可以提高样本未知的测定精度。
（三）把两个正弦波的频率取得较近，另fr=[0.22,0.25]应怎样选择FFT参数
要能分清两个频率，分辨率至少应达到f=0.03.因为此处的数字频率是对采样频率Fs进行归一化后的，几fr最大为1，因此总的样本数至少要达到1/0.03=33
加窗以后可以使频谱函数更加光滑，便于分辨峰值位置和准确的数值，为了提高实际分辨率，应该尽量增加信号的长度n及DF长度N,当受到条件限制不能提高n，则单独提高N可以提高视在分辨率。
当
，
时，16点，32点，64点，128点的频谱运行结果：
（四）谱分析中的误差现象
①混叠误差。
产生混叠误差的原因是：
信号的离散化是通过抽样实现的，而抽样频率再高总是有限的。
除带限信号外，如果信号的最高频率Ω趋于无穷，则实际器件无法满足抽样定理。
而抽样过程如果不满足抽样定理，就会产生频谱的混叠，即混叠误差。
要减少或避免混叠误差，应该提高抽样频率，以设法满足抽样定理，或者采用抗混叠滤波这样的信号预处理措施。
②栅栏效应。
对于非周期信号来说，理论上应具有连续的频谱，但数字谱分析是用DFT来近似的，是用频谱的抽样值逼近连续频谱值，只能观察到有限个频谱值，每一个间隔中的频谱都观察不到。
如同通过“栅栏”观察景物一样，一部分被“栅栏”所阻挡，看不见，把这种现象称为栅栏效应。
连续时间信号只要采用数字谱分析的方法，就必定产生栅栏效应，栅栏效应只能减小而无法避免。
能够感受的频谱最小间隔值，称为频谱分辨率，而频率分辨率与信号截断长度成反比，且频率分辨率越小，插值效果越好。
因此即使连续信号是有限时间的，也应选择一个大的截断长度来改善栅栏效应。
③截断误差。
截断误差是由于对信号进行截断，把无限长的信号限定为有限长，即令有限区间以外的函数值均为零值的近似处理而产生的，这种处理相当于用一个矩形（窗）信号乘待分析的连续时间信号。
必然会引起吉布斯效应（波动），也会把窗函数的频谱引入信号频谱,造成混叠，所以需要考虑其误差的问题。
设计三数字滤波器的设计及实现
【一】设计目的
1.熟悉IIR数字滤波器和FIR数字滤波器的设计原理和方法；
2.学会调用MATLAB信号处理工具箱中的滤波器设计函数设计各种IIR和FIR数字滤波器，学会根据滤波要求确定滤波器指标参数；
3.掌握用IIR和FIR数字滤波器的MATLAB实现方法，并能绘制滤波器的幅频特性、相频特性；
4.通过观察滤波器的输入、输出信号的时域波形及其频谱，建立数字滤波的概念。
【二】设计内容
1.调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，观察st的时域波形和幅频特性曲线；
2.通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率，假定要求滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB；
3.编程调用MATLAB滤波器设计函数分别设计这三个数字滤波器，并绘图显示其幅频特性曲线；
4.用所设计的三个滤波器分别对复合信号st进行滤波，分离出st中的三路不同载波频率的调幅信号，并绘图显示滤波后信号的时域波形和频谱，观察分离效果。
【三】设计过程，调试，分析。
（一）IIR数字滤波设计方法
数字滤波是数字信号处理的一种重要算法，广泛用于对信号的过滤、检测与参数的估计等信号处理中。
数字滤波器按其单位脉冲响应的长度可分为有限脉冲响应（FIR）滤波器和无限脉冲响应（IIR）滤波器两类。
数字滤波器按其通频带分为低通滤波器，高通滤波器，带通滤波器和带阻滤波器。
在本次实验中旨在用椭圆滤波器实现IIR低通，带通，高通滤波。
IIR滤波器基本概念
一个N阶递归型数字滤波器的差分方程
H（Z）=
MEbrz-r
r=0
N
1+K=Eakz-k
从以上的系统函数可知,设计IIR滤波器的任务就是通过计算寻求一个因果、物理上可实现的系统函数H（Z）,使其频率响应H（ejw）满足所希望得到的频域指标,即符合给定的通带截止频率、阻带截止、通带衰减和阻带衰减。
数字滤波器的设计过程
数字滤波器的一般设计过程可分为以下四个步骤:
1）按照实际的需要,确定滤波器的性能要求或指标;
2）用一个因果稳定的离散线性时不变系统,去逼近这一性能指标;
3）用有限精度的运算实现所设计的系统;
4）通过模拟,验证所设计的系统是否符合性能要求
（二）
1：
调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st,画出st的时域波形和幅频特性，三路信号在时域混叠无法在时域分离，但在频域分离，所以可以通过滤波的方法在频域分离。
抑制载波单频调幅信号的数学表达式为
这三路调幅信号的载波频率分别为250Hz、500Hz和1000Hz
2：
通过观察st的幅频特性，分别可以确定三路信号的低通，带通，高通滤波器。
的通带截止频率和阻带截止频率，要求滤波器的通带最大衰减为0.1db,阻带最小衰减为60db
分离混合信号st的滤波器指标参数选取如下：
载波为250Hz的调幅信号，用低通滤波器分离。
指标为：
通带截止频率fp=300HZ,阻带截止频率fs=400Hz,通带最大衰减ap=0.1db,阻带最小衰减as=60db
对载波频率为500Hz的调幅信号，用带通滤波器分离，指标为：
通带截止频率fl=400HZ,fu=500Hz,阻带截止频率fl=600Hz,fu=700Hz,通带最大衰减ap=0.1db,阻带最小衰减as=60db
对载波频率为1000Hz的调幅信号，用高通滤波器分离，其指标为：
通带截止频率fp=900HZ,阻带截止频率fs=800Hz,通带最大衰减ap=0.1db,阻带最小衰减，60db
3：
用椭圆滤波器实现
椭圆滤波器特点：
幅值响应在通带和阻带内都是等波纹的。
对于给定的阶数和给定的波纹要求，椭圆滤波器能获得较其它滤波器更窄的过渡带宽，就这点而言，椭圆滤波器是最优的。
通带和阻带内波纹固定时，阶数越高，过渡带越窄
阶数固定，通带和阻带波纹越小，过渡带越宽
【四】程序清单与运行结果
程序清单：
%1.调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，观察%st的时域
%波形和幅频特性曲线；
functionst=mstg
%产生信号序列st，并显示st的时域波形和频谱
%st=mstg返回三路调幅信号相加形成的混合信号，长度N=800
N=800;%信号长度N=800
Fs=10000;%采样频率Fs=10kHz,Tp为采样时间
T=1/Fs;
Tp=N*T;
t=0:
T:
（N-1）*T;
k=0:
N-1;
f=k/Tp;
fc1=Fs/10;%第一路调幅信号载波频率fc1=1000Hz
fm1=fc1/10;%第一路调幅信号的调制信号频率fm1=100Hz
fc2=Fs/20;%第二路调幅信号载波频率fc2=500Hz
fm2=fc2/10;%第二路调幅信号的调制信号频率fm2=50Hz
fc3=Fs/40;%第三路调幅信号载波频率fc3=250Hz
fm3=fc3/10;%第三路调幅信号的调制信号频率fm3=25Hz
xt1=cos（2*pi*fm1*t）.*cos（2*pi*fc1*t）;%产生第一路调幅信号
xt2=cos（2*pi*fm2*t）.*cos（2*pi*fc2*t）;%产生第二路调幅信号
xt3=cos（2*pi*fm3*t）.*cos（2*pi*fc3*t）;%产生第三路调幅信号
st=xt1+xt2+xt3;%三路信号相加，得到复合信号
fxt=fft（st,N）;%计算信号st的频谱
%以下为绘图命令
figure
（1）
subplot（2,1,1）
plot（t,st）;
grid;
xlabel（'t/s'）;
ylabel（'s（t）'）;
axis（[0,Tp,min（st）,max（st）]）;
title（'（a）s（t）的波形'）
subplot（2,1,2）
stem（f,abs（fxt）/max（abs（fxt））,'.'）;
grid;
title（'（b）s（t）的频谱'）
axis（[0,Fs/8,0,1.2]）;
xlabel（'f/Hz'）;
ylabel（'幅度'）
%低通滤波器
fp=320;fs=350;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;%DF指标（低通滤波器的通、阻带边界频）
[N,wp]=ellipord（wp,ws,rp,rs）;%调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip（N,rp,rs,wp）;%调用ellip计算椭圆带通DF系统函数系数向量B和A
y1t=filter（B,A,st）;%滤波器软件实现
y1=fft（y1t）;
%绘图部分
figure
（2）;
subplot（3,3,1）;
[H,W]=freqz（B,A,800）;
plot（W*Fs/2/pi,abs（H））;
xlabel（'Hz'）;
ylabel（'H（w）'）;
title（'低通滤波器'）;
axis（[0,2000,0,1.2]）;
grid;
subplot（3,3,2）;
plot（t,y1t）;
xlabel（'t'）;
ylabel（'y（t）'）;
title（'分离出的250Hz的波形'）;
subplot（3,3,3）;
stem（f,abs（y1）/max（abs（y1））,'.'）;
xlabel（'Hz'）;
ylabel（'|H|'）;
title（'频谱'）;
axis（[0,1200,0,1]）;
grid;
%带通滤波器
fpl=420;fpu=580;fsl=300;fsu=700;
wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;
[N,wp]=ellipord（wp,ws,rp,rs）;%调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip（N,rp,rs,wp）;%调用ellip计算椭圆带通DF系统函数系数向量B和A
y2t=filter（B,A,st）;%滤波器软件实现
y2=fft（y2t）;
%绘图部分
figure
（2）;subplot（3,3,4）;
[H,W]=freqz（B,A,800）;
plot（W*Fs/2/pi,abs（H））;
xlabel（'Hz'）;
ylabel（'H（w）'）;
title（'带通滤波器'）;
axis（[0,2000,0,1.2]）;
grid;
subplot（3,3,5）;
plot（t,y2t）;
xlabel（'t'）;
ylabel（'y（t）'）;
title（'分离出的500Hz的波形'）;
subplot（3,3,6）;
stem（f,abs（y2）/max（abs（y2））,'.'）;
xlabel（'Hz'）;
ylabel（'|H|'）;
title（'频谱'）;
axis（[0,1200,0,1]）;
grid;
%高通滤波器
fp=870;fs=780;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;%DF指标（低通滤波器的通、阻带边界频）
[N,wp]=ellipord（wp,ws,rp,rs）;%调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip（N,rp,rs,wp,'high'）;%调用ellip计算椭圆带通DF系统函数系数向量B和A
y3t=filter（B,A,st）;%滤波器软件实现
y3=fft（y3t）;
%绘图部分
figure
（2）;subplot（3,3,7）;
[H,W]=freqz（B,A,800）;
plot（W*Fs/2/pi,abs（H））;
xlabel（'Hz'）;
ylabel（'H（w）'）;
title（'高通滤波器'）;
axis（[0,2000,0,1.2]）;
grid;
升级会员 