数字信号处理仿真实验.docx
《数字信号处理仿真实验.docx》由会员分享,可在线阅读,更多相关《数字信号处理仿真实验.docx(28页珍藏版)》请在冰豆网上搜索。
数字信号处理仿真实验
数字信号处理上机实验仿真报告
学院电子工程学院
班级
学号
姓名
实验一:
信号、系统及系统响应
一、实验目的
(1)熟悉连续信号经理想采样前后的频谱变化关系,加深对时域采样定理的理解。
(2)熟悉时域离散系统的时域特性。
(3)利用卷积方法观察分析系统的时域特性。
(4)掌握序列傅里叶变换的计算机实现方法,利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
二、实验原理
采样是连续信号数字处理的第一个关键环节。
对一个连续信号
进行理想采样的过程可用(1.1)式表示。
(1.1)
其中
为
的理想采样,
为周期冲激脉冲,即
(1.2)
的傅里叶变换
为
(1.3)
将(1.2)式代入(1.1)式并进行傅里叶变换,
(1.4)
式中的
就是采样后得到的序列
,即
的傅里叶变换为
(1.5)
比较(1.5)和(1.4)可知
(1.6)
为了在数字计算机上观察分析各种序列的频域特性,通常对
在
上进行M点采样来观察分析。
对长度为N的有限长序列
,有
(1.7)
其中
一个时域离散线性时不变系统的输入/输出关系为
(1.8)
上述卷积运算也可以转到频域实现
(1.9)
三、实验内容及步骤
(1)认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容,阅读本实验原理与方法。
(2)编制实验用主程序及相应子程序。
①信号产生子程序,用于产生实验中要用到的下列信号序列:
xa(t)=Ae-atsin(Ω0t)u(t)
进行采样,可得到采样序列
xa(n)=xa(nT)=Ae-anTsin(Ω0nT)u(n),0≤n<50
其中A为幅度因子,a为衰减因子,Ω0是模拟角频率,T为采样间隔。
这些参数都要在实验过程中由键盘输入,产生不同的xa(t)和xa(n)。
b.单位脉冲序列:
xb(n)=δ(n)
c.矩形序列:
xc(n)=RN(n),N=10
②系统单位脉冲响应序列产生子程序。
本实验要用到两种FIR系统。
a.ha(n)=R10(n);b.hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)
③有限长序列线性卷积子程序,用于完成两个给定长度的序列的卷积。
可以直接调用MATLAB语言中的卷积函数conv。
conv用于两个有限长度序列的卷积,它假定两个序列都从n=0开始。
调用格式如下:
y=conv(x,h)
(3)调通并运行实验程序,完成下述实验内容:
①分析采样序列的特性。
a.取采样频率fs=1kHz,即T=1ms。
b.改变采样频率,fs=300Hz,观察|X(ejω)|的变化,并做记录(打印曲线);进一步降低采样频率,fs=200Hz,观察频谱混叠是否明显存在,说明原因,并记录(打印)这时的|X(ejω)|曲线。
②时域离散信号、系统和系统响应分析。
a.观察信号xb(n)和系统hb(n)的时域和频域特性;利用线性卷积求信号xb(n)通过系统hb(n)的响应y(n),比较所求响应y(n)和hb(n)的时域及频域特性,注意它们之间有无差别,绘图说明,并用所学理论解释所得结果。
b.观察系统ha(n)对信号xc(n)的响应特性。
③卷积定理的验证。
(4)主程序框图
①分析采样序列的特性
②时域离散信号、系统和系统响应分析
③卷积定理的验证
四.实验程序及对应波形
1.子程序
function[XN,n,k]=DFT(xn,N)
n=0:
N-1;
k=-200:
200;
XN=xn*exp(-j*2*pi/N).^(n'*k);%计算DFT[x(n)]
%产生矩形序列
functionx=juxing(n2);
x=[1,ones(1,n2)];
function[x,n]=maichong(n0,n1,n2)
n=(n1:
n2);
x=(n==n0);
%产生信号Xa(n)
functionx=xn(A,a,w,fs)
n=0:
50-1;
x=A*exp((-a)*n/fs).*sin(w*n/fs).*juxing(49);
functionx=u(t);
x=(t>=0);
%产生脉冲信号
function[x,n]=maichong(n0,n1,n2)
n=(n1:
n2);
x=(n==n0);
2.主程序
①分析采样序列的特性
程序代码如下(见后附程序Ex1_1):
A=100;a=200;w0=200;k=-200:
200;
T=0.001;t=0:
T:
0.06;N=50;k1=0:
1:
N;W1max=2*pi*500;W1=W1max*k1/N;w1=W1/pi;
xat=A*exp(-a*t).*sin(w0*t).*u(t);
Xa=xat*exp(-j*t'*W1);
figure
(1);
subplot(2,1,1);
plot(t,xat);
xlabel('t');
ylabel('xa(t)');
title('采样信号a(t)');
axis([0,0.06,-5,35]);
subplot(2,1,2);
plot(w1,abs(Xa));
xlabel('w');
ylabel('X(jw)');title('xa(t)的频谱');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fs=1000Hz
A=100;a=200;w0=200;k=-200:
200;
fs=1000;
w=k/50;
xan=xn(A,a,w0,fs);%产生信号xa(n)
X=DFT(xan,50);
figure
(2)
subplot(2,1,1)
n=0:
49;
stem(n,xan,'.');axis([0,50,-20,50]);
xlabel('n');
ylabel('xa(n)');
title('采样信号fs=1000Hz');
subplot(2,1,2);
plot(w,abs(X));
xlabel('w/pi');
ylabel('X(e^jw)');title('xa(n)的频谱(fs=1000Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fs=300Hz
fs=300;
xan=xn(A,a,w0,fs);%产生信号xa(n)
X=DFT(xan,50);
figure(3)
subplot(2,1,1)
n=0:
49;
stem(n,xan,'.');axis([0,50,-20,50]);
xlabel('n');
ylabel('xa(n)');
title('采样信号fs=300Hz');
subplot(2,1,2);
plot(w,abs(X));
xlabel('w/pi');
ylabel('X(e^jw)');title('xa(n)的频谱(fs=300Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fs=200Hz
fs=200;
xan=xn(A,a,w0,fs);%产生信号xa(n)
X=DFT(xan,50);
figure(4)
subplot(2,1,1)
n=0:
49;
stem(n,xan,'.');axis([0,50,-20,50]);
xlabel('n');
ylabel('xa(n)');
title('采样信号fs=200Hz');
subplot(2,1,2);
plot(w,abs(X));
xlabel('w/pi');
ylabel('X(e^jw)');title('xa(n)的频谱(fs=200Hz)');
仿真结果如下:
由图可见,在折叠频率w=π,即f=fs/2=500Hz处混叠很小。
当fs=300Hz时,存在较明显的混叠失真;当fs=200时,发生严重的混叠失真。
②时域离散信号、系统和系统响应分析
a.观察信号xb(n)和系统hb(n)的时域和频域特性;
实验代码如下:
(见后附程序Ex1_2)
%信号xb
N=64
xb=zeros(1,N)
xb
(1)=1;xn1=[0:
N-1]
figure
(1);subplot(2,1,1)
stem(xn1,xb,'.')
title('信号xb时域特性')
w=linspace(-5*pi,5*pi,1000)
Xb=xb*exp(-j*xn1'*w)
figure
(1)
subplot(2,1,2)
plot(w/pi,abs(Xb),'.')
title('信号xb频域特性')
%ϵͳ信号hb
hb=[1,2.5,2.5,1];nhb=[0:
3]
figure
(2);subplot(2,1,1)
stem(nhb,hb,'.')
title('系统函数ýhb时域特性')
w=linspace(-2.8*pi,2.8*pi,1000)
Hb=hb*exp(-j*nhb'*w)
figure
(2)
subplot(2,1,2)
plot(w/pi,abs(Hb),'.')
title('系统函数ýhb频域特性')
%响应
yb=conv(xb,hb);
nyb=[0:
length(yb)-1]
figure(3);subplot(2,1,1)
stem(yb,'.')
title('响应yb时域特性')
w=linspace(-2.8*pi,2.8*pi,1000)
Yb=yb*exp(-j*nyb'*w)
figure(3)
subplot(2,1,2)
plot(w/pi,abs(Yb),'.')
title('响应Yb频域特性')
仿真结果如下:
b.观察系统ha(n)对信号xc(n)的响应特性。
实验代码如下(见后附程序Ex1_3):
%信号xc
N=10
xc=ones(1,N)
xn1=[0:
N-1]
figure
(1);subplot(2,1,1)
stem(xn1,xc,'.')
title('信号xc时域响应')
w=linspace(-5*pi,5*pi,1000)
Xc=xc*exp(-j*xn1'*w)
figure
(1)
subplot(2,1,2)
plot(w/pi,abs(Xc))
%ϵͳ信号ha
N=10
ha=ones(1,N)
hn1=[0:
N-1]
figure
(2);subplot(2,1,1)
stem(hn1,ha,'.')
title('系统函数ha时域响应')
w=linspace(-2.8*pi,2.8*pi,1000)
Ha=ha*exp(-j*hn1'*w)
figure
(2)
subplot(2,1,2)
plot(w/pi,abs(Ha),'.')
title('系统函数ha频域响应')
%响应
y2=conv(xc,ha);
ny2=[0:
length(y2)-1]
figure(3);subplot(2,1,1)
stem(y2,'.')
title('响应y2时域响应')
w=linspace(-2.8*pi,2.8*pi,1000)
Y2=y2*exp(-j*ny2'*w)
figure(3)
subplot(2,1,2)
plot(w/pi,abs(Y2),'.')
title('响应Yb频域响应')
仿真结果如下:
③卷积定理的验证。
实验代码如下:
(见后附程序Ex1_4):
%时域卷积定理的验证
%信号xc
N=10
xc=ones(1,N)
xn1=[0:
N-1]
w=linspace(-5*pi,5*pi,1000)
Xc=xc*exp(-j*xn1'*w)
%系统信号ha
N=10
ha=ones(1,N)
hn1=[0:
N-1]
w=linspace(-5*pi,5*pi,1000)
Ha=ha*exp(-j*hn1'*w)
%响应
y2=conv(xc,ha);
ny2=[0:
length(y2)-1]
w=linspace(-5*pi,5*pi,1000)
Y2=y2*exp(-j*ny2'*w)
figure
(1)
subplot(2,1,1)
plot(w/pi,abs(Y2))
title('响应Yb=DTFT(conv(xc,ha))频域特性')
Y=Xc.*Ha
subplot(2,1,2)
plot(w/pi,abs(Y))
title('响应Y=Xc.*Ha频域特性')
仿真结果如下:
两种方法得到的谱是一样的,即验证了卷积定理。
五、思考题
(1)在分析理想采样序列特性的实验中,采样频率不同时,相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同?
它们所对应的模拟频率是否相同?
为什么?
答:
由
可知,若采样频率不同,则其周期T不同,相应的数字频率
也不相同;而因为是同一信号,故其模拟频率
保持不变。
(2)在卷积定理验证的实验中,如果选用不同的频域采样点数M值,例如,选M=10和M=20,分别做序列的傅里叶变换,求得
所得结果之间有无差异?
为什么?
答:
有差异。
因为所得
图形由其采样点数唯一确定,由频域采样定理可知,若M小于采样序列的长度N,则恢复原序列时会发生时域混叠现象。
实验二:
用FFT作谱分析
一、实验目的:
(1)进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法,所以FFT的运算结果必然满足DFT的基本性质)。
(2)熟悉FFT算法原理和FFT子程序的应用。
(3)学习用FFT对连续信号和时域离散信号进行谱分析的方法,了解可能出现的分析误差及其原因,以便在实际中正确应用FFT。
二、实验步骤:
(1)复习DFT的定义、性质和用DFT作谱分析的有关内容。
(2)复习FFT算法原理与编程思想,并对照DIT—FFT运算流图和程序框图,读懂本实验提供的FFT子程序。
(3)编制信号产生子程序,产生以下典型信号供谱分析用:
=R4(n)
=cos(pi/4*n)
=sin(pi/8*n)
=cos(pi*8*t)+sin(pi*16*t)+cos(20*pi*t)
应当注意,如果给出的是连续信号xa(t),则首先要根据其最高频率确定采样速率fs以及由频率选择采样点数N,然后对其进行软件采样(即计算x(n)=xa(nT),0<=n<=N-1),产生对应序列x(n)。
对信号x6(t),频率分辨率的选择要以能分辨开其中的三个频率对应的谱线为准则。
对周期序列,最好截取周期的整数倍进行分析,否则有可能产生较大的分析误差。
请实验者根据DFT的隐含周期性思考这个问题。
(4)编写主程序。
下图给出了主程序框图,供参考。
本实验提供FFT子程序和通用绘图子程序。
三.实验程序及对应波形
(1)对2中所给出的信号逐个进行谱分析。
实验代码如下(见后附程序Ex2_1)
%信号x1
x1=[1,1,1,1]
n1=[0:
3]
figure
(1);subplot(3,1,1)
stem(n1,x1)
title('信号x1时域特性')
w=linspace(-5*pi,5*pi,1000)
X1=x1*exp(-j*n1'*w)
figure
(1)
subplot(3,1,2)
plot(w/pi,abs(X1))
title('DTFT分析信号x1频域特性')
X11=fftshift(fft(x1,128))
figure
(1)
subplot(3,1,3)
stem(abs(X11))
title('fft分析信号x1频域特性')
%信号x2
x2=[1,2,3,4,4,3,2,1]
n2=[0:
7]
figure
(2);subplot(3,1,1)
stem(n2,x2)
title('信号x2时域特性')
w=linspace(-5*pi,5*pi,1000)
X2=x2*exp(-j*n2'*w)
figure
(2)
subplot(3,1,2)
plot(w/pi,abs(X2))
title('DTFT分析信号x2频域特性')
X22=fftshift(fft(x2,128))
figure
(2)
subplot(3,1,3)
stem(abs(X22))
title('fft分析信号x2频域特性')
%信号x3
x3=[43211234]
n3=[0:
7]
figure(3);subplot(3,1,1)
stem(n3,x3)
title('信号x3时域特性')
w=linspace(-5*pi,5*pi,1000)
X3=x3*exp(-j*n3'*w)
figure(3)
subplot(3,1,2)
plot(w/pi,abs(X3))
title('DTFT信号x3频域特性')
X33=fftshift(fft(x3,128))
figure(3)
subplot(3,1,3)
stem(abs(X3))
title('fft分析信号x3频域特性')
%信号x4=cos(pi/4*n)
n4=[0:
15]
x4=cos(pi/4*n4)
figure(4);subplot(3,1,1)
stem(n4,x4)
title('信号x4时域特性')
w=linspace(-5*pi,5*pi,1000)
X4=x4*exp(-j*n4'*w)
figure(4)
subplot(3,1,2)
plot(w/pi,abs(X4))
title('DTFT信号x4频域特性')
X44=fftshift(fft(x4,128))
figure(4)
subplot(3,1,3)
stem(abs(X44))
title('fft分析信号x4频域特性')
%信号x5=sin(pi/8*n)
n5=[0:
15]
x5=sin(pi/8*n5)
figure(5);subplot(3,1,1)
stem(n5,x5)
title('信号x5时域特性')
w=linspace(-5*pi,5*pi,1000)
X5=x5*exp(-j*n5'*w)
figure(5)
subplot(3,1,2)
plot(w/pi,abs(X5))
X55=fftshift(fft(x5,128))
figure(5)
subplot(3,1,3)
stem(abs(X55))
title('fft分析信号x4频域特性')
%信号x6=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t)
t=0:
0.0001:
2.5
xa=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t)%原始的连续时间信号xa
figure(6);subplot(3,1,1);
plot(t,xa);title('连续时间信号曲线xa');;axis([02.5-55])%画连续时间信号曲线
T=0.01;n6=0:
2.5/T;%建立离散自变量向量
x6=cos(8*pi*n6*T)+cos(16*pi*n6*T)+cos(20*pi*n6*T)%采样周期为T的离散时间信号x
figure(6);subplot(3,1,2);
stem(n6*T,x6);title('取样信号');axis([02.5-55])%画离散时间信号曲线
w=linspace(-5*pi,5*pi,2500);%均匀取点
X6=x6*exp(-j*n6'*w)%计算DTFT
figure(6);subplot(3,1,3);
plot(w/pi,abs(X6));title('DTFT分析取样信号频谱');
仿真结果如下:
(2)令x(n)=x4(n)+x5(n),用FFT计算8点和16点离散傅里叶变换,X(k)=DFT[x(n)]
(3)令x(n)=x4(n)+jx5(n),重复
(2)。
代码如下:
(见后附程序Ex2_2)
%信号x4=cos(pi/4*n)
n4=[0:
15]
x4=cos(pi/4*n4)
figure
(1);subplot(2,1,1)
stem(n4,x4)
title('信号x4时域特性')
w=linspace(-5*pi,5*pi,1000)
X4=x4*exp(-j*n4'*w)
figure
(1)
subplot(2,1,2)
plot(w/pi,abs(X4))
title('DTFT分析信号x4频域特性')
%信号x5=sin(pi/8*n)
n5=[0:
15]
x5=sin(pi/8*n5)
figure
(2);subplot(2,1,1)
stem(n5,x5)
title('DTFT分析信号x5时域特性')
w=linspace(-5*pi,5*pi,1000)
X5=x5*exp(-j*n5'*w)
figure
(2)
subplot(2,1,2)
plot(w/pi,abs(X5))
%x(n)=x4(n)+x5(n)
xn1=x4+x5
nx1=0:
length(xn1)-1
figure(3);subplot(3,1,1)
stem(nx1,xn1)
title('xn1=x4+x5时域特性')
Xn1_8=fftshift(fft(xn1,8))
n1_8=0:
7
figure(3);subplot(3,1,2)
stem(n1_8,abs(Xn1_8))
title('xn1的8点fft变换')
Xn1_16=fftshift(fft(xn1,16))
n1_16=0:
15
figure(3);subplot(3,1,3)
stem(n1_16,abs(Xn1_16))
title('xn1的16点fft变换')
%x(n)=x4(n)+jx5(n)
xn2=x4+j*x5
nx2=0:
length(xn2)-1
figure(4);subplot(3,1,1)
stem(nx2,xn2)
title('xn2=x4+j*x5时域特性')
Xn2_8=fftshift(fft(xn2,8))
n2_8=0:
7
figure(4);subplot(3,1,2)
stem(n2_8,abs(Xn2_8))
title('xn2的8点fft变换')
Xn2_16=fftshift(fft(xn2,16))
n2_16=0:
15
figure(4);subplot(3,1,3)
stem(n2_16,abs(Xn2_16))
title('xn2的16点fft变换')
仿真结果如下:
四、思考题
(1)在N=8时,x2(n)和x3(n)的幅频特性会相同吗?
为什么?
N=16呢?
答:
在N=8时,x2(n)和x3(n)的幅频特性会相同在N=8时x2(n)和x