信号处理.docx
- 文档编号:11219313
- 上传时间:2023-02-25
- 格式:DOCX
- 页数:45
- 大小:302.01KB
信号处理.docx
《信号处理.docx》由会员分享,可在线阅读,更多相关《信号处理.docx(45页珍藏版)》请在冰豆网上搜索。
信号处理
第九章信号处理
9.1信号的时域分析
9.1.1信号的表示
1.连续时间信号
对于连续时间信号来说,在MATLAB中可以有两种方法来表示。
一种是用符号运算的方式,一种是用数值计算的方式。
(1)用符号对象来表示连续时间信号
【例9.1】产生并显示几种常用信号
symst;
y1=2*sin(pi/4*t+pi/2);%产生正弦信号
y2=exp(3*t);%产生指数信号
y3=sin(t)/t;%产生抽样信号
subplot(311),ezplot(y1);
subplot(312),ezplot(y2);
subplot(313),ezplot(y3)
可以看出,这种表示信号的方法精度很高,但运算和显示速度较慢。
而且对于不能用表达式表示的信号形式无法使用。
MATLAB信号处理工具箱中的大部分函数都不支持符号运算,只有少量几个函数可以使用,但也不一定能够画图。
【例9.2】产生单位冲激函数和单位阶跃函数
symst;
y1=dirac(t);
y2=heaviside(t);
这里由于dirac(0)的值是inf,所以无法显示;而单位阶跃函数中在t=0处无定义。
(2)用向量来表示连续时间信号
可以用连续时间信号的一系列采样值来近似的表示信号,因此,需要采样间隔要足够小。
这时,需要有两个向量来共同描述该信号,一个描述各采样时间,一个描述各采样点的值。
信号的显示可以用plot命令来完成,但这种方法只能描述有限长度的信号。
【例9.2】用向量表示正弦信号和指数信号
dt=0.01;
t=-2*pi:
dt:
2*pi;
y1=2*sin(pi/4*t+pi/2);%产生正弦信号
y2=exp(3*t);%产生指数信号
subplot(211),plot(t,y1),title('2*sin(pi/4*t+pi/2)')
subplot(212),plot(t,y2),title('exp(3*t)');
【例9.3】用向量表示抽样函数sin(t)/t
分析:
由于当t=0时MATLAB会计算出该点的值为NaN,所以当时间向量中有0值时,不能直接使用语句sin(t)./t。
这时可以使用sinc函数来完成。
MATLAB中sinc函数的计算规则是
,显然,只需将时间向量先缩小
倍即可表示抽样函数。
dt=0.01;
t=-2*pi:
dt:
2*pi;
y1=sin(t)./t;%这里t向量中无0值
y2=sinc(t/pi);
subplot(211),plot(t,y1);
subplot(212),plot(t,y2);
【例9.4】单位冲激信号和单位阶跃信号
分析:
描述单位冲激信号时,时间向量中必须有0。
可以使用dirac函数来产生单位冲激信号。
虽然dirac(0)的值是inf而无法显示,但可将信号在t=0处的值指定为1/dt近似地表示单位冲激信号。
可以用heaviside函数来产生单位阶跃信号,也可用1/2+sign(t)/2或其他办法。
dt=0.001;
t=-2:
dt:
2;
y1=dirac(t);
y1(find(y1>0))=1/dt;
y2=heaviside(t);
subplot(211),plot(t,y1);
subplot(212),plot(t,y2);
axis([-3302])
信号处理工具箱中的信号产生函数
函数名
功能
函数名
功能
sawtooth
产生锯齿波或三角波信号
pulstran
产生冲激串
square
产生方波信号
rectpuls
产生非周期的方波信号
sinc
产生sinc函数波形
tripuls
产生非周期的三角波信号
chirp
产生调频余弦信号
diric
产生Dirichlet或周期sinc函数
gauspuls
产生高斯正弦脉冲信号
gmonopuls
产生高斯单脉冲信号
vco
电压控制振荡器
(2)离散时间信号
和连续时间信号不同,离散时间信号不能用符号运算来表示,且时间向量必须为整数。
画图时可以用stem函数来画出离散数据。
注:
离散时间信号和用向量描述的连续时间信号相比,唯一的区别是时间向量,也就是【n】和【n*dt】的区别。
【例9.5】单位冲激序列和单位阶跃序列
n=-5:
5;
y1=zeros(1,length(n));
y1(n==0)=1;
y2=zeros(1,length(n));
y2(n>=0)=1;
subplot(211),stem(n,y1,'.'),ylim([0,2]),title('单位冲激序列');
subplot(212),stem(n,y2,'.'),ylim([0,2]),title('单位阶跃序列');
【例9.6】实指数序列2n与2nU(n)
n=-5:
5;
y1=2.^n;
y2=y1;
y2(n<0)=0;
subplot(211),stem(n,y1,'.');
subplot(212),stem(n,y2,'.');
9.1.2、信号的时域基本运算
1符号对象信号的基本运算
对于用符号对象表示的连续时间信号来说,其各种运算比较简单,均可按照符号运算直接完成。
【例9.7】符号对象信号的基本运算
symst;
x1=sin(t);
x2=cos(t);
y1=x1+x2;%相加
y2=x1*x2;%相乘
y3=subs(x1,t,t-2);%移位
y4=subs(x1,t,-t);%翻折
y5=subs(x1,t,2*t);%尺度变换
y6=-x1;%倒相
注意:
对于连续时间信号的卷积积分
,虽然可以用符号运算“int(subs(x,t,θ)*subs(h,t,t-θ),θ,-inf,inf)”来表示,但式中存在参数t,因此MATLAB往往不能计算出精确结果,而是给出警告提示并返回该函数的原表达式。
【例9.8】符号运算完成卷积积分的运算
symsttheta;
x1=heaviside(t)-heaviside(t-1);
x2=heaviside(t-1)-heaviside(t-3);
y=int(subs(x1,t,theta)*subs(x2,t,t-theta),theta,-inf,inf);
此时,可以通过代入时间向量,将符号运算的解转为数值解,从而形成以向量表示的卷积结果。
如上例,加入以下代码:
tt=1:
0.01:
4;
fork=1:
length(tt)
yy(k)=double(subs(y,t,tt(k)));
end
plot(tt,yy)
2向量信号的基本运算
(1)常用运算
对于向量表示的信号,不论是连续时间信号还是离散时间信号,其各种基本运算都是各数据点按照数值计算的方式来完成的。
对于相加或相乘运算,关键是要判断信号各数据点的位置和长度。
运算时必须通过补零的方法使两信号的各样本点在时间上完全对应,且长度必须相同。
y(n)=x1(n)+x2(n)MATLAB实现:
y=x1+x2
y(n)=x1(n)x2(n)MATLAB实现:
y=x1.*x2
其他基本运算:
y(n)=x(n-k)MATLAB实现:
ny=nx+k
y(n)=x(-n)MATLAB实现:
ny=fliplr(-nx),y=fliplr(x)
y(n)=-x(n)MATLAB实现:
y=-x
(2)卷积
对于离散时间信号来说,卷积运算表现为卷积和的形式:
对此可以直接调用conv函数,但需考虑两序列卷积后的结果位置。
对于连续时间信号的卷积积分,由于
因此
显然,可以用conv函数计算出两信号的离散卷积,再乘以采样间隔,就可以得到卷积积分的采样值。
【例9.9】编写卷积积分的计算函数
function[f,t]=ctsconv(f1,f2,t1,t2,Ts)
%计算连续时间信号的卷积积分:
f=f1*f2
%f1、f2:
输入信号的样值向量
%t1、t2:
输入信号的时间向量
%Ts:
采样间隔
%t1、t2的采样间隔必须等于Ts
%--------------------------------------------------------------------
f=conv(f1,f2);%计算序列f1与f2的卷积和f
f=f*Ts;%计算卷积积分信号f(t)离散样值
t0=t1
(1)+t2
(1);%计算序列f非零样值的起点位置
l=length(t1)+length(t2)-2;%计算卷积积分f的非零样值的宽度
t=t0:
Ts:
(t0+l*Ts)%确定卷积积分f非零样值的时间向量
【例9.10】已知x(t)=u(t)-u(t-1),h(t)=u(t-1)-u(t-3),求y=x(t)*h(t)
clear;
p=0.01;
t1=0:
p:
1;
f1=ones(size(t1));
t2=1:
p:
3;
f2=ones(size(t2));
[f,t]=ctsconv(f1,f2,t1,t2,p);
subplot(2,2,1)
plot(t1,f1)
title('x(t)')
xlabel('t1')
subplot(2,2,2)
plot(t2,f2)
title('h(t)')
xlabel('t2')
subplot(2,2,3)
plot(t,f);
p=get(gca,'position');
p(3)=2.5*p(3);%p(3)代表坐标轴的宽度(width)set(gca,'position',p)%将第三个子图的横坐标范围扩为原来的2.5倍
title('y(t)=x(t)*h(t)')
xlabel('t')
执行后即可实现连续信号的卷积积分运算。
说明:
gca%获取当前坐标轴句柄
创建坐标轴句柄:
h_axes=axes(‘position’,[left,bottom,width,height])定义轴的位置和大小
(3)序列的循环移位
根据循环移位的表达式:
可以看出,循环移位的实现只需对时间向量进行移位在求余即可。
【例9.11】编写循环移位函数
functiony=CircularShift(x,m,N)
%y=CircularShift(x,m,N)
%yistheputputsequenceafterthecircularshift
%xistheinputsequence
%misthesamplesoftheshift
%Nisthesizeofcircularbuffer
iflength(x)>N
error('Nmustbelargerorequaltothelengthofx')
end
x=[x,zeros(1,N-length(x))];
n=0:
N-1;
n=mod(n-m,N);
y=x(n+1);
(4)有限长序列的循环卷积
根据循环卷积的计算公式:
可以写出矩阵形式:
显然,中间的卷积矩阵可以用循环移位来完成。
【例9.12】编写循环卷积函数
functiony=CircularConv(x,h,N)
%yistheoutputsequenceofthecircularconvolution
%x,haretheinputsequencesindividually
%Nisthesizeofcircularconvolutionbuffer
x=[x,zeros(1,N-length(x))];
h=[h,zeros(1,N-length(h))];
k=0:
N-1;
h=h(mod(-k,N)+1);
Z=zeros(N,N);
forn=1:
N
Z(n,:
)=CircularShift(h,n-1,N);
end
y=x*Z';
实现循环卷积的方法有多种,除了上述以外,还可以用循环卷积和DFT的关系来完成:
也可以利用线性卷积和循环卷积的关系来完成,即有限长序列
和
的线性卷积
以L为周期的周期延拓序列的主值序列即为这两个有限长序列的循环卷积。
9.2信号的频域分析
9.2.1连续时间周期信号的傅里叶级数分解与合成
设周期为T的周期信号f(t)满足狄利赫利条件,则f(t)可以由三角函数的傅里叶级数线性组合表示如下:
其中:
n=1,2,3,….
n=1,2,3,….
可以看出,傅里叶级数系数的计算主要就是定积分的计算,因此可使用多种方法完成。
【例9.13】符号运算求解积分的方法求解傅里叶级数系数
function[an,bn]=fs(x,t,T,K)
%求解周期信号的傅里叶级数系数
%x:
周期信号在【-T/2,T/2】内的符号表达式
%t:
符号表达式中的自变量
%T:
周期信号的周期
%K:
傅里叶级数分解的谐波次数
%---------------------------------------------------------
symsn;
an=zeros(K+1,1);
bn=zeros(K+1,1);
A0=int(x,t,-T/2,T/2)/T;
An=2*int(x*cos(2*pi*n*t/T),t,-T/2,T/2)/T;
Bn=2*int(x*sin(2*pi*n*t/T),t,-T/2,T/2)/T;
an
(1)=double(A0);
fork=1:
K
an(k+1)=double(subs(An,n,k));
bn(k+1)=double(subs(Bn,n,k));
end
假设周期信号为周期为1,脉宽为0.3的周期矩形脉冲信号,进行3次傅里叶级数分解:
clear;
T=1;
tao=0.3;
K=3;
symst
%s1=strcat('heaviside(t+',num2str(tao/2),')');%构造u(t+tao/2)
%s2=strcat('heaviside(t-',num2str(tao/2),')');%构造u(t-tao/2)
%x=sym(s1)-sym(s2);%构造一个周期的脉冲信号u(t+tao/2)-u(t-tao/2)
x=heaviside(t+tao/2)-heaviside(t-tao/2);
[a,b]=fs(x,t,T,K)
继续输入以下代码,可以方便的图示出信号的分解与合成情况:
dt=T/1000;
tt=-T/2:
dt:
T/2-dt;
y=zeros(length(a)+1,length(tt));
y(1,:
)=double(subs(x,t,tt));
fork=1:
length(a)
y(k+1,:
)=a(k)*cos(2*pi*(k-1)*tt/T)+b(k)*sin(2*pi*(k-1)*tt/T);
end
yy=sum(y);
subplot(211),plot(tt,y),title('原信号与谐波分解')
subplot(212),plot(tt,yy),title('合成信号')
注:
strcat()和strvcat()两函数的区别在于,前者是将字符串沿横向连接成更长的字符串,而后者是将字符串沿纵向连接成二维字符数组。
【例9.14】梯形法求解积分的方法求解傅里叶级数系数
function[an,bn]=fs(x,t,T,K)
%求解周期信号的傅里叶级数系数
%x:
周期信号在一个周期内的样值向量
%t:
x所对应的的时间向量
%T:
周期信号的周期
%K:
傅里叶级数分解的谐波次数
%---------------------------------------------------------
an=zeros(K+1,1);
bn=zeros(K+1,1);
an
(1)=trapz(t,x)/T;
fork=1:
K
an(k+1)=trapz(t,x.*cos(2*pi*k*t/T))*2/T;
bn(k+1)=trapz(t,x.*sin(2*pi*k*t/T))*2/T;
end
【例9.15】矩形法求解积分的方法求解傅里叶级数系数
function[an,bn]=fs(x,t,T,K)
%求解周期信号的傅里叶级数系数
%x:
周期信号在一个周期内的样值向量
%t:
x所对应的的时间向量
%T:
周期信号的周期
%K:
傅里叶级数分解的谐波次数
%---------------------------------------------------------
an=zeros(K+1,1);
bn=zeros(K+1,1);
dt=T/length(t);
an
(1)=sum(x)*dt/T;
fork=1:
K
an(k+1)=sum(x.*cos(2*pi*k*t/T))*dt*2/T;
bn(k+1)=sum(x.*sin(2*pi*k*t/T))*dt*2/T;
end
此外,也可以利用函数句柄,并使用quad泛函命令来完成积分计算。
利用相同的方法,可以自行编写复指数形式的傅里叶级数的计算函数。
当然,也可以利用
的系数关系进行转换。
9.2.2连续时间非周期信号的傅里叶变换
对于符号对象描述的信号,fourier变换和反变换可以利用积分函数int来实现,也可以直接使用fourier或ifourier函数实现。
【例9.16】求解x(t)=u(t+1)-u(t-1)的傅里叶变换,并画出幅频特性。
symstw;
x=heaviside(t+1)-heaviside(t-1);
y=fourier(x,t,w)
ezplot(abs(y),[-10,10])
注意:
1.符号运算中没有求相角的函数!
2.如果积分时不能给出精确解,或者要求相频特性时,可以转化为向量数据,使用plot函数作图。
对于用向量描述的信号,其傅里叶变换
仍然是连续谱,因此只能对
轴进行采样来近似
,即
写成矩阵形式:
所以,可根据频率分辨率的要求确定频域采样间隔,然后通过矩阵运算来近似原信号的傅里叶变换。
【例9.17】求解x(t)=u(t+1)-u(t-1)的傅里叶变换,并画出幅频特性和相频特性。
dt=0.01;
t=-3:
dt:
3;
x=heaviside(t+1)-heaviside(t-1);
dw=0.001;
w=-10:
dw:
10;
f=x*exp(-j*t’*w)*dt;
subplot(311),plot(t,x),title(‘原信号’);
subplot(312),plot(w,abs(f)),title(‘幅度’);
subplot(313),plot(w,angle(f)),title(‘相位’)
9.2.3离散时间非周期信号的傅里叶变换
由于
是连续且以2π为周期,所以根据频率分辨率的要求对其进行采样:
写成矩阵形式:
可以看出,这里与连续时间傅里叶变换的区别有两点:
没有
的增益
时间向量是【n】,而不是【n*dt】
【例9.18】求解x(n)=u(n+1)-u(n-1)的傅里叶变换,并画出幅频特性和相频特性。
n=-3:
3;
x=zeros(1,length(n));
x((n>=-1)&(n<=1))=1;
dw=0.001;
w=-10:
dw:
10;
f=x*exp(-j*n’*w);
subplot(311),stem(n,x),title(‘原信号’);
subplot(312),plot(w,abs(f)),title(‘幅度’);
subplot(313),plot(w,angle(f)),title(‘相位’)
从结果中可以看出DTFT的周期性,因此大多数情况下只需描述出一个周期内的频率特性即可。
比如此例中将w的范围变为:
w=-pi:
dw:
pi;或者w=0:
dw:
2*pi;
9.2.4离散傅里叶变换DFT
可以看出,DFT的计算和DTFT的计算基本相同,只是一下区别:
(1)时间向量【n】总是从0开始,N结束;
(2)频率向量【k】总是从0开始,N结束,但实际描述的频率是【k*2pi/N】,所以实际频率范围是0~2*pi,这一点在谱分析时要注意;
(3)时间向量的长度和频率向量的长度均为N。
因此只需注意时间向量和频率向量的赋值即可,其他均与DTFT一样。
【例9.19】编写计算DFT的函数。
function[Xk]=dft(xn,N)
%计算DFT
%[Xk]=dft(xn,N)
%Xk是频域的DFT序列
%xn是时域的有限长序列
%N是DFT的长度
n=[0:
1:
N-1];
k=[0:
1:
N-1];
if(N>length(xn))
xn=[xn,zeros(1,N-length(xn))];
else
xn=xn(1:
N);
end
Xk=xn*exp(-j*2*pi/N*n’*k);
根据DFT反变换公式
可以看出IDFT的计算和DFT的计算基本相同,只需略加修改即可。
【例9.20】编写计算IDFT的函数。
function[xn]=idft(Xk,N)
%计算IDFT
%[xn]=idft(Xk,N)
%Xk是频域的DFT序列
%xn是时域的有限长序列
%N是IDFT的长度
n=[0:
1:
N-1];
k=[0:
1:
N-1];
if(N>length(Xk))
Xk=[Xk,zeros(1,N-length(Xk))];
else
Xk=Xk(1:
N);
end
xn=Xk*exp(j*2*pi/N*n’*k)/N;
其实,可以直接调用fft函数和ifft函数来完成DFT和IDFT的计算。
9.2.5离散傅里叶级数DFS
DFS的完成很简单,根据DFS和DFT的关系,只需对时域周期序列的主值序列进算DFT,然后对结果进行延拓即可。
9.2.6利用DFT(FFT)进行谱分析
1.对有限长的时域离散序列进行谱分析
根据DFT的物理意义,如果有限长序列的长度为M,则该序列的N(
)点DFT就是它的频谱函数在
区间上的N点等间隔采样。
显然,DFT的变换点数N越大,频域的抽样就越密集,其DFT结果就越能反映频谱曲线的变化情况。
但实际应用中不可能选择无限大的N,而是根据频率分辨率的要求确定N的值。
N点DFT的频域采样间隔为
,也就是说N点DFT能够实现的频率分辨率是
。
如果要求的频率分辨率是D弧度,则N必须满足
,即
。
【例9.21】假设
,试用DFT分析其频谱,要求频率分辨率为
。
分析:
根据频率分辨率要求,有
,即
,因此取
N=100;
n=0:
9;
xn=0.5.^n;
Xk=fft(xn,N);
k=0:
N-1;
wk=2*k/N;
subplot(1,2,1);plot(wk,abs(Xk));xlabel('ω/π');ylabel('幅度
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 信号 处理