数字信号matlab试验程序以及实验报告.docx
- 文档编号:11898440
- 上传时间:2023-04-08
- 格式:DOCX
- 页数:30
- 大小:283.44KB
数字信号matlab试验程序以及实验报告.docx
《数字信号matlab试验程序以及实验报告.docx》由会员分享,可在线阅读,更多相关《数字信号matlab试验程序以及实验报告.docx(30页珍藏版)》请在冰豆网上搜索。
数字信号matlab试验程序以及实验报告
数字信号处理课程设计实验报告
基础实验部分
实验一离散时间系统及离散卷积
一、实验目的
(1)熟悉MATLAB软件的使用方法。
(2)熟悉系统函数的零极点分布、单位脉冲响应和系统频率响应等概念。
(3)利用MATLAB绘制系统函数的零极点分布图、系统频率响应和单位脉冲响应。
(4)熟悉离散卷积的概念,并利用MATLAB计算离散卷积。
二、实验原理
本实验利用MATLAB的库函数,实验冲击响应零极点、幅度相位谱的绘制。
熟悉MATLAB的使用。
并利用离散卷积的概念编写一个计算离散卷积的程序。
三、实验内容及要求
1.离散时间系统的单位冲击响应
(1)选择一个离散时间系统;
(2)用笔进行差分方程的递推计算;
(3)编制差分方程的递推计算程序;
(4)在计算机上实现递推运算;
(5)将程序计算结果与笔算的计算结果进行比较,验证程序运行的正确性;
2.离散系统的幅频、相频的分析方法
(1)离散系统的幅频、相频的分析方法
(2)用笔计算几个特殊的幅频、相频的值,画出示意曲线图;
(3)编制离散系统的幅频、相频的分析程序;
(4)在计算机上进行离散系统的幅频、相频特性计算,并画出曲线;
(5)通过比较,验证程序的正确性;
3.离散卷积的计算
(1)选择两个有限长序列,用笔计算其线性卷积;
(2)编制有限长序列线性卷积程序;
(3)利用计算程序对
(1)选择的有限长序列进行卷积运算;
(4)比较结果验证程序的正确性。
四、实验过程
(1)离散时间系统的单位脉冲响应:
a=[1,-1,0.91]; %系统的系数
b=[1,2,3];
n0=0;
n=-20:
120;
x=[(n-n0)==0]; % 利用filter函数进行滤波,从而求得系统的单位冲击响应
h=filter(b,a,x);
figure
(1)
stem(n,h); %利用stem绘制曲线图
title('冲击响应');
xlabel('n');
ylabel('h(n)');
h
figure
(2)
[z,p,g]=tf2zp(b,a); %求传递函数
zplane(z,p) %求零极点
(2)离散系统的幅频、相频的分析:
b=[0.6,1.2,0.4,0.9];%分子的系数
a=[1.0,-1.7,1.1,-0.5]; %分母的系数
m=0:
length(b)-1;
l=0:
length(a)-1;
K=500;
k=1:
K;
w=pi*k/K;
H=(b*exp(-j*m'*w))./(a*exp(-j*l'*w)); %求传递函数
magH=abs(H); %求传递函数的幅值
angH=angle(H); %求传递函数的相位
figure
(1)
subplot(2,1,1);
plot(w/pi,magH); %定义坐标轴的区间
grid;
axis([0,1,0,1]);
xlabel('w(pi)');
ylabel('|H|');
title('·幅度,相位响应');
subplot(2,1,2);
plot(w/pi,angH);
grid;
xlabel('w(pi)');
ylabel('angle(H)');
(3)离散卷积的计算
functionpr13()
n=-5:
50;
u1=stepseq(0,-5,50);
u2=stepseq(10,-5,50);
x=u1-u2;
h=((0.9).^n).*u1;
figure
(1)
subplot(3,1,1);
stem(n,x);
axis([-5,50,0,2]);
title('输入序列');
xlabel('n');
ylabel('x(n)');
subplot(3,1,2);
stem(n,h);
axis([-5,50,0,2]);
title('冲击响应序列');
xlabel('n');
ylabel('h(n)');
[y,ny]=conv_m(x,n,h,n);
subplot(3,1,3);
stem(ny,y);
axis([-5,50,0,8]);
title('输出响应');
xlabel('n');
ylabel('y(n)');
function[x,n]=stepseq(n0,n1,n2)
n=n1:
n2;
x=[(n-n0)>=0];
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、加深理解离散傅立叶变换及快速傅立叶变换概念;
2、学会应用FFT对典型信号进行频谱分析的方法;
3、研究如何利用FFT程序分析确定性时间连续信号;
4、熟悉应用FFT实现两个序列的线性卷积的方法。
二、实验原理
在各种信号序列中,有限长序列信号处理占有很重要地位,对有限长序列,我们可以使用离散Fouier变换(DFT)。
这一变换不但可以很好的反映序列的频谱特性,而且易于用快速算法在计算机上实现,当序列x(n)的长度为N时,它的DFT定义为
反变换为
有限长序列的DFT是其Z变换在单位圆上的等距采样,或者说是序列Fourier变换的等距采样,因此可以用于序列的谱分析。
FFT并不是与DFT不同的另一种变换,而是为了减少DFT运算次数的一种快速算法。
它是对变换式进行一次次分解,使其成为若干较短序列的组合,从而减少运算量。
常用的FFT是以2为基数的,其长度
。
它的效率高,程序简单,使用非常方便,当要变换的序列长度不等于2的整数次方时,为了使用以2为基数的FFT,可以用末位补零的方法,使其长度延长至2的整数次方。
(一)、在运用DFT进行频谱分析的过程中可能产生三种误差:
(1) 混叠
序列的频谱是被采样信号频谱的周期延拓,当采样速率不满足Nyquist定理时,就会发生频谱混叠,使得采样后的信号序列频谱不能真实的反映原信号的频谱。
避免混叠现象的唯一方法是保证采样速率足够高,使频谱混叠现象不致出现,即在确定采样频率之前,必须对频谱的性质有所了解。
在一般情况下,为了保证不出现频谱混叠,在采样前,先进行抗混叠滤波。
(2) 泄漏
实际中我们往往用截短的序列来近似很长的甚至是无限长的序列,这样可以使用较短的DFT来对信号进行频谱分析,这种截短等价于给原信号序列乘以一个矩形窗函数,也相当于在频域将信号的频谱和矩形窗函数的频谱卷积,所得的频谱是原序列频谱的扩展。
泄漏不能与混叠完全分开,因为泄漏导致频谱的扩展,从而造成混叠。
为了减少泄漏的影响,可以选择适当的窗函数使频谱的扩散减至最小。
(3) 栅栏效应
DFT是对单位圆上Z变换的均匀采样,所以它不可能将频谱视为一个连续函数,就一定意义上看,用DFT来观察频谱就好像通过一个栅栏来观看一个图景一样,只能在离散点上看到真实的频谱,这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所拦住,不能别我们观察到。
减小栅栏效应的一个方法就是借助于在原序列的末端填补一些零值,从而变动DFT的点数,这一方法实际上是人为地改变了对真实频谱采样的点数和位置,相当于搬动了每一根“尖桩栅栏”的位置,从而使得频谱的峰点或谷点暴露出来。
(二)、用FFT计算线性卷积
用FFT可以实现两个序列的圆周卷积。
在一定的条件下,可以使圆周卷积等于线性卷积。
一般情况,设两个序列的长度分别为N1和N2,要使圆周卷积等于线性卷积的充要条件是FFT的长度
N≥N1+N2
对于长度不足N的两个序列,分别将他们补零延长到N。
当两个序列中有一个序列比较长的时候,我们可以采用分段卷积的方法。
通常,我们选择重叠相加法。
重叠相加法:
将长序列分成与短序列相仿的片段,分别用FFT对它们作线性卷积,再将分段卷积各段重叠的部分相加构成总的卷积输出。
三、实验内容
1、用离散傅立叶变换程序处理时间抽样信号,并根据实序列离散傅立叶变换的对称性,初步判定程序的正确性;
2、观察三角波和反三角波序列的时域和幅频特性,用N=8点FFT分析信号序列
和
的幅频特性,观察两者的序列形状和频谱曲线有什么异同?
绘出两序列及其幅频特性曲线。
三角波序列
反三角波序列
3、已知余弦信号如下
当信号频率
,采样间隔
,采样长度
时,对该信号进行傅立叶变换。
用FFT程序分析正弦信号,分别在以下情况下进行,并且分析比较结果
(1)F=50,N=32,T=0.000625;
(2)F=50,N=32,T=0.005;
(3)F=50,N=32,T=0.0046875;
(4)F=50,N=32,T=0.004;
(5)F=50,N=64=0.000625
4、选定某一时间信号进行N=64点离散傅立叶变换,并且,对同一信号进行快速傅立叶变换,并比较它们的速度。
四、实验要求
1、调试实验程序,并且,给参考程序加注释;
2、完成实验内容2,并对结果进行分析。
实验中的信号序列
和
,在单位圆上的Z变换频谱
和
会相同吗?
如果不同,你能说出哪一个低频分量更多一些吗?
为什么?
3、完成实验内容3,并对结果进行分析;
4、利用编制的计算卷积的计算程序,分别给出一下三组函数的卷积结果。
(1)
(2)
(3)
、
五、实验过程
functionshiyan2_1()
%此程序用于观察三角波和放三角波的时域和幅频特性;
%xc(n)为三角波,xd(n)为反三角波;
clearall;
clc;
N=64;
n=1:
N;
xc=zeros(1,N);
xd=zeros(1,N);
%生成三角波和反三角波
fori=1:
4
xc(i+28)=i;
end
forj=5:
8
xc(j+28)=9-j;
end
fori=1:
4
xd(i+28)=5-i;
end
forj=5:
8
xd(j+28)=j-4;
end
figure
(1);
subplot(2,1,1),stem(n,xc);
title('三角波的时域波形');
subplot(2,1,2),stem(n,xd);
title('反三角波的时域波形');
%FFT变换
XC=fft(xc,N);
XD=fft(xd,N);
AXC=abs(XC);
AXD=abs(XD);
figure
(2);
subplot(2,1,1),plot(AXC);
title('三角波的幅度谱波形');
subplot(2,1,2),plot(AXD);
title('反三角波的幅度谱波形');
functionshiyan2_2()
%此程序用来对函数x(t)=cos(2*pi*F*t),进行不同频率F、不同采样间隔T、不同采样点N的FFT分析;
%T<=2*(1/F) 0.04时满足抽样定理;
clearall;
clc;
%F=50,N=64,T=0.000625
F=50;
N=64;
T=0.000625;
n=1:
N;
x=cos(2*pi*F*n*T);
figure
(1);
subplot(2,2,1);
plot(n,x);
title('x(n)F=50,N=64,T=0.000625');
xlabel('n');
X=fft(x);
magX=abs(X);
subplot(2,2,3);
plot(n,X);
title('FFT|X|');
xlabel('f(pi)');
%F=50,N=32,T=0.000625
F=50;
N=32;
T=0.000625;
n=1:
N;
x=cos(2*pi*F*n*T);
subplot(2,2,2);
plot(n,x);
title('x(n)F=50,N=32,T=0.000625');
xlabel('n');
X=fft(x);
magX=abs(X);
subplot(2,2,4);
plot(n,X);
title('FFT|X|');
xlabel('f(pi)');
%F=50,N=32,T=0.005
F=50;
N=32;
T=0.005;
n=1:
N;
x=cos(2*pi*F*n*T);
figure
(2);
subplot(2,2,1);
plot(n,x);
title('x(n)F=50,N=32,T=0.005');
xlabel('n');
X=fft(x);
magX=abs(X);
subplot(2,2,3);
plot(n,X);
title('FFT|X|');
xlabel('f(pi)');
%F=50,N=32,T=0.0046875
F=50;
N=32;
T=0.0046875;
n=1:
N;
x=cos(2*pi*F*n*T);
subplot(2,2,2);
plot(n,x);
title('x(n)F=50,N=32,T=0.0046875');
xlabel('n');
X=fft(x);
magX=abs(X);
subplot(2,2,4);
plot(n,X);
title('FFT|X|');
xlabel('f(pi)');
%F=50,N=32,T=0.004
F=50;
N=32;
T=0.004;
n=1:
N;
x=cos(2*pi*F*n*T);
figure(3);
subplot(2,1,1);
plot(n,x);
title('x(n)F=50,N=32,T=0.004');
xlabel('n');
X=fft(x);
magX=abs(X);
subplot(2,1,2);
plot(n,X);
title('FFT|X|');
xlabel('f(pi)');
functionshiyan2_3(N,yun)
%此程序用于比较离散傅立叶变换与同一信号进行快速傅立叶变换处理时间;
%N为测试的数据点数,yun为0/1关闭/开启绘图
clc;
%离散傅立叶变换处理时间
F=50;
%N=64;
T=0.000625;
n=1:
N;
x=cos(2*pi*F*n*T);
ifyun==1
subplot(3,1,1);
stem(n,x);
title('x(n)');
xlabel('n');
end
t1=cputime;
X=dft(x,N);
tt1=cputime-t1;
fprintf('测试的离散数据点为:
\n');
disp(N);
disp('离散傅立叶变换耗时为');
disp(tt1);
magX=abs(X);
ifyun==1
subplot(3,1,2);
stem(n,X);
title('DFT|X|离散傅立叶变换');
xlabel('f(pi)');
end
%快速傅立叶变换处理时间
t2=cputime;
FFTX=fft(x,N);
tt2=cputime-t2;
disp('快速傅立叶变换耗时为');
disp(tt2);
magfftx=abs(FFTX);
ifyun==1
subplot(3,1,3);
stem(n,FFTX);
title('DFT|X|快速傅立叶变换');
xlabel('f(pi)');
end
%dft的子程序
function[Xk]=dft(xn,N)
n=0:
N-1;
k=0:
N-1;
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
functionshiyan2_4()
%此程序用于计算两组卷积;
%用FFT计算卷积,首先将时域的x(n)与h(n)变换到频域;
%根据卷积定理时域卷积对于频域成绩,然后将结果反变换回去即可;
clearall;
clc;
%第一组;
N1=32;
n1=1:
14;
x1n=ones(1,14);
h1n=(4/5).^n1;
figure
(1);
X1=fft(x1n,N1);
H1=fft(h1n,N1);
Y1=X1.*H1;
y1=ifft(Y1);
stem(abs(y1));
title('第一组y(n)=x(n)*h(n))N=32点的卷积');
%第二组
N2=32;
x2n([1:
10])=1;
n2=1:
20;
h2n=0.5*sin(0.5*n2);
figure
(2);
X2=fft(x2n,N2);
H2=fft(h2n,N2);
Y2=X2.*H2;
y2=ifft(Y2);
stem(abs(y2));
title('第二组y(n)=x(n)*h(n))N=32点的卷积');
实验三 IIR数字滤波器设计
一、实验原理
(1)脉冲响应不变法
用数字滤波器的单位脉冲响应序列
模仿模拟滤波器的冲激响应
让
正好等于
的采样值,即
其中T为采样间隔,如果以
及
分别表示
的拉氏变换及
的Z变换,则
(2)双线性变换法
s平面与z平面之间满足以下映射关系:
s平面的虚轴单值地映射于z平面的单位圆上,s平面的左半平面完全映射到z平面的单位圆内。
双线性变换不存在混叠问题。
双线性变换时一种非线性变换
,这种非线性引起的幅频特性畸变可通过预畸而得到校正。
以低通数字滤波器为例,将设计步骤归纳如下:
1. 确定数字滤波器的性能指标:
通带临界频率fp、阻带临界频率fs;通带内的最大衰减Rp;阻带内的最小衰减As;采样周期T;
2. 确定相应的数字角频率,ωp=2πfpT;ωs=2πfrT;
3. 计算经过预畸的相应模拟低通原型的频率,
根据Ωp和Ωs计算模拟低通原型滤波器的阶数N,并求得低通原型的传递函数Ha(s);
4. 用上面的双线性变换公式代入Ha(s),求出所设计的传递函数H(z);
5. 分析滤波器特性,检查其指标是否满足要求。
二、实验目的
1、学习模拟-数字变换滤波器的设计方法;
2、掌握双线性变换数字滤波器设计方法;
3、掌握实现数字滤波器的具体方法。
三、实验步骤
1、设计一个巴特沃思数字低通滤波器,设计指标如下:
通带内
幅度衰减不大于1dB;阻带
幅度衰减不小于15dB;
2、编制计算设计的数字滤波器幅度特性和相位特性的程序,并进行实验验证。
3、编制实现该数字滤波器程序并且实现数字滤波
(1)分别让满足所设计的滤波器的通带、过渡带、阻带频率特性的正弦波通过滤波器,验证滤波器性能;
(2)改变正弦抽样时间,验证数字低通滤波器的模拟截止频率实抽样时间的函数。
四、实验过程
functionpr7()
%31IIRfilter
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
T=1;
Fs=1/T;
OmegaP=(2/T)*tan(wp/2);
OmegaS=(2/T)*tan(ws/2);
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep.^2));
Attn=1/10^(As/20);
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS)));
OmegaC=OmegaP/((10.^(Rp/10)-1).^(1/(2*N)));
[cs,ds]=u_buttap(N,OmegaC);
[b,a]=bilinear(cs,ds,Fs);
[mag,db,pha,w]=freqz_m(b,a);
subplot(3,1,1);
plot(w/pi,mag);
title('幅度响应');
xlabel('w(pi)');
ylabel('H');
axis([0,1,0,1.1]);
set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);
set(gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1]);
grid;
subplot(3,1,2);
plot(w/pi,db);
title('幅度响应(dB)');
xlabel('w(pi)');
ylabel('H');
axis([0,1,-40,5]);
set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);
set(gca,'YTickmode','manual','YTick',[-50,-15,-1,0]);
grid;
subplot(3,1,3);
plot(w/pi,pha);
title('相位响应');
xlabel('w(pi)');
ylabel('piunit');
%axis([0,1,0,1.1]);
set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);
set(gca,'YTickmode','manual','YTick',[-1,0,1]);
grid;
function[b,a]=u_buttap(N,OmegaC)
[z,p,k]=buttap(N);
p=p*OmegaC;
k=k*OmegaC.^N;
B=real(poly(z));
b0=k;
b=k*B;
a=real(poly(p));
function[mag,db,pha,w]=freqz_m(b,a)
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:
501))';
w=(w(1:
501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
实验四 FIR数字滤波器设计
一、实验原理
窗函数法设计线性相位FIR滤波器步骤
●确定数字滤波器的性能要求:
临界频率{ωk},滤波器单位脉冲响应长度N;
●根据性能要求,合理选择单位脉冲响应
的奇偶对称性,从而确定理想频率响应
的幅频特性和相频特性;
●求理想单位脉冲响应
,在实际计算中,可对
按M(M远大于N)点等距离采样,并对其求IDFT得
,用
代替
;
●选择适当的窗函数
,根据
求所需设计的FIR滤波器单位脉冲响应;
●求
,分析其幅频特性,若不满足要求,可适当改变窗函数形式或长度N,重复上述设计过程,以得到满意的结果。
窗函数的傅式变换
的主瓣决定了
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 matlab 试验 程序 以及 实验 报告