数字信号实验报告1.docx
《数字信号实验报告1.docx》由会员分享，可在线阅读，更多相关《数字信号实验报告1.docx（23页珍藏版）》请在冰豆网上搜索。
数字信号实验报告1
数字信号处理导论实验报告
姓名：
金涛
学号：
201005090209
实验一信号、系统及系统响应
当系统的输入输出差分方程为：
Y（n）-0.8y（n-1）-0.5y（n-2）=0.7x（n）+0.3x（n-1）
通过MATLAB编程实现并考虑如下问题：
（1）当系统的输入为单位冲激函数时，分别利用filter函数和impz函数得到的系统单位冲激响应曲线。
（2）当系统的输入为单位阶跃函数时，分别利用filter函数和impz函数得到的系统单位阶跃响应曲线。
（3）对于不同的输入，系统的输出有什么变化，试讨论之。
（1）实验原理
对一个给定的LSI系统，其转移函数H（z）的定义和表示形式为：
.习惯上，令H（z）=B（z）/A（z）.在MATLAB中，因为数组的下标不能为零（当然也不能为负值），因此，可重新表示为H（z）=
在有关MATLAB的系统分析的文件中。
分子和分母的系数被定义为向量，即
b=[b
（1）,b
（2）,b（3）,...,b（
+1）]
a=[a
（1）,a
（2）,...,a（
）]
并要求a
（1）=1，如果a
（1）≠1,则程序自动将其规划为1
（2）实验内容
源程序
x=ones（100）;t=1:
100;%产生单位阶跃序列
b=[.7,.3];%b向量
a=[1,.8,.5];%a向量
y=filter（b,a,x）;%实现实验1，图
（1）
plot（t,x,'g.',t,y,'k-'）;
[h,t]=impz（b,a,40）;%求出单位抽样响应。
图
（2）
stem（t,h,'.'）;gridon;
t=0:
20;
x=[1,zeros（1,20）];%产生单位采样序列
b=[.7,.3];%形成b向量
a=[1,.8,.5];%形成a向量
y=filter（b,a,x）;%filter函数图（3）
plot（t,x,'g.',t,y,'k-'）;
t=0:
20;
x=[1,zeros（1,20）];
b=[.7,.3];
a=[1,.8,.5];
[h,t]=impz（b,a,40）;%impz函数图（4）
stem（t,h,'.'）;gridon;
（3）实验结果
（1）
（2）
（3）
（4）
（4）分析结果
输入为离散是，输出为连续。
相反，输入为连续，输出为离散。
对于单位阶跃和单位抽样输入来说，输出没有变化
实验二使用FFT作频谱分析
（1）使用FFT对MATLAB中的数据noissin.dat进行谱分析。
（2）使用FFT对语音数据noisyspeech.wav进行谱分析。
（1）实验原理
（1）离散傅里叶变换（DFT）公式为：
X（k）=∑x（n）W^nk;
x（n）=∑X（k）W^-nk;
其中w=e^（2∏nk/N）,N为离散序列的长度。
（2）快速傅里叶变换（FFT）是利用w因子的取值特点，减少DFT的复数乘法的次数。
其中一种是时间抽取基2算法，它将时间按奇偶逐级分开，直到两点的DFT。
MATLAB提供了fft函数可用于计算FFT，器调用形式为;X=fft（x）或X=fft（x，N），N为2的整数次幂，若x的长度小于N，则补零，若超过则舍去N后的数据。
（3）函数形式[y,fs,bits]=wavread（'Blip',[N1N2]）;用于读取语音，采样值放在向量y中，fs表示采样频率（Hz），bits表示采样位数。
[N1N2]表示读取从N1点到N2点的值（若只有一个N的点则表示读取前N点的采样值）。
函数形式s=noissin（n1，n2）用于读取MALAB的噪声信号。
（2）实验内容
xx=wavread（'noisyspeech.wav'）;
fs=100;N=128;
x=xx（1:
N）;
n=1:
N;
X=abs（fft（x,N））;
subplot（221）;plot（n,x）;
xlabel（'n'）;ylabel（'x（n）'）;
gridon;
subplot（222）;plot（n,X）;
xlabel（'k'）;ylabel（'|X（k）|'）;
gridon;
loadnoissin;
s=noissin（1:
20）;
S=fft（s）;
subplot（223）;plot（abs（s））;
xlabel（'n'）;ylabel（'|s（n）|'）;
gridon;
subplot（224）;plot（abs（S））;
xlabel（'k'）;ylabel（'|S（k）|'）;
gridon;
（3）实验结果及分析
实验三使用双线性Z变换设计IIR滤波器
使用双线性Z变换法设计一个低通数字IIR滤波器，给定的数字滤波器的技术指标为fp=100Hz,fs=300Hz,ap=3dB,as=20dB,抽样频率Fs=1000Hz
（1）实验原理
1）设计滤波器就是要设计一个系统是其能让一定频率的波段通过或滤去，对IIR滤波器，器转移函数是：
H（Z）=（∑bz^（-r））/（1+∑az^（-k））。
（2）设计的一般原则：
若使滤波器拒绝某个频率，应在单位园上相应的频率处设置一个零点，反之则设置一个极点。
（3）低通数字IIR滤波器设计步骤：
给出数字低通滤波器的技术要求→映射为模拟低通的技术要求→归一化为模拟低通滤波器的技术要求→设计出G（P）→G（S）→映射到数字滤波器的转移函数G（Z）。
（3）双线性Z变换，即S平面到Z平面的映射关系：
S=2（Z-1）/Ts（Z+1）。
（2）实验内容
fp=100;fs=300;Fs=1000;%给定要设计的低通滤波器的设计频率；
rp=3;rs=20;%给定需要的衰减，单位为db；
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;Fs=Fs/Fs;%令Fs=1
wap=tan（wp/2）;was=tan（ws/2）;%实现由数字滤波器到模拟滤波器的转换；
[n,wn]=buttord（wap,was,rp,rs,'s'）;%求模拟滤波器的阶次；
[z,p,k]=buttap（n）;%设计模拟低通巴特沃斯滤波器。
得到极零点；
[bp,ap]=zp2tf（z,p,k）;%由极零点得到滤波器分子分母多项式的系数；
[bs,as]=lp2lp（bp,ap,wap）;%归一化低通到实际低通的转换；
[bz,az]=bilinear（bs,as,Fs/2）;%实现s到z的转换；
[h,w]=freqz（bz,az,256,Fs*1000）;%求出滤波器的频率响应；
plot（w,abs（h））;gridon;
（3）实验结果
（4）结果分析
由图可知，数字滤波器的幅频曲线完全满足技术要求，而且在f>fp后，数字滤波器的幅频响应的衰减比较快，这正是我们希望的
实验四窗函数法设计FIR数字滤波器
根据下列指标，设计一个FIR数字低通滤波器：
ω=2.0∏，ω=4.0∏
Ap=0.25dB，As=50dB。
（1）分别考查矩形窗、三角窗、Hanning窗、Hamming窗并分析这些不同的窗函数对滤波器的幅度响应有什么影响。
（2）选择一个合适的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。
（一）实验目的：
了解常用的几种窗函数，能正确选择适当的窗函数进行滤波器设计；掌握窗函数法设计数字低通滤波器。
（二）实验原理
常用的窗函数：
1、矩形窗函数为boxcar和rectwin，调用格式：
w=boxcar（N）w=rectwin（N）
其中N是窗函数的长度，返回值w是一个N阶的向量。
2、三角窗函数为triang，调用格式：
w=triang（N）
3、汉宁窗函数为hann，调用格式：
w=hann（N）
4、海明窗函数为hamming，调用格式：
w=hamming（N）
各个窗函数的性能比较
窗函数
第一瓣相对于主瓣衰减（dB）
主瓣宽
阻带最小衰减
矩形窗
-13
4∏/N
21
三角窗
-25
8∏/N
25
汉宁窗
-31
8∏/N
44
海明窗
-41
8∏/N
53
（三）实验内容
1、题一：
生成四种窗函数。
矩形窗、三角窗、汉宁窗、海明窗，并观察其频率响应。
源程序：
clearall;
n=20;
window=rectwin（n）;
[h,w]=freqz（window,1）;
subplot（4,2,1）;
stem（window）;
title（'矩形窗'）;
subplot（4,2,2）;
plot（w/pi,20*log（abs（h））/abs（h
（1）））;
title（'矩形窗的频率响应'）;
holdon;
window=triang（n）;
[h,w]=freqz（window,1）;
subplot（4,2,3）;
stem（window）;
title（'三角窗'）;
subplot（4,2,4）;
plot（w/pi,20*log（abs（h））/abs（h
（1）））;
title（'三角窗的频率响应'）;
holdon;
window=hann（n）;
[h,w]=freqz（window,1）;
subplot（4,2,5）;
stem（window）;
title（'汉宁窗'）;
subplot（4,2,6）;
plot（w/pi,20*log（abs（h））/abs（h
（1）））;
title（'汉宁窗的频率响应'）;
holdon;
window=hamming（n）;
[h,w]=freqz（window,1）;
subplot（4,2,7）;
stem（window）;
title（'海明窗'）;
subplot（4,2,8）;
plot（w/pi,20*log（abs（h））/abs（h
（1）））;
title（'海明窗的频率响应'）;
holdon;
实验结果：
分析实验结果：
通过对窗函数的调用，使4个窗函数的时域及频域图放置在一起，比较直观地看出各个窗函数的特性。
2、题二：
根据下列技术指标，设计一个FIR数字低通滤波器：
wp=0.2π，ws=0.4π，ap=0.25dB，as=50dB，选择一个适当的窗函数，确定单位冲激响应，绘出所设计的滤波器的幅度响应。
提示：
根据窗函数最小阻带衰减的特性表，可采用海明窗可提供大于50dB的衰减，其过渡带为6.6π/N，因此具有较小的阶次。
源程序：
clearall;
wp=0.2*pi;%通带截止频率
ws=0.4*pi;%阻带截止频率
tr_wdith=ws-wp;%过渡带宽度
N=ceil（6.6*pi/tr_wdith）+1;%得到N
n=0:
1:
N-1;
wc=（ws+wp）/2;%理想低通滤波器的截止频率
hd=ideal_lp1（wc,N）;%理想低通滤波器的单位冲激响应
w_ham=（hamming（N））';%海明窗
h=hd.*w_ham;%截取得到实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2（h,1）;%计算实际滤波器的幅度响应
subplot（221）;
stem（n,hd,'k'）;
title（'理想低通滤波器的单位脉冲响应hd（n）'）;
subplot（222）;
stem（n,w_ham,'k'）;
title（'海明窗w（n）'）;
subplot（223）;
stem（n,h,'k'）;
title（'实际低通滤波器的单位脉冲响应h（n）'）;
subplot（224）;
plot（w/pi,db,'k'）;
title（'实际低通滤波器的幅度响应（dB）'）;
axis（[0,1,-100,10]）;
子函数freqz_m2（b,a）：
function[db,mag,pha,w]=freqz_m2（b,a）;
%滤波器的幅值响应（相对、绝对）、相位响应
%db：
相对幅值响应；
%mag：
绝对幅值响应；
%pha：
相位响应；
%w：
采样频率；
%b:
系统函数H（Z）的分子项（对FIR,b=h）;
%a:
系统函数H（Z）的分母项（对FIR,a=1）;
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
1:
501））';
w=（w（1:
1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）
子函数ideal_lp1（wc,N）：
functionhd=ideal_lp1（wc,N）;
alpha=（N-1）/2;
n=0:
1:
N-1;
m=n-alpha+eps;
hd=sin（wc*m）./（pi*m）;
实验结果：
分析实验结果：
用窗函数的方法设计FIR低通数字滤波器，其性能十分逼近与理想的IIR数字低通滤波器，并且具有较好的衰减特性。
如图所示。
实验五用C++编程实现FFT并作频分析
（一）源程序
#include
#include
#include
constintN=1024;
constfloatPI=3.1416;
inlinevoidswap（float&a,float&b）
{
floatt;
t=a;
a=b;
b=t;
}
voidbitrp（floatxreal[],floatximag[],intn）
{
//位反转置换Bit-reversalPermutation
inti,j,a,b,p;
for（i=1,p=0;i{
p++;
}
for（i=0;i{
a=i;
b=0;
for（j=0;j{
b=（b<<1）+（a&1）;//b=b*2+a%2;
a>>=1;//a=a/2;
}
if（b>i）
{
swap（xreal[i],xreal[b]）;
swap（ximag[i],ximag[b]）;
}
}
}
voidFFT（floatxreal[],floatximag[],intn）
{
//快速傅立叶变换，将复数x变换后仍保存在x中，xreal,ximag分别是x的实部和虚部
floatwreal[N/2],wimag[N/2],treal,timag,ureal,uimag,arg;
intm,k,j,t,index1,index2;
bitrp（xreal,ximag,n）;
//计算1的前n/2个n次方根的共轭复数W'j=wreal[j]+i*wimag[j],j=0,1,...,n/2-1
arg=-2*PI/n;
treal=cos（arg）;
timag=sin（arg）;
wreal[0]=1.0;
wimag[0]=0.0;
for（j=1;j{
wreal[j]=wreal[j-1]*treal-wimag[j-1]*timag;
wimag[j]=wreal[j-1]*timag+wimag[j-1]*treal;
}
for（m=2;m<=n;m*=2）
{
for（k=0;k{
for（j=0;j{
index1=k+j;
index2=index1+m/2;
t=n*j/m;//旋转因子w的实部在wreal[]中的下标为t
treal=wreal[t]*xreal[index2]-wimag[t]*ximag[index2];
timag=wreal[t]*ximag[index2]+wimag[t]*xreal[index2];
ureal=xreal[index1];
uimag=ximag[index1];
xreal[index1]=ureal+treal;
ximag[index1]=uimag+timag;
xreal[index2]=ureal-treal;
ximag[index2]=uimag-timag;
}
}
}
}
voidIFFT（floatxreal[],floatximag[],intn）
{
//快速傅立叶逆变换
floatwreal[N/2],wimag[N/2],treal,timag,ureal,uimag,arg;
intm,k,j,t,index1,index2;
bitrp（xreal,ximag,n）;
//计算1的前n/2个n次方根Wj=wreal[j]+i*wimag[j],j=0,1,...,n/2-1
arg=2*PI/n;
treal=cos（arg）;
timag=sin（arg）;
wreal[0]=1.0;
wimag[0]=0.0;
for（j=1;j{
wreal[j]=wreal[j-1]*treal-wimag[j-1]*timag;
wimag[j]=wreal[j-1]*timag+wimag[j-1]*treal;
}
for（m=2;m<=n;m*=2）
{
for（k=0;k{
for（j=0;j{
index1=k+j;
index2=index1+m/2;
t=n*j/m;//旋转因子w的实部在wreal[]中的下标为t
treal=wreal[t]*xreal[index2]-wimag[t]*ximag[index2];
timag=wreal[t]*ximag[index2]+wimag[t]*xreal[index2];
ureal=xreal[index1];
uimag=ximag[index1];
xreal[index1]=ureal+treal;
ximag[index1]=uimag+timag;
xreal[index2]=ureal-treal;
ximag[index2]=uimag-timag;
}
}
}
for（j=0;j{
xreal[j]/=n;
ximag[j]/=n;
}
}
voidFFT_test（）
{
charinputfile[]="input.txt";//从文件input.txt中读入原始数据
charoutputfile[]="output.txt";//将结果输出到文件output.txt中
floatxreal[N]={},ximag[N]={};
intn,i;
FILE*input,*output;
if（!
（input=fopen（inputfile,"r"）））
{
printf（"Cannotopenfile."）;
exit
（1）;
}
if（!
（output=fopen（outputfile,"w"）））
{
printf（"Cannotopenfile."）;
exit
（1）;
}
i=0;
while（（fscanf（input,"%f%f",xreal+i,ximag+i））!
=EOF）
{
i++;
}
n=i;//要求n为2的整数幂
while（i>1）
{
if（i%2）
{
fprintf（output,"%disnotapowerof2!
",n）;
exit
（1）;
}
i/=2;
}
FFT（xreal,ximag,n）;
fprintf（output,"FFT:
irealimag"）;
for（i=0;i{
fprintf（output,"%4d%8.4f%8.4f",i,xreal[i],ximag[i]）;
}
fprintf（output,"================================="）;
IFFT（xreal,ximag,n）;
fprintf（output,"IFFT:
irealimag"）;
for（i=0;i{
fprintf（output,"%4d%8.4f%8.4f",i,xreal[i],ximag[i]）;
}
if（fclose（input））
{
printf（"Filecloseerror."）;
exit
（1）;
}
if（fclose（output））
{
printf（"Filecloseerror."）;
exit
（1）;
}
}
intmain（）
{
FFT_test（）;
return0;
}
（二）实验结果
/**************************************************
////////////////////////////////////////////////
//sample:
input.txt
////////////////////////////////////////////////
升级会员 