DSP实验报告word版
实验一 离散系统的时域分析
一、实验目的
1、掌握离散时间信号的MATLAB 表示；
2、信号运算；
3、差分方程的求解；
4、离散时间信号的卷积运算。
二、实验原理
1、离散时间信号
离散时间信号只在某些离散的瞬时给出函数值，而在其他时刻无定义。它是时间上不连续按一定先后次序排列的一组数的集合，称为时间序列，用x(n)表示，n 取整数代表时间的离散时刻。
在matlab 中用向量来表示一个有限长度的序列。
2、序列的类型
为了分析的方便，在数字信号处理中规定了一些基本的序列。 a) 单位采样序列 function [x,n]=impseq(n1,n2,n0) n=[n1:n2];
x=[(n-n0)==0]; 调用该函数 [x,n]=impseq(-2,8,2);
stem(n,x)
00
10
()00
1()0n n n n n
n n n n δδ =?=
? ≠?
=?-?
≠?
单位采样序列的另一种生成方法
n0=-2; n=[-10:10]; nc=length(n); x=zeros(1,nc); for i=1:nc if n(i)==n0 x(i)=1 end end stem(n,x)
b) 单位阶跃序列
function [x,n]=stepseq(n1,n2,n0) n=[n1:n2];
x=[(n-n0)>=0];
调用该函数 [x,n]=stepseq(-2,8,2); stem(n,x)
000
10()001()
0n n n n n n n n n
εε >=?=? =?-?
c) 实数指数序列
x(n)=an (运算符“.^”)
n=[0:10]; x=0.9.^n; stem(n,x)
d) 复数指数序列 n=[-10:10]; alpha=-0.1+0.3*j; x=exp(alpha*n);
real_x=real(x); image_x=imag(x); mag_x=abs(x); phase_x=angle(x); subplot(2,2,1); stem(n,real_x) subplot(2,2,2); stem(n,image_x) subplot(2,2,3); stem(n,mag_x) subplot(2,2,4); stem(n,phase_x)
()()j n
x n e αω+=(0.1j0.3)n x(n)e (10n 10)
-+= -<<
e) 正弦和余弦序列 n=[0:10];
x=3*cos(0.1*pi*n+pi/3); stem(n,x)
f)
随机序列
rand(1,N)产生其元素在[0，1]之间均匀分布长度为N 的随机序列。 randn(1,N)产生均值为0，方差为1，长度为N 的高斯随机序列。 3、信号运算
a) 信号时移 b) 信号倒置 c) 信号加 d) 信号乘 e) 信号微分 f)
信号积分等
0()sin()
x n n ωθ=+0()cos()
x n n ωθ=+
4、序列运算
在matlab 中进行序列运算要求参与运算的序列的长度要一样，如果出现长度不一样或者长度相同但采样位置不同时，不能直接利用加减运算符，利用matlab 中的下标运算使其具有相同的长度。
function [y,n]=sigadd(x1,n1,x2,n2)
n=min(min(n1),min(n2)) : max(max(n1),max(n2)); y1=zeros(1,length(n)); y2=y1;
y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; y=y1+y2;
移位
function [y,n]=sigshift(x,m,n0) n=m+n0; y=x;
例： 求出下列波形
x1(n)=2x(n-5)-3x(n+4)
n=[-2:10]; x=[1:7,6:-1:1];
[x11,n11]=sigshift(x,n,5); [x12,n12]=sigshift(x,n,-4);
[x1,n1]=sigadd(2*x11,n11,-3*x12,n12); stem(n1,x1)
5、差分方程
在《信号与系统》和《数字信号处理》课程中，我们知道描述线性移不变离散时间系统的数学模型是常系数差分方程，它与系统的结构流图之间可以互相推导。迭代解法（也称递推解法）是求解差分方程的最简单也最适用的方法，也是实现数字滤波器的一种基本方法。
(){1,2,3,4,5,6,7,6,5,4,3,2,1}x n
差分方程通式为：
x(n)与y(n)分别为系统的激励和响应。 6、差分方程MATLAB 实现
MATLAB 以函数filter(num , den , x)，来计算在给定输入和差分方程系数时求差分方程的数值解。
num,den 分别为系统方程的系数向量。 x 是输入序列。 已知某一系统方程为：
y[n]-y[n-1]+0.9y[n-2]=x[n]
计算并画出脉冲响应h(n),n=(-20,100)
n=[-20:100];
num=[1]; den=[1 -1 0.9]; x=impseq(-20,100,0); h=filter(num,den,x); stem(n,h)
xlabel('时间序号N'); ylabel('脉冲响应h'); title('脉冲响应');
离散卷积的计算公式如下： 序列x(n)、h(n) 可以是有限长或无限长，但为了在计算机上绘图观察方便，我们主要讨论有限长序列。如果x(n) 和h(n)长度分别为M 和N ，则响应序列y(n) 也为有限长序列，其长度为 L=M+N-1 。于是，上式可以“形象”地描述为两个有限长序列的反褶、移位、相乘、累加过程，这使计算机编程十分方便。 卷积函数conv(a,b)
实现两个序列a,b 的卷积。 例：假定两个序列。
x=[3,11,7,0,-1,4,2]; y=[2,3,0,-5,2,1]; c=conv(x,y); stem(c)
[]()[]()
N N
k r a k y n k b r x n r ==-=-∑∑
()()*()()()
m y n x n h n x m h n m +∞
=-∞
==-∑
将函数conv 稍加扩展为函数conv-m ，它可以对任意的序列求卷积。格式如下：
function [y,ny]= conv_m(x,nx,h,nh,p) %信号处理的改进卷积程序 nyb=nx(1)+nh(1);
nyc=nx(length(x))+nh(length(h)); ny=[nyb:p:nyc]; y=conv(x , h);
? 已知
试求卷积C(t)=f1(t)*f2(t),并绘制出f1、f2及卷积以后的波形。
p=0.1;
t1= [0:p:1]; f1=t1.*(t1>0); t2= [-1:p:2];
f2=t2.*exp(-t2).*(t2>=0)+exp(t2).*(t2<0); [y,ny]=conv_m(f1,t1,f2,t2,p); Subplot(3,1,1); stem(t1,f1) Subplot(3,1,2); stem(t2,f2) Subplot(3,1,3); stem(ny,y)
12()()01
,0()12
,0t
t f t t t t te t f t t e t ε-= ≤≤?>=?= -≤≤?
四、实验报告要求
1、简述实验目的和实验原理。
2、用笔算求出你选定的序列x(n)、h(n)的卷积结果并与计算机计算结果相比较。
实验二 离散傅立叶变换
一、实验目的
1. 掌握离散傅里叶变换的有关性质。
2. 利用matlab 验证有关性质。
3. 利用傅立叶变换进行相关运算。
二、实验原理及方法
在工程技术的许多分支中，要掌握的基本内容之一就是正确理解时域和频域的关系。对于数字系统来说，就是要精通离散傅立叶变换，因此离散傅立叶变换在数字信号处理中占有十分重要的地位。在实际应用中，有限长序列有相当重要的地位，由于计算机容量的限制，只能对过程进行逐段分析。由于有限长序列，引入DFT(离散付里叶变换)。 傅里叶变换
建立以时间t 为自变量的“信号”与以频率f 为自变量的“频率函数”(频谱)之间的某种变换关系。所以“时间”或“频率”取连续还是离散值, 就形成各种不同形式的傅里叶变换对。
四种不同傅里叶变换对
1. 傅里叶级数(FS)：连续时间, 离散频率的傅里叶变换。周期连续时间信号傅里叶级
数(FS)得到非周期离散频谱密度函数。
2. 傅里叶变换(FT)：连续时间, 连续频率的傅里叶变换。非周期连续时间信号通过连
续付里叶变换(FT)得到非周期连续频谱密度函数。
3. 序列的傅里叶变换(DTFT):离散时间,连续频率的傅里叶变换。非周期离散的时间信
号(单位园上的Z 变换(DTFT))得到周期性连续的频率函数。 4. 离散傅里叶变换(DFT):离散时间, 离散频率的傅里叶变换。
上面讨论的前三种傅里叶变换对,都不适用在计算机上运算, 因为至少在一个域( 时域或频域)中, 函数是连续的。因为从数字计算角度我们感兴趣的是时域及频域都是离散的情况, 这就是第四种离散傅里叶变换。 离散傅里叶级数(DFS)
设 为周期为N 的周期序列, 则其离散傅里叶级数(DFS) 变换对为: 正变换 逆变换 其中 利用MATLAB 实现傅立叶级数计算
编写函数实现DFS 计算
function xk=dfs(xn,N) n=[0:1:N-1]; k=n;
WN=exp(-j*2*pi/N);
21100()[()]()()---=====∑∑
N N j nk
nk N N n n X k DFS x n x n e x n W π21100
11()[()]()()---=====∑
∑
N N j nk nk
N N n k x n IDFS X k X k e X k W N N π2j N N W e π-=()x n
nk=n'*k; WNnk=WN.^nk; xk=xn* WNnk;
例：xn=[0,1,2,3]，N=4
xn=[0,1,2,3]; N=4; xk=dfs(xn,N)'
逆运算IDFS
function xn=idfs(xk,N) n=[0:1:N-1]; k=n;
WN=exp(-j*2*pi/N); nk=n'*k; WNnk=WN.^(-nk); xn=xk*WNnk/N;
离散傅立叶变换的正、逆变换定义为：
比较正、逆变换的定义式可以看出，只要把DFT 公式中的系数 改为 ，并最后乘以1/N ，那么，DFT 的计算程序就可以用来计算IDFT 。
例：已知序列 试绘制序
列及其傅立叶变换幅度谱
N=100; n=0:N-1;
xn=cos(0.48*pi*n)+cos(0.52*pi*n); xk=dfs(xn,N); magxk=abs(xk); subplot(2,1,1) plot(n,xn) subplot(2,1,2) k=0:length(magxk)-1; plot(k,magxk)
210()[()]()N j nk N
n X k DFT x n x n e
π--===∑
2101()[()]()N j kn N k x n IDFT x k X k e N π-===∑
x(n)cos(0.48n)cos(0.52n),(0n 100)=π+π ≤≤2j nk N e π-2j nk N
e π
DFT 的应用
DFT 在数字滤波、功率谱分析、仿真、系统分析、通讯理论方面有广泛的应用。 DFT 的特性
1. 周期性
2. 对称性
3. 线性
4. 时移
5. 频移
6. 共轭
7. 折叠
8. 实序列的对称性 9. 卷积
利用MATLAB 对DFT 的特性进行验证
例： 分析：因为x(n)是复指数，它满足周期性，我们将在两个周期中的401个频点上作计算来观察其周期性。
n=0:10;
x=(0.9*exp(j*pi/3)).^n; k=-200:200; w=(pi/100)*k;
X=x*(exp(-j*pi/100)).^(n'*k); magX=abs(X); angX=angle(X); subplot(2,1,1); plot(w/pi,magX); subplot(2,1,2); plot(w/pi,angX/pi);
()(0.9exp(/3)),010
n x n j n π=≤≤
检验频移特性
乘以复数指数对应于一个频移
令
n=0:100; x=cos(pi*n/2); k=-100:100; w=(pi/100)*k; X=x*(exp(-j*pi/100)).^(n'*k); y=exp(j*pi*n/4).*x;
Y=y*(exp(-j*pi/100)).^(n'*k); subplot(2,2,1);plot(w/pi,abs(X)); axis([-1,1,0,60]);
subplot(2,2,2);plot(w/pi,angle(X)/pi); axis([-1,1,-1,1]);
subplot(2,2,3);plot(w/pi,abs(Y)); axis([-1,1,0,60]);
subplot(2,2,4);plot(w/pi,angle(Y)/pi); axis([-1,1,-1,1]);
从图中可以看出幅值和相位均沿频率轴平移了
()cos(/2),0100x n n n π= ≤≤/4()()
j n y n e x n π=4
π
从差分方程求频率响应当LTI 系统用差分方程表示如下：
上式做变换 消去共有项 得 例：一个LTI 系统的差分方程如下： y(n)=0.8y(n-1)+x(n) 求H(ejw) ，求出并画出它对输入 的稳态响应 把差分方程改写成
y(n)-0.8y(n-1)=x(n) 利用上面分析的公式，可得
将系统的输入x(n)带入 因此 输出端信号放大4.0928倍并移位3.42个采样周期
函数filter
对给定输入和差分方程系数时求解差分方程的数值解。 格式
y=filter(b,a,x)
其中b,a 为差分方程的系数向量，x 是输入序列。 输出y 和输入x 的长度一致。
b=1; a=[1,-0.8]; n=0:100; x=cos(0.05*pi*n); y=filter(b,a,x);
subplot(2,1,1); stem(n,x)
xlabel('n');ylabel('x(n)');title('输入序列'); subplot(2,1,2); stem(n,y)
xlabel('n');ylabel('y(n)');title('输出序列');
10
()()()
N M
l m l m y n a y n l b x n m ==+-=-∑∑
()()10
()()N M j j n j j n l j n m l m l m H e e a H e e b e ωωωωω--==+=∑∑
j n e ω0
1
()1M
j m m j m N j l
l l b e H e a e ωωω-
=-==+∑
∑
()cos(0.05)()x n n u n π=1()10.8j j H e e ωω-=
-0.050.5377
0.051() 4.092810.8j j j H e e
e ππ
-==-() 4.0928cos(0.050.5377) 4.0928cos[0.05( 3.42)]
rs y n n n ππ=-=-
四、实验报告要求
1.简述实验目的和实验原理。
2.总结实验中的主要结论，你的收获和体会。
实验三快速傅立叶变换(FFT)及其应用
一、实验目的
1.了解计算DFT算法存在的问题及改进途径。
2.掌握几种DFT算法（时间抽取算法DIT算法，频率抽取算法DIF算法，线性调频Z
变换即CZT法）。
3.学习并掌握FFT的应用。
二、实验原理
有限长序列通过离散傅里叶变换(DFT)将其频域离散化成有限长序列。但其计算量太大(与N的平方成正比）,，很难实时地处理问题,，因此引出了快速傅里叶变换(FFT)。
FFT并不是一种新的变换形式，它只是DFT的一种快速算法。并且根据对序列分解与选取方法的不同而产生了FFT的多种算法。
DFT的快速算法—FFT是数字信号处理的基本方法和基本技术，是必须牢牢掌握的。
时间抽选FFT算法的理论推导和流图详见《数字信号处理》教材。该算法遵循两条准则：（1）对时间奇偶分；（2）对频率前后分。
这种算法的流图特点是：
1)基本运算单元都是蝶形
任何一个长度为N=2M的序列，总可通过M次分解最后成为2点的DFT计算。如图所示：
?WNk称为旋转因子
?计算方程如下：
Xm+1(p)=Xm(p)+WNkXm(q)
Xm+1(q)=Xm(p)-WNkXm(q)
2)同址（原位）计算
这是由蝶形运算带来的好处，每一级蝶形运算的结果Xm+1(p)无须另外存储，只要再存入Xm(p)中即可，Xm+1(q)亦然。这样将大大节省存储单元。
3)变址计算
输入为“混序”（码位倒置）排列，输出按自然序排列，因而对输入要进行“变址”计算（即码位倒置计算）。“变址”实际上是一种“整序”的行为，目的是保证“同址”。FFT的应用
凡是利用付里叶变换来进行分析、综合、变换的地方，都可以利用FFT算法来减少其计算量。
FFT 主要应用在 1. 快速卷积 2. 快速相关 3. 频谱分析
快速傅立叶变换的MATLAB 实现
提供fft 函数计算DFT 格式
X=fft(x) X=fft(x,N)
如果x 的长度小于N ，则在其后填零使其成为N 点序列，若省略变量N ，则DFT 的长度即为x 的长度。
如果N 为2的幂，则得到高速的基-2FFT 算法；若N 不是2的乘方，则为较慢的混合算法。
如果x 是矩阵，则X 是对矩阵的每一列向量作FFT 。
例：已知信号由15Hz 幅值0.5的正弦信号和40Hz 幅值2的正弦信号组成，数据采样频率为100Hz,试绘制N=128点DFT 的幅频图。 由题目可得
fs=100; N=128; n=0:N-1; t=n/fs;
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); y=fft(x,N);
f=(0:length(y)-1)'*fs/length(y); mag=abs(y); stem(f,mag); title('N=128点')
利用FFT 进行功率谱的噪声分析
? 已知带有测量噪声信号 其中f1=50Hz,f2=120Hz, 为均值为零、方差为1的随机信号，采样频率为1000Hz ，
数据点数N=512。试绘制信号的频谱图和功率谱图。
12()sin(2)sin(2)2()x t f t f t
t ππω=++()t ω