数字信号处理实验1.docx
- 文档编号:23284948
- 上传时间:2023-05-15
- 格式:DOCX
- 页数:49
- 大小:370.24KB
数字信号处理实验1.docx
《数字信号处理实验1.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验1.docx(49页珍藏版)》请在冰豆网上搜索。
数字信号处理实验1
实验一离散系统时域和z域分析
一、实验目的
掌握两个有限长度序列的卷积和计算;
掌握差分方程解法;
学习和掌握离散系统频率特性。
二、实验原理
1.已知
直接计算卷积和
利用卷积定理求卷集和
若
,
则
2.LTI离散系统也能用下列形式的线性常系数差分方程来描述:
我们知道线性常系数差分方程有两种解的形式。
一种形式包括通解和特解,另一种形式包括零输入响应和零状态响应。
使用Z变换,提供了一种得到这两种形式的方法。
另外,我们也将讨论暂态响应和稳态响应。
在数字信号处理中,差分方程通常按正n方向给出。
因此,解的时间范围是n≥0为此,我们定义双边Z变换的一个变形,叫做单边Z变换。
已知差分方程,我们可以通过编程求出它对输入
的响应。
单边Z变换定义
3.系统函数定义
频响函数
定义
与频响函数
的情况类似,
可以定义Z域函数H(z),称为系统函数。
但是,与
不同的是,对于一个不稳定的系统,系统函数H(z)仍然存在。
给定一因果系统,通过使用函数pzplotz、freqz进行编程,可以画出零极点示意图
和
,即系统的幅频响应和相位响应。
三、实验内容
1、已知两个有限长度序列直接计算和利用卷积定理计算卷积和;
求解差分方程;
2、已知一个由差分方程描述的因果系统:
y(n)=a1*y(n-2)+b1*x(n)+b2*x(n-2)
求单位脉冲响应;
求频率响应函数并画出在0≤ω≤π区间的幅频和相频曲线;
求系统函数并画出零极点图。
四、实验说明
改进卷积函数:
function[y,ny]=conv_m(x,nx,h,nh)
[y,ny]=卷积结果
[x,nx]=第一个信号
[h,nh]=第二个信号
五、实验报告要求
1.阅读程序并运行程序,加上详细的注释;
2.运行所给例题给出机算的结果,并和笔算结果进行比较;
3.自选一道题编程求出结果并进行分析。
六、附程
信号处理的改进卷积函数:
function[y,ny]=conv_m(x,nx,h,nh)
nyb=nx
(1)+nh
(1);nye=nx(length(x))+nh(length(h));
ny=[nyb:
nye];
y=conv(x,h);
例1:
已知x1=[1,2,3],n1=[-1:
1],x2=[2,4,3,5],n2=[-2:
1],求x1(n)*x2(n)
程序如下:
x1=[1,2,3];n1=[-1:
1];
subplot(3,1,1);stem(n1,x1)
title('x1(n)')
xlabel('n');ylabel('x1(n)')
axis([-5,5,0,10])
set(gca,'XTickMode','manual','XTick',[-5:
1:
5])
x2=[2,4,3,5];n2=[-2:
1];
subplot(3,1,2);stem(n2,x2)
title('x2(n)')
xlabel('n');ylabel('x2(n)')
axis([-5,5,0,10])
set(gca,'XTickMode','manual','XTick',[-5:
1:
5])
[x3,n3]=conv_m(x1,n1,x2,n2)
subplot(3,1,3);stem(n3,x3)
title('x1(n)*x2(n)')
xlabel('n');ylabel('x3(n)')
axis([-5,5,0,30])
set(gca,'XTickMode','manual','XTick',[-5:
1:
5])
运行结果:
x3=
2817231915
n3=
-3-2-1012
例2:
求解差分方程:
n≥0
其中
程序如下:
s='求解差分方程y(n)=c1*x(n)+c2*x(n-1)+c3*x(n-2)+d1*y(n-1)+d2*y(n-2)';disp(s);
d='其中x(n)=cos(π*n/a)';disp(d);
b=[111]/3;
Y=[-2-3];
a=[1-0.950.9025];
X=[11];
%b=input('请输入序列的系数矢量x(n)=');Y=input('及其初始值[y(-1),y(-2)]=');
%a=input('请输入序列的系数矢量y(n)=');X=input('及其初始值[x(-1),x(-2)]=');
xic=filtic(b,a,Y,X)
bxplus=[1,-0.5];axplus=[1,-1,1];
ayplus=conv(a,axplus)
byplus=conv(b,bxplus)+conv(xic,axplus)
[R,p,C]=residuez(byplus,ayplus)
Mp=abs(p)
Ap=angle(p)/pi
n=[0:
50];
x=cos(pi*n/3);
y=filter(b,a,x,xic);
stem(n,y);
grid
运行结果:
初始条件转化结果:
xic=1.47422.1383
Yplus(z)的分母:
ayplus=1.0000-1.95002.8525-1.85250.9025
Yplus(z)的分子:
byplus=1.80750.8308-0.49751.9717
部分分式展开系数:
R=0.0584-3.9468i0.0584+3.9468i0.8453+2.0311i0.8453-2.0311i
p=0.5000+0.8660i0.5000-0.8660i0.4750+0.8227i0.4750-0.8227i
C=[]
部分分式展开系数极坐标形式:
Mp=1.00001.00000.95000.9500
Ap=0.3333-0.33330.3333-0.3333
则:
例3:
已知一个由差分方程描述的因果系统:
单位脉冲响应、阶跃响应;
求频率响应函数并画出在0≤ω≤π区间的幅频和相频曲线
求系统函数并画出零极点图
程序如下:
b=[10-1];
a=[10-0.81];
figure
(1)
subplot(2,1,1)
dimpulse(b,a,50)
gtext('impulseresponse')
subplot(2,1,2)
dstep(b,a,50)
gtext('stepresponse')
[R,p,C]=residuez(b,a)
pause
figure
(2)
w=[0:
1:
500]*pi/500;
freqz(b,a,w);
pause
figure(3)
zplane(b,a)
运行结果:
R=-0.1173-0.1173
p=-0.90000.9000
C=1.2346
得到:
实验二离散傅里叶变换和快速傅里叶变换
一、实验目的
学习和掌握DFT算法;
掌握时间抽取基-2FFT算法,了解优点。
二、实验原理
1.已知长度外为N的序列
则其离散傅立叶变换为
2.时间抽取基-2FFT算法原理
设序列长度N是2的整数幂次方,N=2M。
这样我们首先将序列分成两组:
偶数项一组,奇数项一组,得到2个N/2的子序列,即
,
对两个子序列分别进行DFT,可得
对于原序列
有
,
由
的性质可得
,
由于
仍然是偶数,可类似地将N/2点序列再分解为两个N/4点的序列,……,继续这样的分解过程,直到最后是2点序列的DFT。
三、实验内容
已知一个有限长度为N的序列
1、直接计算其离散傅立叶变换,并画出图形;
2、用FFT方法和IFFT方法计算其离散傅立叶变换和逆变换,并画出图形;
3、比较直接法和FFT方法的结果,研究二种方法的运算速度
4、比较添零和增加取样点,对频谱的影响。
四、实验说明
1.计算离散傅立叶变换函数:
function[Xk]=dft(xn,N)
Xk=在0≤n≤N-1间的DFT系数数组
Xn=N点有限持续时间序列
N=DFT长度
五、实验报告要求
阅读程序,并且加上详细的注释。
给出机算结果并针对各种结果进行分析。
六、附程
1.计算离散傅立叶变换函数:
function[Xk]=dft(xn,N)
n=[0:
1:
N-1];
k=[0:
1:
N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
例1:
已知一个有限长度为N的序列
,N=1000
直接计算和用时间抽取基-2FFT方法计算其离散傅立叶变换并画出图形,并比较计算所需时间。
直接计算离散傅立叶变换程序:
t0=clock;
figure
(1)
N=1000;
n=[0:
1:
999];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(2,1,1);plot(n,x);title('signalx(n),0<=n<=999');xlabel('n')
axis([0,1000,-2.5,2.5])
Xk=dft(x,N);magX=abs(Xk);phaX=angle(Xk);
k=0:
length(magX)-1;
subplot(2,1,2);plot(k,magX);title('DTFTMagnitude');
xlabel('frequencyinpiunits')
axis([0,1000,0,1000])
tc=etime(clock,t0)
时间单位(秒)
运行结果:
tc=
9.6100
时间抽取基-2FFT方法计算离散傅立叶变换程序:
t0=clock;
figure
(1)
n=[0:
1:
999];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(2,1,1);plot(n,x);title('signalx(n),0<=n<=999');xlabel('n')
axis([0,1000,-2.5,2.5])
X=fft(x);magX=abs(X);
k=0:
1:
999;
subplot(2,1,2);plot(k,magX);title('DTFTMagnitude');
xlabel('frequencyinpiunits')
axis([0,1000,0,1000])
tf=etime(clock,t0)
运行结果:
tf=
0.1100
例2:
已知信号
求DFT,对其结果进行IDFT,并将IDFT的结果和原信号进行比较。
程序如下:
clf
fs=100;
N=128;
n=0:
N-1;
t=n/fs;
x=sin(2*pi*40*t)+sin(2*pi*15*t);
subplot(2,2,1);plot(t,x);
xlabel('t(s)');ylabel('x(t)')
title('Originalsignal');grid
y=fft(x,N);
mag=abs(y);
f=(0:
length(y)-1)'*fs/length(y);
subplot(2,2,2);plot(f(1:
N/2),mag(1:
N/2));
xlabel('Frequency(Hz)');ylabel('magnitude');
title('FFTtooriginalsignal');grid
xifft=ifft(y);
magx=real(xifft);
ti=[0:
length(xifft)-1]/fs;
subplot(2,2,3);plot(ti,magx);
xlabel('t(s)');ylabel('x(t)');
title('SignalfromIFFT');grid
yif=fft(xifft,N);
mag=abs(yif);
f=(0:
length(y)-1)'*fs/length(y);
subplot(2,2,4);plot(f(1:
N/2),mag(1:
N/2));
xlabel('Frequency(Hz)');ylabel('Magnitude');
title('FFTtosignalfromIFFT');grid
运行结果:
例3:
已知
1.取x(n)(0≤n≤10)时,求x(n)的DFTX(k)
2.将
(1)中的x(n)以补零方式使x(n)加长到0≤n≤100,求X(k)
3.取x(n)(0≤n≤100)时,求X(k)
程序如下:
figure
(1)
n=[0:
1:
9];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(2,1,1);stem(n,x);title('signalx(n),0<=n<=9');xlabel('n')
axis([0,10,-2.5,2.5])
Xk=fft(x);magX=abs(Xk(1:
1:
6));
k=0:
1:
5;w1=2*pi/10*k;
subplot(2,1,2);stem(w1/pi,magX);title('DTFTMagnitude');
xlabel('frequencyinpiunits')
axis([0,1,0,10])
pause
figure
(2)
n=[0:
1:
99];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
y=[x(1:
1:
10),zeros(1,90)];
subplot(2,1,1);stem(n,y);title('signalx(n),0<=n<=9+90zeros');xlabel('n')
axis([0,100,-2.5,2.5])
Xk=fft(y);magX=abs(Xk(1:
1:
51));
k=0:
1:
50;w1=2*pi/100*k;
subplot(2,1,2);plot(w1/pi,magX);title('DTFTMagnitude');
xlabel('frequencyinpiunits')
axis([0,1,0,10])
pause
figure(3)
n=[0:
1:
99];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
subplot(2,1,1);stem(n,x);title('signalx(n),0<=n<=99');xlabel('n')
axis([0,100,-2.5,2.5])
Xk=fft(x);magX=abs(Xk(1:
1:
51));
k=0:
1:
50;w1=2*pi/100*k;
subplot(2,1,2);plot(w1/pi,magX);title('DTFTMagnitude');
xlabel('frequencyinpiunits')
axis([0,1,0,60])
运行结果:
实验三IIR数字滤波器实现
一、实验目的
1.学习和掌握系统传递函数映射的两种方法—冲击响应不变法和双线性变换法;
2.掌握和实现IIR数字滤波器的经典设计法。
二、实验原理
1.冲击响应不变法
就是使数字滤波器的冲击响应序列
等于模拟滤波器的冲击响应
的采样值,即
式中,T为采样周期。
已知
求
的方法是:
求模拟滤波器的冲击响应
:
求数字滤波器的冲击响应序列
求数字滤波器的系统函数:
2.双线性变换法
由于s平面和z平面的单值双线性映射关系为:
;
式中,T为采样周期。
则已知
求
的方法是:
3.IIR数字滤波器的经典设计法
根据给定的性能指标和方法不同,首先对设计性能指标中的频率指标,如边界频率进行转换,转换后的频率指标作为模拟滤波器原型设计的性能指标;
设计模拟低通原型滤波器;
由冲击响应不变法或双线性变换法获得IIR数字滤波器。
三、实验内容
1.编程实现由冲击响应不变法和双线性变换法将模拟滤波器变换为数字滤波器;
2.已知要设计IIR数字滤波器的性能指标采用经典设计法设计。
用冲击响应不变法和双线性变换法设计巴特沃思原型和切贝雪夫原型数字低通滤波器(Ⅰ型):
(1)求出其并联型结构;
(2)画出振幅、相位、振幅衰减、群延迟图。
四、实验说明
完成上述实验内容主要采用以下扩展函数:
1.巴特沃思型模拟低通滤波器设计函数:
function[b,a]=afd_butt(Wp,Ws,Rp,As)
其中:
b=Ha(s)分子的系数
a=Ha(s)分母的系数
Wp=以弧度/秒为单位的通带边缘频率;Wp>0
Ws=以弧度/秒为单位的阻带边缘频率;Ws>Wp>0
Rp=通带中的振幅波动的+dB数;(Rp>0)
As=阻带衰减的+dB数;(As>0)
2.未归一化巴特沃思型模拟低通滤波器原型设计函数:
function[b,a]=u_buttap(N,Omegac);
其中:
b=Ha(s)分子的系数
a=Ha(s)分母的系数
N=滤波器的阶数
Omegac=以弧度/秒为单位的截止频率
3.切贝雪夫Ⅰ型模拟低通滤波器设计函数:
function[b,a]=afd_chbl(Wp,Ws,Rp,As);
其中:
b=Ha(s)分子的系数
a=Ha(s)分母的系数
Wp=以弧度/秒为单位的通带边缘频率;Wp>0
Ws=以弧度/秒为单位的阻带边缘频率;Ws>Wp>0
Rp=通带中的振幅波动的+dB数;(Rp>0)
As=阻带衰减的+dB数;(As>0)
4.未归一化切贝雪夫Ⅰ型模拟低通滤波器原型设计函数:
function[b,a]=u_chblap(N,Rp,Omegac);
其中:
b=Ha(s)分子的系数
a=Ha(s)分母的系数
N=滤波器的阶数
Omegac=以弧度/秒为单位的截止频率
五、实验报告要求
阅读程序,并且加上详细的注释;
分别给出机算和笔算的结果;
分析所得结果并给出结论;
自行设计一符合某性能要求的滤波器,分析所得结果并给出结论。
六、附程
1.巴特沃思型模拟低通滤波器设计函数:
function[b,a]=afd_butt(Wp,Ws,Rp,As);
ifWp<=0
error('通带边缘必须大于0')
end
ifWs<=0
error('阻带边缘必须大于通带边缘')
end
if(Rp<=0)|(As<0)
error('通带波动或阻带衰减必须大于0')
end
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)));
fprintf('\n***ButterworthFilterOrder=%2.0f\n',N)
OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N)))
[b,a]=u_buttap(N,OmegaC);
2.未归一化巴特沃思型模拟低通滤波器原型设计函数:
function[b,a]=u_buttap(N,Omegac);
[z,p,k]=buttap(N);
p=p*Omegac;
k=k*Omegac^N;
B=real(poly(z));
bo=k;
b=k*B;
a=real(poly(p));
3.切贝雪夫Ⅰ型模拟低通滤波器设计函数:
function[b,a]=afd_chbl(Wp,Ws,Rp,As);
ifWp<=0
error('通带边缘必须大于0')
end
ifWs<=Wp
error('阻带边缘必须大于通带边缘')
end
if(Rp<=0)|(As<0)
error('通带波动或阻带衰减必须大于0')
end
ep=sqrt(10^(Rp/10)-1);
A=10^(As/20);
OmegaC=Wp;
OmegaR=Ws/Wp;
g=sqrt(A*A-1)/ep;
N=ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1)));
fprintf('\n***切比雪夫I型滤波器阶次=%2.0f\n',N)
[b,a]=u_chblap(N,Rp,OmegaC);
4.未归一化切贝雪夫Ⅰ型模拟低通滤波器原型设计函数:
function[b,a]=u_chblap(N,Rp,Omegac);
[z,p,k]=cheb1ap(N,Rp);
a=real(poly(p));
aNn=a(N+1);
p=p*Omegac;
a=real(poly(p));
aNu=a(N+1);
k=k*aNu/aNn;
b0=k;
B=real(poly(z));
b=k*B;
5.计算滤波器的幅值响应(绝对、相对)、相位响应、群延迟函数:
function[db,mag,pha,grd,w]=freqz_m(b,a);
[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);
grd=grpdelay(b,a,w);
6.将结果绘制成图形函数
functiony=ht(b,a,Wp,Ws,Attn,Ripple)
figure
(1);
[db,mag,pha,grd,w]=freqz_m(b,a);
subplot(2,2,1);plot(w/pi,mag);title('MagnitudeResponse')
xlabel('frequencyinpiunits');ylabel('magnitude');axis([0,1,0,1.1])
set(gca,'XTickMode','manual','XTick',[0,Wp/pi,Ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,Ripple,1]);grid
subplot(2,2,3);plot(w/pi,db);title('MagnitudeindB')
xlabel('frequencyinpiunits');ylabel('decibels');axis([0,1,-40,5])
set(gca,'XTickMode','manual','XTick',[0,Wp/pi,Ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[-40,-15,-5,0]);grid
set(gca,'YTickLabelMode','manual','YTickLabel',['40';'15';'5';'0'])
subplot(2,2,2);plot(w/pi,pha/pi);title('PhaseResponse')
xlabel('frequen
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 实验