Matlab基础2015-3

——离散时间信号与系统MATLAB基础课程安排MATLAB简介连续时间信号与系统离散时间信号与系统MATLAB实验内容离散时间系统的时域分析z变换、z域分析离散时间系统的频响特性内容离散时间系统的时域分析1单位样值(单位冲激)响应常系数线性差分方程的求解卷积、解卷积常系数线性差分方程的求解•MATLAB提供数值解法（即迭代法）计算差分方程的完全解。即给定传递函数、激励序列和边界条件后，用filter函数得到输出序列•[y,wf]=filter(b,a,x,wi),wi和wf分别表示系统的初始状态和终止状态常系数线性差分方程的求解ZN(n)=bNx(n)-aNy(n)ZN-1(n)=ZN(n-1)+bN-1x(n)-aN-1y(n)……Z1(n)=Z2(n-1)+b1x(n)-a1y(n)y(n)=Z1(n-1)+b0x(n)wi=[Z1(-1),Z2(-1),...,ZN(-1)]Ta=[1,-1,0.24];%定义差分方程左侧系数矢量ab=[1,-1];%定义差分方程右侧系数矢量bn=[0:20]‘;%生成仿真序列时间nx=n.^2;%生成激励信号x(n)wi=[-a(3)*(-2)-a(2)*(-1),-a(3)*(-1)];%定义初始状态wi[ywf]=filter(b,a,x,wi);%滤波得到输出序列常系数线性差分方程的求解常系数线性差分方程的求解a=[1,-0.9,0.3];b=0.05;n=[0:20]';u=(n=0);%生成激励序列u(n)uz=zeros(size(u));%生成零输入信号wi1=[-a(3);0];%定义两个初始状态wi2=[-a(2);-a(3)];[yzswf]=filter(b,a,[u,u]);%滤波得到零状态响应[yziwf]=filter(b,a,[uz,uz],[wi1,wi2]);%滤波得到零输入响应[ywf]=filter(b,a,[u,u],[wi1,wi2]);%滤波得到完全响应常系数线性差分方程的求解常系数线性差分方程的求解•对处理一维信号（矢量）的函数，如果以二维信号（矩阵）作输入，那么MATLAB会把它当成多个一维信号同时进行处理•和理论中常用的表达方法一样，MATLAB默认信号以列矢量形式表示•其他数学和逻辑函数也是如此单位样值(单位冲激)响应•可以用filter函数实现离散时间系统的单位样值响应•鉴于单位样值响应的重要性，MATLAB提供了用impz函数直接实现单位样值响应的方法a=[1,-0.5,0.6];b=[1,0,-0.3];n=[0:20]';[hi,t]=impz(b,a,n);%用impz函数计算单位样值响应hi(n)x=(n==0);%以单位样值序列为激励信号hf=filter(b,a,x);%经filter函数滤波得到hf(n)单位样值(单位冲激)响应•两种方法得到的结果完全相同卷积、解卷积•conv函数可以计算离散时间卷积和根据筛选特性，序列通过线性系统就是序列和系统单位样值响应进行卷积，因而卷积运算也可以用filter函数实现，即将两个待卷积序列分别理解为系统单位样值响应和激励序列a=0.8;%定义系数a=0.8N=6;%定义激励长度N=6n=[0:20]‘;%定义时间h=(a.^n).*(n=0);%生成单位样值响应h(n)x=(n=0&nN);%生成激励信号x(n)yc=conv(h,x);%用conv计算卷积，注意yc长度变化（方法1：直接卷积）yf1=filter(x,1,h);%用h(n)作为激励（方法2：把激励和单位样值响应调换）yf2=filter(1,[1,-a],x);%用x(n)作为激励（方法3：通过z变换得到h(n)的差分方程表示b=1,a=[1,–a]）卷积、解卷积卷积、解卷积•MATLAB提供[q,r]=deconv(b,a)函数实现解卷积其中b=conv(a,q)+r，即a和q卷积后再加上余量r得到b，而解卷积就要根据b和a解出q和rn=[0:20]‘;%定义抽样时间x=[1;0.5];%定义发射信号y=0.5.^n;%定义回波信号[h,r]=deconv(y,x);%解卷积求单位样值响应h0=0.5.^n.*(mod(n,2)==0);%定义理论结果做比较卷积、解卷积内容离散时间系统的时域分析z变换、z域分析离散时间系统的频响特性内容z变换、z域分析2利用z变换解差分方程z变换和逆z变换离散系统的系统函数z变换和逆z变换•用符号运算方法进行z变换，具体用ztrans函数实现symsn;%定义n为符号x1=(1/2)^n;%定义x1X1=ztrans(x1)%计算x1的z变换X1输出：X1=2*z/(2*z-1)x2=n*(n-1)/2;%定义x2X2=ztrans(x2)%计算x2的z变换X2X2s=simplify(X2)%将X2公式进行化简得到X2s输出：X2=1/2*z*(z+1)/(z-1)^3-1/2*z/(z-1)^2X2s=z/(z-1)^3z变换和逆z变换•同z变换相似，可以用符号运算函数iztrans实现逆z变换symsz;%符号法X=z^2/(z^2-1.5*z+0.5);x=iztrans(X)输出：x=-(1/2)^n+2a=[1,-1.5,0.5];%标准形式的分母系数b=1;[r,p,k]=residuez(b,a)%部分分式法求得r=[2,-1]’p=[1,0.5]’k=[]查z变换对应表得利用z变换解差分方程•零状态symsnabz;%定义符号n，a，b，zx=a^n*heaviside(n);%定义激励信号x(n)X=ztrans(x);%计算激励信号的变换X(z)H=1/(1-b*z^(-1));%由差分方程直接写出系统z变换式H(z)Y=H*X;%计算输出的变换式Y(z)=H(z)X(z)y1=iztrans(Y);%计算输出时域表达式y1(n)y=simplify(y1)%将y1(n)简化成y(n)•完全响应先取z变换，再由位移性质得到代入初始条件得到求逆变换即可得到响应y(n)序列的傅里叶变换（DTFT）•离散时间傅里叶变换——DTFT•离散傅里叶变换——DFT计算[-π,π]区间中的N个抽样值•快速傅里叶变换——FFT计算复杂度N2Nlog2N用FFT实现傅里叶变换的数值计算离散时间系统的频响特性•系统的差分方程为y(n)-a1y(n-1)-a2y(n-2)=b1x(n-1)•a1=1.1;a2=-0.7;b1=1;%为参数赋值•a=[1,-a1,-a2];%定义左侧系数•b=[0,b1];%定义右侧系数•figure;%生成新图框•subplot(2,1,1),zplane(b,a);%在第一个子图内绘制零、极点分布图•subplot(2,1,2),impz(b,a);%在第二个子图内绘制单位样值响应•figure,freqz(b,a);%生成新图框并绘制频率特性离散时间系统的频响特性离散时间系统的频响特性•系统差分方程为y(n)+b1y(n-1)+b2y(n-2)=b2x(n)+b1x(n-1)+x(n-2)•b1=-1.1;b2=0.6;%为参数赋值•a=[1,b1,b2];•b=[b2,b1,1];•figure;•subplot(2,1,1),zplane(b,a);%绘制零、极点分布图•subplot(2,1,2),impz(b,a);%绘制单位样值响应•figure,freqz(b,a);%绘制频率响应•n=[0:40]‘;%生成时间点•x1=sin(0.1*pi*n);%生成单频信号•x2=0*n;%准备方波信号•x2(mod(n,10)5)=1;%生成方波信号，周期是10•y1=filter(b,a,x1);%分别对两个信号滤波•y2=filter(b,a,x2);离散时间系统的频响特性傅里叶变换应用于通信系统[t,omg,FT,IFT]=prefourier([0,5],1000,[-250,250],1000);g=3*cos(10*t)+2*cos(20*t);%由定义生成基带信号f=g.*cos(100*t);%调制g0=f.*cos(100*t);%解调G0=FT*g0;%解调输出的频谱H=(omg30&omg-30);%定义低通滤波器H(jw)G1=G0.*H;%在频域进行低通滤波g1=IFT*G1;%逆变换得到时域输出G=FT*g;%计算其他信号的频谱F=FT*f;傅里叶变换应用于通信系统傅里叶变换应用于通信系统傅里叶变换应用于通信系统新建一个m文件，编写以下函数function[t,omg,FT,IFT]=prefourier(Trg,N,OMGrg,K)T=Trg(2)-Trg(1);t=linspace(Trg(1),Trg(2)-T/N,N)';OMG=OMGrg(2)-OMGrg(1);omg=linspace(OMGrg(1),OMGrg(2)-OMG/K,K)';FT=T/N*exp(-j*kron(omg,t.'));IFT=OMG/2/pi/K*exp(j*kron(t,omg.'));end输入参数：Trg:二维矢量，两个元素分别表示时域信号的起止时间；N:时域抽样数量；OMGrg:二维矢量，两个元素分别表示频谱的起止频率；K:频域抽样数量。输出参数：t:抽样时间；omg:抽样频率；FT:实现傅里叶变换的矩阵~U~及系数；IFT:实现傅里叶逆变换的矩阵~V~及系数。MATLAB实验：图像复原0,00,1...0,11,01,1...1,1,............1,01,1...1,1fffNfffNfxyfMfMfMN112(//)00112(//)001(,)(,)0,1,2,1;0,1,2,1(,)(,)0,1,2,1;0,1,2,1MNjuxMvyNxyMNjuxMvyNuvFuvfxyeMNuMvNfxyFuvexMyNMATLAB实验：图像复原(,)(,)*(,)(,)(,)(,)(,)(,)gxyhxyfxyxyGuvHuvFuvNuv空域：频域：f(x,y)是输入图像函数，表示该点的强度或灰度h(x,y)又称为点扩散函数(PSF)η(x,y)是噪声函数，这次作业我们将噪声设置为高斯噪声MATLAB实验：图像复原直接逆滤波(,)(,)ˆ(,)(,)(,)(,)GuvNuvFuvFuvHuvHuv维纳滤波（最小均方误差滤波）寻找一个使统计误差的函数最小维纳滤波ˆf22ˆ[()]eEff*222(,)(,)ˆ(,)(,)(,)|(,)|(,)|(,)|1(,)(,)|(,)|(,)/(,)fffHuvSuvFuvGuvSuvHuvSuvHuvGuvHuvHuvSuvSuv|H(u,v)|2=H*(u,v)H(u,v)Sη(u,v)=|N(u,v)|2，噪声功率Sf(u,v)=|F(u,v)|2，原图功率谱MATLAB实验：图像复原L-R迭代算法⊗和⊕分别为卷积运算和相关运算k为迭代次数1(,)ˆˆ(,)(,)(,)ˆ(,)(,)kkkgxyfxyfxyhxyfxyhxyMATLAB实验：图像复原评价标准，灰度平均梯度法GMG越大，表明图像越清晰222200ˆˆˆˆ[(1,)(,)][(,1)(,)]2(1)(1)MNijfijfijfijfijGMGMNMATLAB实验：图像复原可能用到的函数：fspecial，用于生成退化函数imnoise，用于加噪声imfilter，实现线性空间滤波fft2，ifft2，二维离散傅里叶变换