秋季数字信号处理实验讲义修改版.docx
- 文档编号:24458788
- 上传时间:2023-05-27
- 格式:DOCX
- 页数:25
- 大小:111.05KB
秋季数字信号处理实验讲义修改版.docx
《秋季数字信号处理实验讲义修改版.docx》由会员分享,可在线阅读,更多相关《秋季数字信号处理实验讲义修改版.docx(25页珍藏版)》请在冰豆网上搜索。
秋季数字信号处理实验讲义修改版
实验一序列的产生及绘图
一、实验目的
1.熟悉信号处理软件MATLAB的使用。
2.离散信号的基本运算实现。
3.了解基本序列及复杂序列的产生方法。
4.运用卷积方法观察系统的时域特性。
5.掌握线性时不变系统的频域表示方法。
二、实验内容
1.熟悉扩展函数
2.运行例题程序
3.编程实现下列内容
(1)利用扩展函数产生序列并画图
(a)
-5<=n<=5
(b)
和
0<=n<=50
w(n)为白噪声函数为w=randn(size(n))
(2)设线性移不变系统的抽样响应为
输入序列为
求系统输出y(n)并画图
提示:
输出为输入和抽样响应的卷积
三、实验报告要求
1.记录例题程序的实验结果、图形。
2.写出自己编写的程序并记录结果、图形。
注:
以下程序中所有以%开头的行均为注释,所有汉字均为注释,%后的内容不用写入程序
%如果要了解哪个函数的应用方法请用help命令如helpzreos
%本软件中*表示乘法,卷积用函数conv或修改后的卷积conv_m
%以下是7个扩展函数
%扩展函数1~7的用法和该软件自带函数用法一致,即在调用时要将实参代入
%例:
应用扩展函数3需要输入x1(n),x2(n)的值。
%在CommandWindow(命令窗口)中输入
%n1=1:
5;
%n2=2:
6;
%x1=[13579];
%x2=[246810];
%[y,n]=sigadd(x1,n1,x2,n2)即可得两序列相加的结果
%7个扩展函数要分别存到不同的文件中,并且文件名要和该扩展函数的函数名一致
%如产生单位取样序列的函数所存文件的文件名必须为impseq
%1.单位取样序列x(n)=delta(n-n0)要求n1<=n0<=n2
function[x,n]=impseq(n0,n1,n2)
n=[n1:
n2];
x=[(n-n0)==0];
%2.单位阶跃序列x(n)=u(n-n0)要求n1<=n0<=n2
function[x,n]=stepseq(n0,n1,n2)
n=[n1:
n2];
x=[(n-n0)>=0];
%3.信号加y(n)=x1(n)+x2(n)
%find函数:
找出非零元素的索引号
%x1:
第一个序列的值,n1:
序列x1的索引号
%x2:
第二个序列的值,n2:
序列x2的索引号
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;
%4.信号乘y(n)=x1(n)*x2(n)
function[y,n]=sigmult(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;
%5.移位y(n)=x(n-n0)
function[y,n]=sigshift(x,m,n0)
n=m+n0;y=x;
%6.翻褶y(n)=x(-n)
function[y,n]=sigfold(x,n)
y=fliplr(x);n=-fliplr(n);
%7.修改后的卷积(系统自带卷积函数(conv)只能得到输出值,而无法表示取值范围(n),修改后的卷
%积既给出了卷积值,也给出了它的取值范围.)
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序列的基本运算
nx1=1:
5;
nx2=2:
6;
x1=[13579];
x2=[246810];
[y1,n1]=sigadd(x1,nx1,x2,nx2)%序列相加
[y2,n2]=sigmult(x1,nx1,x2,nx2)%序列相乘
n0=3;
[y3,n3]=sigshift(x1,nx1,n0)%序列移位,n0为移动的位数
[y4,n4]=sigfold(x2,nx2)%序列的翻褶
y5=sum(x1)%序列和
y6=prod(x1)%序列积
ex=sum(x1.*conj(x1))%或ex=sum(abs(x).^2);
%信号能量
%例2产生序列并画图
%
0<=n<=20
%subplot函数:
将一个图形分为几个小区域,每次选择一个作为当前区域画图
%subplot(m,n,p)将图形分为m行,n列的区域,p为当前区域
%axis函数:
控制坐标系的刻度和形式title函数:
图形标题
%xlabel函数:
X轴标记ylabel函数:
Y轴标记
%stem函数:
画离散序列图
n=[0:
20];
x1=n.*(stepseq(0,0,20)-stepseq(10,0,20));
x2=10*exp(-0.3*(n-10)).*(stepseq(10,0,20)-stepseq(20,0,20));
x=x1+x2;
subplot(2,1,1);stem(n,x);
xlabel('n');ylabel('x(n)');axis([0,20,-1,11]);
%例3产生复信号
-10<=n<=10
%并画出复序列的实部、虚部、幅值和相位离散图
figure
(2);clf
n=[-10:
10];alpha=-0.1+0.3j;x=exp(alpha*n);
subplot(2,2,1);stem(n,real(x));title('real');xlabel('n')
subplot(2,2,2);stem(n,imag(x));title('imag');xlabel('n')
subplot(2,2,3);stem(n,abs(x));title('magtitude');
xlabel('n')
subplot(2,2,4);stem(n,(180/pi)*angle(x));title('phase');
xlabel('n')
%例4线性时不变系统的频域表示
b=[0.20.1];
a=[1-0.4-0.5];%系统函数的系数
h=impz(b,a);%系统的单位取样响应
figure
(1)
stem(h)%画出单位取样响应
title('h(n)')
figure
(2)
fs=1000;
[H,f]=freqz(b,a,256,fs);%求出系统的频率响应
mag=abs(H);%幅度响应
ph=angle(H);%相位响应
ph=ph*180/pi;
subplot(2,1,1),plot(f,mag);grid%画出幅度响应
xlabel('frequency(Hz)');
ylabel('magnitude');
subplot(2,1,2);plot(f,ph);grid%画出相位响应
xlabel('frequency(Hz)');
ylabel('phase');
figure(3)
zr=roots(b)%求出系统的零点
pk=roots(a)%求出系统的极点
zplane(b,a);%zplane函数画出零极点图
实验二离散傅立叶变换及谱分析
一、实验目的
1.掌握离散傅里叶变换的计算机实现方法。
2.检验实序列傅里叶变换的性质。
3.掌握计算序列的圆周卷积的方法。
4.熟悉连续信号经理想采样前后的频谱变化关系,加深对时域采样定理的理解。
5.学习用DFT对连续信号和时域离散信号进行谱分析的方法,了解可能出现的分析误差,以便在实际中正确应用DFT。
二、实验内容
1.实现离散傅里叶变换。
2.计算序列圆周卷积。
3.计算实序列傅里叶变换并检验DFT性质。
4.实现连续信号傅里叶变换以及由不同采样频率采样得到的离散信号的傅里叶变换。
5.实现补零序列的傅里叶变换。
6.实现高密度谱和高分辨率谱,并比较二者的不同。
三、实验报告要求
见各程序要求
%以下为4个扩展函数
%
(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;
%
(2)逆离散傅立叶变换
function[xn]=idft(Xk,N)
n=[0:
1:
N-1];
k=[0:
1:
N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(-nk);
xn=(Xk*WNnk)/N;
%(3)实序列的分解
%实序列可分解为共扼对称分量
%和共扼反对称分量
function[xec,xoc]=circevod(x)
N=length(x);
n=0:
(N-1);
xec=0.5*(x+x(mod(-n,N)+1));%根据上面的公式计算,mod函数为取余
xoc=0.5*(x-x(mod(-n,N)+1));
%(4)序列的循环移位
functiony=cirshftt(x,m,N)
iflength(x)>N
error('Nmustbe>=thelengthofx')%要求移位周期大于信号长度
end
x=[xzeros(1,N-length(x))];
n=[0:
1:
N-1];
n=mod(n-m,N);
y=x(n+1);
%例1本例检验实序列的性质DFT[xec(n)]=Re[X(k)]DFT[xoc(n)]=Im[X(k)]
%设x(n)=10*(0.8).^n0<=n<=10将x(n)分解为共扼对称及共扼反对称部分
%实验报告要求:
(1)将实验结果画出
(2)实验结果说明什么
n=0:
10;
x=10*(0.8).^n;
[xec,xoc]=circevod(x);
subplot(2,1,1);stem(n,xec);%画出序列的共扼对称分量
title('Circular-evencomponent')
xlabel('n');ylabel('xec(n)');axis([-0.5,10.5,-1,11])
subplot(2,1,2);stem(n,xoc);%画出序列的共扼反对称分量
title('Circular-oddcomponent')
xlabel('n');ylabel('xoc(n)');axis([-0.5,10.5,-4,4])
figure
(2)
X=dft(x,11);%求出序列的DFT
Xec=dft(xec,11);%求序列的共扼对称分量的DFT
Xoc=dft(xoc,11);%求序列的共扼反对称分量的DFT
subplot(2,2,1);stem(n,real(X));axis([-0.5,10.5,-5,50])
title('Real{DFT[x(n)]}');xlabel('k');%画出序列DFT的实部
subplot(2,2,2);stem(n,imag(X));axis([-0.5,10.5,-20,20])
title('Imag{DFT[x(n)]}');xlabel('k');%画出序列DFT的虚部
subplot(2,2,3);stem(n,Xec);axis([-0.5,10.5,-5,50])
title('DFT[xec(n)]');xlabel('k');
subplot(2,2,4);stem(n,imag(Xoc));axis([-0.5,10.5,-20,20])
title('DFT[xoc(n)]');xlabel('k');
%例2本例为计算序列的圆周卷积程序
%运行之前应在命令窗口输入x1,x2,N的值
%实验报告要求:
自己选择2个序列进行计算,将实验结果写出
iflength(x1)>N
error('Nmustbe>=thelengthofx1')
end
iflength(x2)>N
error('Nmustbe>=thelengthofx2')
end
x1=[x1zeros(1,N-length(x1))];%将x1,x2补0成为N长序列
x2=[x2zeros(1,N-length(x2))];
m=[0:
1:
N-1];
x2=x2(mod(-m,N)+1);%该语句的功能是将序列翻褶,延拓,取主值序列
H=zeros(N,N);
forn=1:
1:
N%该for循环的功能是得到x2序列的循环移位矩阵
H(n,:
)=cirshftt(x2,n-1,N);%和我们手工计算圆周卷积得到的表是一致的
end
y=x1*H'%用矩阵相乘的方法得到结果
%例3本例验证采样定理
%令
,绘制其傅立叶变换
。
用不同频率对其进行采样,分别画出离
%散时间傅立叶变换
。
已给出采样频率为
时的的程序
%实验报告要求:
(1)请写出
时的程序
(2)将实验结果画出
(3)实验结果说明什么(采样频率不同结果有何不同)。
Dt=0.00005;%步长为0.00005s
t=-0.005:
Dt:
0.005;
xa=exp(-1000*abs(t));%取时间从-0.005s到0.005s这段模拟信号
Wmax=2*pi*2000;%信号最高频率为2
*2000
K=500;%频域正半轴取500个点进行计算
k=0:
1:
K;
W=k*Wmax/K;%
求模拟角频率
Xa=xa*exp(-j*t'*W)*Dt;%计算连续时间傅立叶变换(利用矩阵运算实现)
Xa=real(Xa);%取实部
W=[-fliplr(W),W(2:
501)];%将角频率范围扩展为从-到+
Xa=[fliplr(Xa),Xa(2:
501)];
subplot(2,2,1);
plot(t*1000,xa);%画出模拟信号,横坐标为时间(毫秒),纵坐标为幅度
xlabel('time(millisecond)');ylabel('xa(t)');
title('anologsignal');
subplot(2,2,2);
plot(W/(2*pi*1000),Xa*1000);%画出连续时间傅立叶变换
xlabel('frequency(kHZ)');%横坐标为频率(kHz)
ylabel('xa(jw)');%纵坐标为幅度
title('FT');
%下面为采样频率5kHz时的程序
Ts=0.0002;%采样间隔为
n=-25:
1:
25;
x=exp(-1000*abs(n*Ts));%离散时间信号
K=500;k=0:
1:
K;w=pi*k/K;%w为数字频率
X=x*exp(-j*n'*w);%计算离散时间傅立叶变换(序列的傅立叶变换)
X=real(X);
w=[-fliplr(w),w(2:
K+1)];
X=[fliplr(X),X(2:
K+1)];
subplot(2,2,3);
stem(n*Ts*1000,x);%画出采样信号(离散时间信号)
xlabel('time(millisecond)');
gtext('Ts=0.2ms');%该语句可以将引号中的内容放置在figure中的任何地方,只需
%将十字的中心放在想放置内容的地方,然后按鼠标即可。
ylabel('x1(n)');
title('discretesignal');
subplot(2,2,4);
plot(w/pi,X);%画出离散时间傅立叶变换
xlabel('frequency(radian)');%横坐标为弧度
ylabel('x1(jw)');title('DTFT');
%例4本例说明补零序列的离散傅立叶变换
%序列
,已给出序列的傅立叶变换程序和将原序列补零到10长序列的DFT
%实验报告要求:
(1)编写将序列补零到20长后计算DFT的程序
(2)将实验结果画出
(3)实验结果说明什么(即序列补零后进行DFT,频谱有何变化)
n=0:
4;
x=[ones(1,5)];%产生矩形序列
k=0:
999;w=(pi/500)*k;
X=x*(exp(-j*pi/500)).^(n'*k);%计算离散时间傅立叶变换
Xe=abs(X);%取模
subplot(3,2,1);stem(n,x);ylabel('x(n)');%画出矩形序列
subplot(3,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|');%画出离散时间傅立叶变换
N=10;x=[ones(1,5),zeros(1,N-5)];%将原序列补零为10长序列
n=0:
1:
N-1;
X=dft(x,N);%进行DFT
magX=abs(X);
k=(0:
length(magX)'-1)*N/length(magX);
subplot(3,2,3);stem(n,x);ylabel('x(n)');%画出补零序列
subplot(3,2,4);stem(k,magX);%画出DFT结果
axis([0,10,0,5]);ylabel('|X(k)|');
%例5本题说明高密度谱和高分辨率谱之间的区别,高密度谱是信号补零后得到的,虽然谱线相当
%密但是因为信号有效长度不变,所以其分辨率也不变,因此还是很难看出信号的频谱成分。
高分辨
%率谱是将信号有效长度加长,因此分辨率提高,可以看出信号的成分。
%有一个序列为
(该序列周期计算可得40)
%
(1)下面给出有10个有效采样点序列的DFT程序
%
(2)请写出将第一问中的10长序列补零到40长,计算其DFT
%(3)采样n=0:
39,计算有40个有效采样点的序列的DFT
%实验报告要求:
(1)请编写将有10个有效采样点的序列补零到40长后计算DFT的程序
(2)请编写计算有40个有效采样点的序列的DFT的程序
(3)将实验结果画出并分析实验结果说明什么
M=10;
n=0:
M-1;
x=2*cos(0.35*pi*n)+cos(0.5*pi*n);
subplot(2,1,1);stem(n,x);title('没有足够采样点的信号');
Y=dft(x,M);
k1=0:
M-1;w1=2*pi/M*k1;
subplot(2,1,2);stem(w1/pi,abs(Y));title('信号的频谱');
实验三快速傅立叶变换
一、实验目的
1.掌握fft函数的用法。
2.了解FFT的用途。
二、实验内容
1.利用FFT进行信号检测。
分析信号频谱所对应频率轴的数字频率和频率轴的频率范围。
2.对例2,进一步增加截取长度和FFT点数,如N加大到256,观察信号频谱的变化,分析产生这一变化的原因。
在截取长度不变的条件下改变采样频率,观察信号频谱的变化,分析产生这一变化的原因。
3.对例3,加大噪声到2*randn(1,N)和8*randn(1,N),画出并比较不同噪声下时域波形和频谱。
4.比较DFT和FFT的运算时间。
(计时函数tic,toc)
5.对给定语音信号进行谱分析,写出采样频率,画出语音信号的波形及频谱,并分析语音信号的频率分布特点。
三、实验报告要求
1.记录例题程序的实验结果、图形。
2.写出自己编写的程序并记录结果、图形。
3.对实验结果进行分析。
%fft一维快速傅立叶变换函数
%格式:
y=fft(x)FFT算法计算矢量x的离散傅立叶变换,当x为矩阵时,y为矩阵x的每一列
%的FFT。
%当x的长度为2的幂次方时,则fft采用基2的FFT算法,否则采用稍慢的混合基算法。
%y=fft(x,n)实现n点FFT。
当x的长度小于n时,fft函数在x的尾部补零,以构成n点数据;%当x的长度大于n时,fft函数会截断序列x。
当x为矩阵时,fft函数按类似的方法处理列长度。
%ifft一维快速傅立叶反变换
%格式:
y=ifft(x)用于计算矢量x的IFFT。
%y=ifft(x,n)计算n点IFFT。
%例1求信号功率谱
t=0:
0.001:
0.6;
x=sin(2*pi*50*t)+sin(2*pi*120*t)+randn(1,length(t));
Y=fft(x,512);
P=Y.*conj(Y)/512;
f=1000*(0:
255)/512;
plot(f,P(1:
256))
%例2模拟信号
,以
进行取样,求N点DFT的幅值谱
%N分别为
(1)N=45
(2)N=50(3)N=55(4)N=60。
subplot(2,2,1)
N=45;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
plot(q,abs(y));title('FFTN=45')
subplot(2,2,2)
N=50;n=0:
N-1;t=0.01*n;q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
plot(q,abs(y));title('FFTN=50')
subplot(2,2,3)
N=55;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
plot(q,abs(y));title('FFTN=55')
subplot(2,2,4)
N=60;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);
y=fft(x,N);
plot(q,abs(y));title('FFTN=60')
%例3参数同上,N取64,并在信号中加入噪声w(t),比较有无噪声时的信号谱
%由结果可以看出这种噪声不影响信号检测
figure
(2)
subplot(2,1,1)
N=64;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin(4*pi*t)+5*cos(8*pi*t);y=fft(x,N);
plot(q,abs(y));title('FFTN=64')
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 秋季 数字信号 处理 实验 讲义 修改