红眼消除-多媒体作业(附matlab程序)
摘要
红眼效应是指用闪光灯拍摄人物照片时,在人眼瞳孔中央形成的红眼现象。本文主要用matlab实现了一个红眼消除系统,该系统需要手动选定红眼区域，然后对红眼部分进行红眼修正和平滑操作。
绪论
红眼效应是指照相机在闪光灯模式下拍摄人像照片时，在照片中的人眼瞳孔处呈现红色斑点的现象。其成因是人的瞳孔在环境光比较暗时会放大，近距离闪光灯的强光经过放大的瞳孔，照在视网膜后的微血管组织上，反射回红色的光线，造成实际成像的照片呈现“红眼”状。“红眼”和一般人们所认知的眼睛颜色差别很大，降低了照片的质量，给摄影对象留下了遗憾。由于照相的机会往往无法重复获得，因此，消除照片中的红眼现象很有必要。
红眼主要受环境亮度（环境光较暗时红眼现象更明显）、对象年龄（小孩的红眼现象更明显）、闪光灯光线反射入镜头的角度（角度越小，红眼效应越强）和特定的人群或人种（在白种人中出现红眼现象的机会更多）的影响。人们根据红眼的成因，采用了一些方法来消除其对照片的影响，如提高环境亮度、使用外置闪光灯、缩短与被摄对象的距离、使被摄对象不直视照相机镜头等等。不少照相机都具有红眼减弱（Red-eye Reduction）功能，其工作方式就是在成像闪光之前加闪一次，使被摄对象的瞳孔在预闪后缩小，成像时红眼效应就会减轻一些。但该功能往往不能确保消除红眼现象，而且要求被摄对象在预闪时必须直视照相机，年龄较小的孩子注意力不易集中，常常达不到预期效果。
本文给出了一种消除红眼的思路。红眼图像取自数码相片，为RGB彩色，把它转换到HSI模型下，检测到红眼时，调低饱和度S，这样使得红眼效果减弱甚至被消除。在这里有必要介绍一下这两种颜色模型，以及它们之间的转换关系。
1 颜色模型
为了科学地定量描述和使用颜色，人们提出了各种颜色模型。目前常用的颜色模型按用途可分为三类：计算颜色模型、视觉颜色模型和工业颜色模型。
计算颜色模型用于进行有关颜色的理论研究。常见的RGB模型、CIE XYZ 模型、Lab模型等均属于此类型。
视觉颜色模型是指与人眼对颜色感知的视觉模型相似的模型，它主要用于色彩的理解，常见的有HSI模型、HSV模型和HSL模型。
工业颜色模型侧重于实际应用，包括彩色系统、彩色传输系统及电视传输系统等。如印刷中用的CMYK模型、电视系统中用的YUV模型、用于彩色图像压缩的YCbCr模型。
1.1RGB模型
RGB模型也称为加色法混色模型。它是以RGB三色光互相叠加来实现混色的方法，因而适合于显示器等发光体的显示。其混色规律是：以等量的红、绿、蓝基色光混合。
设颜色传感器把数字图像上的一个像素编码成（R,G,B），每个分量量化范围为[0,255]共256级。因此，RGB模型可以大约表示成1670万种颜色。这足以表示自然界的任一颜色，故又称其为24位真彩色。
RGB模型与显示器等设备有着很好的对应关系，RGB显示器中，有三种荧光粉能够分别发出红、绿、蓝三种颜色，三个相邻的荧光点构成一个像素，这些荧光点受到三束强度分别为c1 、c2、c3的电子束的轰击，会发生不同的亮度，通过物理上的叠加或混合，使可显示出相应的颜色。
1.2HSI模型
HSI模型是美国色彩学家蒙塞尔于1915年提出的，它反映了人的视觉系统感知彩色的方式，以色调、饱和度和强度三种基本特征量来感知颜色。色调H （Hue）与光波的波长有关，它表示人的感官对不同的颜色的感受，如红色、绿色、蓝色等，它可以表示一定的范围的颜色，如暖色、冷色等。饱和度S(Saturatuin)
表示颜色的纯度，纯光谱色是完全饱和的，加入白光会稀释饱和度。饱和度越大，颜色看起来就会越鲜艳，反之亦然。强度I(Intensity)，对应成像亮度和图像灰度，是颜色的明亮程度。
HSI模型的建立是基于两个重要的事实：a、I分量与图像的彩色信息无关；
b、H和S分量与人感受颜色的方式是紧密相联的。这些特点使得HIS模型非常适合彩色特性检测与分析。
HSI色彩空间和RGB色彩空间只是同一物理量的不同表示法,因而它们之间存在着如下转换关系。
(l)RGB转换到HSI
对任何3个[0,l]范围内的RGB值,其对应的HSI空间中的ISH分量
的计算公式如式（1-1）下所示：
式（1-1）
式中
式中计算出的H值的范围为[0°，180°],对应于G≥B.在G (2)HSI转换到RGB 假设s、I的值在[0,1]之间,R、G、B的值也在[0,l]之间,则HSI转换到RGB 的计算公式如式(1-2)、式(1-3)和式(1-4)所示(分成三个区间以利用对称性)。 当H在[0°,120°]之间: 式（1-2） 当H在[120°,240°]之间: 式（1-3） 当H在[240°,360°]之间: 式（1-4） 2.自动红眼消除 本实验的思路是先将获取自数码相机中的红眼图像，由RGB模型转换到HSI 模型下，再检测出红色通道的值，当它大于某个设定值后，让s通道的值为0，处理后的图像再转到RGB模型下。下面图2.1至图2.4是由这个方法试验出来的效果： 图2.1 小狗处理前图2.2小狗处理后 从上述实验的例子可以看出，本方法虽然简单，处理某些图像有点效果，但仅仅只限于红眼红色值远远大于其它部位的情况，所以处理之前需要对RGB图像的R通道进行检测。所以此方法还远达不到所谓的自动处理的要求。 3.手动选定区域红眼处理 本方法需要人工手动选定照片红眼区域，需要操作细心，可分多次选定区域。选定区域后，把RGB图像转换到HSI下进行处理，跟前面方法一样，对于红眼部分S通道为0，再把处理过的图像转换到RGB模型下。下面图3.1至图3.4是 该方法的实验效果。 图3.1 小狗处理前图3.2 小狗处理后 图3.3 人脸处理前图3.4 人脸处理后 上述实验可以看出，本方法效果明显，但操作麻烦，由于选取红眼区域全靠手工，而选取的方式是确定一些点，用这些点围成一个圈，来选取红眼区域。有时候选定一个红眼区域往往需要几遍才能成功。红眼区域一般是圆的，或者椭圆状的，这些手工选的点围成的圈是个多边形，所以想把红眼全包进去肯定会有一定的误差，带来一些不必要的处理。 4.总结 通过这次红眼处理，使我对图像处理有了新的认识，图像处理看似简单，实则不易。本次实验没用到比较难的算法，只是简单的实现了手动去红眼。本次设 计用matlab编程，这样我对matlab的相关功能又有了一次熟悉和训练的机会，不但把一些以前运用过的知识温习了一遍，更是在这次去红眼编程实现的训练中学习到了更多的函数，感受到了matlab的强大之处。 本次设计还有很多不足之处，现在已经出现了很多自动去红眼的解决方案。以后需要多加学习别人的先进经验，了解最新解决思路，争取以后我也能提出一种自动去红眼的方法。 附程序 function hsi=rgb2hsi(rgb) % RGB转换到HSI rgb=im2double(rgb); r=rgb(:,:,1); g=rgb(:,:,2); b=rgb(:,:,3); %Implement the conversion equations. num=0.5*((r-g)+(r-b)); den=sqrt((r-g).^2+(r-b).*(g-b)); theta=acos(num./(den+eps)); H=theta; H(b>g)=2*pi-H(b>g); H=H/(2*pi); num=min(min(r,g),b); den=r+g+b; den(den==0)=eps; S=1-3.*num./den; H(S==0)=0; I=(r+g+b)/3; %Combine all three results into an hsi image. hsi=cat(3,H,S,I); function rgb=hsi2rgb(hsi) % HSI转换到RGB H=hsi(:,:,1)*2*pi; S=hsi(:,:,2); I=hsi(:,:,3); %Implement the conversion equations. R=zeros(size(hsi,1),size(hsi,2)); G=zeros(size(hsi,1),size(hsi,2)); B=zeros(size(hsi,1),size(hsi,2)); % RG sector (0<=H<2*pi/3). idx=find((0<=H)&(H<2*pi/3)); B(idx)=I(idx).*(1-S(idx)); R(idx)=I(idx).*(1+S(idx).*cos(H(idx))./cos(pi/3-H(idx))); G(idx)=3*I(idx)-(R(idx)+B(idx)); %BG sector (2*pi/3<=H<4*pi/3). idx=find((2*pi/3<=H)&(H<4*pi/3)); R(idx)=I(idx).*(1-S(idx)); G(idx)=I(idx).*(1+S(idx).*cos(H(idx)-2*pi/3)./cos(pi-H(idx))); B(idx)=3*I(idx)-(R(idx)+G(idx)); %BR sector. idx=find((4*pi/3<=H)&(H<=2*pi)); G(idx)=I(idx).*(1-S(idx)); B(idx)=I(idx).*(1+S(idx).*cos(H(idx)-4*pi/3)./cos(5*pi/3-H(idx))); R(idx)=3*I(idx)-(G(idx)+B(idx)); rgb=cat(3,R,G,B); rgb=max(min(rgb,1),0); %本程序的功能是去除用户指定区域的红眼 clear all;clc;%内存清理 [file, path] = uigetfile('*.bmp', 'Image');%获取图像名及存储位置 if isequal(file,0) %未选择图像，退出 disp('取消选择'); else file=[path,file]; %获取图像路径 [filename,map]=imread(file); %读取图像 if ndims(filename)~=3 %是否为灰度图像,是退出并显示图像 disp('非彩色图像'); imshow(filename); else %去红眼操作 goagain=1; %操作控制变量，1进行操作 change=2; %参数更改控制变量，1更改预设参数 h1=1/6; %预设对rol区域内-pi/3 h2=1/6; s=0.2; fprintf('预设对rol区域内-pi/3 while goagain==1 fprintf('是否修改H、S预设值：1、是；2、否\n') change=input('请选择：'); %是否更改预设参数 while change~=1&change~=2 fprintf('错误命令！请重新输入~~\n'); change=input('请选择：'); end if change==1 %预设参数更改 h1=input('请输入H下限值（-2*pi/3~0）：')/pi/2; while (h1>0)||(h1<-1/3) fprintf('错误命令！请重新输入~~\n'); h1=input('请输入H下限值（-2*pi/3~0）：')/pi/2; end; h2=input('请输入H上限值（0~2*pi/3）：')/pi/2; while (h2<0)||(h2>1/3) fprintf('错误命令！请重新输入~~\n'); h2=input('请输入H上限值（0~2*pi/3）：')/pi/2 end; s=input('请输入S值(0~1)：'); while (s>1)||(s<0) fprintf('错误命令！请重新输入~~\n'); s=input('请输入S值(0~1)：'); end; change=2; end; B=roipoly(filename); %绘制要去红眼的rol区域，返回与filename大小相同的逻辑矩阵B，rol内像素值为1 %用户通过鼠标绘制任意多边形，点左键绘制顶点，点右键绘制最后的一 %个顶点 [hsi,H,S,I]=rgb2hsi(filename); %rgb->hsi rol=find(B~=0); %获取rol内像素数 [L,M]=size(rol); for i=1:L %遍历rol内像素，判定-1/6 if (((H(rol(i))>0)&(H(rol(i)) if S(rol(i))>s S(rol(i))=0; I(rol(i))=0; end; end; end; hsi=cat(3,H,S,I); %生成去红眼后的hsi图像 rgb=hsi2rgb(hsi); %hsi->rgb imshow(rgb); %显示去红眼后的rgb图像 fprintf('是否继续操作：1、是；2、否\n');%提示是否继续操作 goagain=input('请选择：'); while goagain~=1&goagain~=2 fprintf('错误命令！请重新输入~~\n'); goagain=input('请选择：'); end if goagain==1 filename=rgb; clf; end; end; end end; 几种常见窗函数及其MATLAB程序实现 2013-12-16 13:58 2296人阅读评论(0) 收藏举报 分类： Matlab（15） 数字信号处理中通常是取其有限的时间片段进行分析，而不是对无限长的信号进行测量和运算。具体做法是从信号中截取一个时间片段，然后对信号进行傅里叶变换、相关分析等数学处理。信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的。在FFT分析中为了减少或消除频谱能量泄漏及栅栏效应，可采用不同的截取函数对信号进行截短，截短函数称为窗函数，简称为窗。 泄漏与窗函数频谱的两侧旁瓣有关，对于窗函数的选用总的原则是，要从保持最大信息和消除旁瓣的综合效果出发来考虑问题，尽可能使窗函数频谱中的主瓣宽度应尽量窄，以获得较陡的过渡带；旁瓣衰减应尽量大，以提高阻带的衰减，但通常都不能同时满足这两个要求。 频谱中的如果两侧瓣的高度趋于零，而使能量相对集中在主瓣，就可以较为接近于真实的频谱。不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。信号的加窗处理，重要的问题是在于根据信号的性质和研究目的来选用窗函数。图1是几种常用的窗函数的时域和频域波形，其中矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低，如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用矩形窗，例如测量物体的自振频率等；布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。表1 是几种常用的窗函数的比较。 如果被测信号是随机或者未知的，或者是一般使用者对窗函数不大了解，要求也不是特别高时，可以选择汉宁窗，因为它的泄漏、波动都较小，并且选择性也较高。但在用于校准时选用平顶窗较好，因为它的通带波动非常小，幅度误差也较小。 《Matlab 编程训练》 作业 专 业 学生姓名 班级 学 号 指导教师 完成日期 实训一 MATLAB 语言介绍和数值计算 1．先求下列表达式的值，然后显示MATLAB 工作空间的使用情况并保存变量。 12 2sin851z e =+ . 2. 已知 1234413134787,2033657327A B --???? ????==???? ????-???? ，求下列表达式的值： (1) A+6*B 和A-B+I （其中I 为单位矩阵） A+6*B: A-B+I: (2)A*B和A.*B A*B程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] c=A*B 结果： A.*B程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] D=A.*B 结果： (3)A^3和A.^3 A^3程序： A=[12 34 -4;34 7 87;3 65 7] E=A^3 结果： A．^3程序： A=[12 34 -4;34 7 87;3 65 7] C=A.^3 (4)A/B及B\A A/B程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] C=A/B 结果： B\A程序： A=[12 34 -4;34 7 87;3 65 7] B=[1 3 -1;2 0 3;3 -2 7] D=B\A 结果： (5)将矩阵C=B\A的右下角2*2子矩阵赋给D, 并（3）保存变量（mat文件）程序： A=[12 34 -4;34 7 87;3 65 7]; B=[1 3 -1;2 0 3;3 -2 7]; C=B*inv(A); D=C(2:3,2:3) 结果： clear all U=形式背景 [m,n]=size(U); waiyan={}; waiyan{1}=1:m; temp1=[]; temp={}; for i=1:n temp1(1)=1; for j=2:m if(U(j,i)==1) temp1=[temp1,j];%得到形式背景 end end temp{i}=temp1; temp1=[]; end a=intersect(cell2mat(temp(1)),cell2mat(waiyan(1))); b{1}=a; waiyan=[waiyan,b]; b=[]; for i=2:n m=size(waiyan); for j=1:m(2) a=intersect(cell2mat(temp(i)),cell2mat(waiyan(j)));%求子集 if(~isziji(a,waiyan));%若原来不存在该概念 waiyan=[waiyan,a];%则加入该概念 end end end disp('概念格构造为：'); for i=1:max(size(waiyan)) t=num2str(cell2mat(waiyan(i))); disp(strcat('{',t,'}')); end ----------------------------------------------------------------------------- function logic=isziji(a,waiyan) m=size(waiyan); m1=max(size(a)); for i=1:max(m) m2=max(size(cell2mat(waiyan(i)))); 信息与随机性报告 随机走动 （1）随机走动回到零点的概率 a.一维随机走动：假设有一只青蛙，它处在一维坐标系的零点处，有1/2的概率向左跳，有1/2的概率往右跳。向左跳，坐标减1，向右跳，坐标加1。 进行10000次试验，青蛙走的最大步数为10000。 程序; clear all clc; b=0; for i=1:10000; a=0; for j=1:10000 x=rand; if x>0.5 a=a+1; else a=a-1; end if a==0; pp=j; b=b+1; break; end end end return1=b/10000;%返回的概率 运行结果： 返回的概率为99.12%，因此可以认为，一维随机走动一定会回到原点。 b.二维随机走动：假设青蛙处在二维坐标系中，每一次走动它向上向下向左向右移动的概率均为1/4，考虑它能回到原点的概率。 进行1000次试验，青蛙走的最大步数为1000000。 程序： clear all clc; total=0; for i=1:1000; a=0; b=0; for j=1:1000000 x=rand; y=rand; if x>0.5; x=1; else x=-1; end if y>0.5 a=a+x; else b=b+x; end if a==0 && b==0; pp=j; total=total+1; break; end end end return2=total/1000;%返回的概率 运行结果： 可以看到，青蛙回到原点的概率为97.63%，因此可以认为在二维随机走动中，青蛙一定是可以回到原点的。 c.三维随机走动：假设青蛙处在三维坐标系中，每一次走动它移动的方向有八个，每个方向的概率为1/8，考虑它能回到原点的概率。 进行1000次试验，青蛙走的最大步数为100000。 程序： clear all clc; total=0; for i=1:1000; 《现代设计方法》课程 关于黄金分割法和二次插值法的Matlab语言实现在《现代设计方法》的第二章优化设计方法中有关一维搜索的最优化方法的 一节里，我们学习了黄金非分割法和二次插值法。它们都是建立在搜索区间的优先确定基础上实现的。 为了便于方便执行和比较，我将两种方法都写进了一个程序之内，以选择的方式实现执行其中一个。下面以《现代设计方法》课后习题为例。见课本70页，第2—7题。原题如下： 求函数f(x)=3*x^2+6*x+4的最优点，已知单谷区间为[-3,4],一维搜索精度为0.4。 1、先建立函数f(x),f(x)=3*x^2+6*x+4。函数文件保存为：lee.m 源代码为：function y=lee(x) y=3*x^2+6*x+4; 2、程序主代码如下，该函数文件保存为：ll.m clear; a=input('请输入初始点'); b=input('请输入初始步长'); Y1=lee(a);Y2=lee(a+b); if Y1>Y2 %Y1>Y2的情况 k=2; Y3=lee(a+2*b); while Y2>=Y3 %直到满足“大，小，大”为止 k=k+1; Y3=lee(a+k*b); end A=a+b;B=a+k*b; elseif Y1 Matlab程序设计 班级 姓名 学号 《MATLAB程序设计》作业 1、考虑如下x-y 一组实验数据： x=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10] y=[1.2, 3, 4, 4, 5, 4.7, 5, 5.2, 6, 7.2] 分别绘出plot的原始数据、一次拟合曲线和三次拟合曲线，给出MATLAB代码和运行结果。 代码如下： x=[1,2,3,4,5,6,7,8,9,10]; y=[1.2,3,4,4,5,4.7,5,5.2,6,7.2]; plot(x,y); title('原始数据'); p=polyfit(x,y,1); q=polyval(p,x); figure,plot(x,q); title('一次拟合'); p=polyfit(x,y,2); q=polyval(p,x); figure,plot(x,q); title('二次拟合'); 运行结果如下： 1 2 3 4 5 6 7 8 9 10 12 3 4 5 6 7 8 原始数据 123 456789 102 2.5 3 3.54 4.5 55.56 6.57一次拟合 123456789 101 2 3 4 5 6 7 二次拟合 2、在[0，3π]区间，绘制y=sin(x)曲线（要求消去负半波，即(π，2π)区间内的函数值置零），求出曲线y 的平均值，以及y 的最大值及其最大值的位置。给出执行代码和运行结果。 代码如下： clear clc x=(0:0.01:3*pi); y=sin(x); plot(x,y); y1=(y>=0).*y; figure,plot(x,y1); 频率特性类题目 1 一个系统的开环传递函数为 ，试绘制其当K=5、30时系统的开环频率特性Nyquist 图，并判断系统的稳定性。 w=linspace(0.5,5,1000)*pi; sys1=zpk([],[0 -10 -2],100) sys2=zpk([],[0 -10 -2],600) figure(1) nyquist(sys1,w) title('system nyquist charts with k=5') figure(2) nyquist(sys2,w) title('system nyquist charts with k=30') 由图可知K=5时，开环Nyquist 曲线没有包围（-1，j0）点，所以系统稳定。 K=30时，开环Nyquist 曲线包围（-1，j0）点，所以系统不稳定。 2系统开环传递函数为 ，建立其零极点增益模型， 然后分别绘制当K=5、K=30时系统的开环频率特性Bode 图，并判断系统的稳定性。 sys1=zpk([],[0 -10 -2],100) sys2=zpk([],[0 -10 -2],600) figure(1) [Gm1,Pm1,Wcg1,Wcp1]=margin(sys1) bode(sys1) title('system bode charts with k=5'),grid figure(2) [Gm2,Pm2,Wcg2,Wcp2]=margin(sys2) bode(sys2) title('system bode charts with k=30'),grid 因为K=5时，Wcg>Wcp,所以系统稳定。 K=10时，Wcg MATLAB主要命令汇总 MATLAB函数参考 附录1.1 管理用命令 函数名功能描述函数名功能描述 addpath 增加一条搜索路径 rmpath 删除一条搜索路径 demo 运行Matlab演示程序 type 列出.M文件 doc 装入超文本文档 version 显示Matlab的版本号 help 启动联机帮助 what 列出当前目录下的有关文件 lasterr 显示最后一条信息 whatsnew 显示Matlab的新特性 lookfor 搜索关键词的帮助 which 造出函数与文件所在的目录 path 设置或查询Matlab路径 附录1.2管理变量与工作空间用命令 函数名功能描述函数名功能描述 clear 删除内存中的变量与函数 pack 整理工作空间内存 disp 显示矩阵与文本 save 将工作空间中的变量存盘 length 查询向量的维数 size 查询矩阵的维数 load 从文件中装入数据 who,whos 列出工作空间中的变量名 附录1.3文件与操作系统处理命令 函数名功能描述函数名功能描述 cd 改变当前工作目录 edit 编辑.M文件 delete 删除文件 matlabroot 获得Matlab的安装根目录 diary 将Matlab运行命令存盘 tempdir 获得系统的缓存目录 dir 列出当前目录的内容 tempname 获得一个缓存(temp)文件 ! 执行操作系统命令 附录1.4窗口控制命令 函数名功能描述函数名功能描述 echo 显示文件中的Matlab中的命令 more 控制命令窗口的输出页面format 设置输出格式 附录1.5启动与退出命令 函数名功能描述函数名功能描述 matlabrc 启动主程序 quit 退出Matlab环境 startup Matlab自启动程序 附录2 运算符号与特殊字符附录 2.1运算符号与特殊字符 函数名功能描述函数名功能描述 作业：用Matlab写C/A码生成器程序，并画生成码的方波图。 C/A码生成原理 C/A 码是用m 序列优选对组合形成的Gold 码。Gold码是由两个长度相同而互相关极大值为最小的m 序列逐位模2 相加所得到的码序列。它是由两个10 级反馈移位寄存器组合产生的，其产生原理如图1 所示。 图1 C/A码生成原理 发生器的抽头号为3和10，发生器的抽头号为2、3、6、8、9、10；发生器的第10位输出的数字即为码，而码是由的两个抽头的输出结果进行模2相加得到。 卫星的PRN码与延时的量是相关联的,对C/A码来说,每颗卫星都有特别的延时,如第1颗GPS卫星的G2 抽为2、6,第2颗为3、7,第3 颗为4、8,第4 颗为5、9 等，如图2所示。通过G2 相位选择可以产生结构不同的伪随机码,从而可以实现不同卫星之间的码分多址技术与卫星识别。 图2 prn序号与G2抽头、时延对应关系 基于MATLAB的GPS信号实现 编写成“codegen”程序，输入[ca_used]=codegen(svnum)，其中svnum为卫星号，ca_used 为得到的C/A码序列。程序具体实现流程如下： 在程序中定义一个数组，使得卫星号与G2的码片延时一一对应。 gs2=[5;6;7;8;17;18;139;140;141;251;252;254;255;256;257;258;469;470;471;472;473;474;509;512 ;513;514;515;516;859;860;861;862]; 定义两个1×1 023 的数组g1、g2 用来存放生成的Gold 码。定义一个全1 的10 位数组，作为移位寄存器，相当于G1、G2 生成模块的初值均置为全“1”。按原理式 1. p138 第6题在同一坐标轴中绘制下列两条曲线并标注两曲线交叉点。 >> t=0:0.01:pi; >> x1=t; >> y1=2*x1-0.5; >> x2=sin(3*t).*cos(t); >> y2=sin(3*t).*sin(t); >> plot(x1,y1,'r-',x2,y2,'g-') >> axis([-1,2,-1.5,1]) >> hold on >> s=solve('y=2*x-0.5','x=sin(3*t)*cos(t)','y=sin(3*t)*sin(t)'); >> plot(double(s.x),double(s.y),'*'); 截图： p366 第4题绘制极坐标曲线，并分析对曲线形状的影响。 function [ output_args ] = Untitled2( input_args ) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here theta=0:0.01:2*pi; a=input('请输入a的值：'); b=input('请输入b的值：'); n=input('请输入n的值：'); rho=a*sin(b+n*theta); polar(theta,rho,'k'); end 下面以a=1，b=1，n=1的极坐标图形为基础来分析a、b、n的影响。 对a的值进行改变：对比发现a只影响半径值的整倍变化 对b的值进行改变：对比发现b的值使这个圆转换了一定的角度 对n的值进行改变：对比发现当n>=2时有如下规律 1、当n为整数时，图形变为2n个花瓣状的图形 2、当n为奇数时，图形变为n个花瓣状的图形 分别让n为2、3、4、5 MATLAB简单程序大全 求特征值特征向量 A=[2 3 4;1 5 9;8 5 2] det(A) A' rank(A) inv(A) rref(A) eig(A)%求特征值和特征向量 卫星运行问题 h=200,H=51000,R=6378; a=(h+H+2*R)/2; c=(H-h)/2; b=(a^2-c^2)^(1/2); e=c/a; f=sqrt(1-exp(2).*cos(t)^2); l=int(f,t,0,pi/2) L=4*a.*l 动态玫瑰线 n=3;N=10000; theta=2*pi*(0:N)/N; r=cos(n*theta); x=r.*cos(theta); y=r.*sin(theta); comet(x,y) 二重积分 syms x y f=x^2*sin(y); int(int(f,x,0,1),y,0,pi) ezmesh(f,[0,1,0,pi]) 函数画图 syms x;f=exp(-0.2*x)*sin(0.5*x); ezplot(f,[0,8*pi]) 玫瑰线 theta=0:0.01:2*pi; r=cos(3*theta); polar(theta,r,'r') 求x^2+y^2=1和x^2+z^2=1所围成的体积 syms x y z R r=1; Z=sqrt(1-x^2); y0=Z; V=8*int(int(Z,y,0,y0),x,0,1) 求导数及图像 f='1/(5+4*cos(x))'; subplot(1,2,1);ezplot(f) f1=diff(f) subplot(1,2,2);ezplot(f1) 绕x轴旋转 t=(0:20)*pi/10; r=exp(-.2*t).*sin(.5*t); theta=t; x=t'*ones(size(t)); y=r'*cos(theta); z=r'*sin(theta); mesh(x,y,z) colormap([0 0 0]) 某年是否闰年 year=input('input year:='); n1=year/4; n2=year/100; n3=year/400; if n1==fix(n1)&n2~=fix(n2) disp('是闰年') elseif n1==fix(n1)&n3==fix(n3) disp('是闰年') else 基于m a t l a b程序实现 人脸识别 TYYGROUP system office room 【TYYUA16H-TYY-TYYYUA8Q8- 基于m a t l a b程序实现人脸识别 1.人脸识别流程 基于YCbCr颜色空间的肤色模型进行肤色分割。在YCbCr色彩空间内对肤色进行了建模发现，肤色聚类区域在Cb—Cr子平面上的投影将缩减，与中心区域显着不同。采用这种方法的图像分割已经能够较为精确的将人脸和非人脸分割开来。 人脸识别流程图 2.人脸识别程序 （1）人脸和非人脸区域分割程序 function result = skin(Y,Cb,Cr) %SKIN Summary of this function goes here % Detailed explanation goes here a=; b=; ecx=; ecy=; sita=; cx=; cy=; xishu=[cos(sita) sin(sita);-sin(sita) cos(sita)]; %如果亮度大于230，则将长短轴同时扩大为原来的倍 if(Y>230) a=*a; b=*b; end %根据公式进行计算 Cb=double(Cb); Cr=double(Cr); t=[(Cb-cx);(Cr-cy)]; temp=xishu*t; value=(temp(1)-ecx)^2/a^2+(temp(2)-ecy)^2/b^2; %大于1则不是肤色，返回0；否则为肤色，返回1 if value>1 result=0; else result=1; end end （2）人脸的确认程序 function eye = findeye(bImage,x,y,w,h) %FINDEYE Summary of this function goes here % Detailed explanation goes here part=zeros(h,w); %二值化 for i=y:(y+h) for j=x:(x+w) if bImage(i,j)==0 part(i-y+1,j-x+1)=255; else part(i-y+1,j-x+1)=0; end end end [L,num]=bwlabel(part,8); %如果区域中有两个以上的矩形则认为有眼睛 if num<2 eye=0; Matlab程序代码： clc; clear; N=20; T=0.1 t=0:T:N m=length(t) syms x1 x2 x3 fx=[0;x1+x2^2;x1-x2] gx=[exp(x2);exp(x2);0] hx=x3; R=10*eye(1) Q=[10 0 0;0 1 0;0 0 1] A=[0 1 0;0 0 1;0 0 0] B=[0;0;1] SS=B*inv(R)*B' [p1,p2,lamp,perr,wellposed,P]=aresolv(A,Q,SS) z1=hx z2=[diff(hx,x1) diff(hx,x2) diff(hx,x3)]*fx z3=[diff(z2,x1), diff(z2,x2), diff(z2,x3)]*fx ax=[diff(z3,x1), diff(z3,x2), diff(z3,x3)]*fx bx=[diff(z3,x1), diff(z3,x2), diff(z3,x3)]*gx z=[z1;z2;z3] k=inv(R)*B'*P %diff(z)=A*z+B*v=(A-B*K)*Z %x(0)=[1;0;0] abk=A-B*k x1(1)=1 x2(1)=0 x3(1)=1 z1(1)=x3(1) z2(1)=x1(1)-x2(1) z3(1)=-(x1+x2^2) for i=2:m z1(i)=z1(i-1)+T*(abk(1,1)*z1(i-1)+abk(1,2)*z2(i-1)+abk(1,3)*z3(i-1)) z2(i)=z2(i-1)+T*(abk(2,1)*z1(i-1)+abk(2,2)*z2(i-1)+abk(2,3)*z3(i-1)) z3(i)=z3(i-1)+T*(abk(3,1)*z1(i-1)+abk(3,2)*z2(i-1)+abk(3,3)*z3(i-1)) 遗传算法的原理及MATLAB程序实现 1 遗传算法的原理 1.1 遗传算法的基本思想 遗传算法(genetic algorithms，GA)是一种基于自然选择和基因遗传学原理，借鉴了生物进化优胜劣汰的自然选择机理和生物界繁衍进化的基因重组、突变的遗传机制的全局自适应概率搜索算法。 遗传算法是从一组随机产生的初始解(种群)开始，这个种群由经过基因编码的一定数量的个体组成，每个个体实际上是染色体带有特征的实体。染色体作为遗传物质的主要载体，其内部表现(即基因型)是某种基因组合，它决定了个体的外部表现。因此，从一开始就需要实现从表现型到基因型的映射，即编码工作。初始种群产生后，按照优胜劣汰的原理，逐代演化产生出越来越好的近似解。在每一代，根据问题域中个体的适应度大小选择个体，并借助于自然遗传学的遗传算子进行组合交叉和变异，产生出代表新的解集的种群。这个过程将导致种群像自然进化一样，后代种群比前代更加适应环境，末代种群中的最优个体经过解码，可以作为问题近似最优解。 计算开始时，将实际问题的变量进行编码形成染色体，随机产生一定数目的个体，即种群，并计算每个个体的适应度值，然后通过终止条件判断该初始解是否是最优解，若是则停止计算输出结果，若不是则通过遗传算子操作产生新的一代种群，回到计算群体中每个个体的适应度值的部分，然后转到终止条件判断。这一过程循环执行，直到满足优化准则，最终产生问题的最优解。图1-1给出了遗传算法的基本过程。 1.2 遗传算法的特点 1.2.1 遗传算法的优点 遗传算法具有十分强的鲁棒性，比起传统优化方法，遗传算法有如下优点: 1. 遗传算法以控制变量的编码作为运算对象。传统的优化算法往往直接利用控制变量的实际值的本身来进行优化运算，但遗传算法不是直接以控制变量的值，而是以控制变量的特定形式的编码为运算对象。这种对控制变量的编码处理方式，可以模仿自然界中生物的遗传和进化等机理，也使得我们可以方便地处理各种变量和应用遗传操作算子。 2. 遗传算法具有内在的本质并行性。它的并行性表现在两个方面，一是遗传 开始 初始化，输入原始参 数及给定参数，gen=1 染色体编码，产生初始群体 计算种群中每个个体的适应值 终止条件的判断, N gen=gen+1 选择 交叉 Y 变异 新种群 输出结果 结束 图1-1 简单遗传算法的基本过程 实验二 MATLAB 程序设计 一、 实验目的 1.掌握利用if 语句实现选择结构的方法。 2.掌握利用switch 语句实现多分支选择结构的方法。 3.掌握利用for 语句实现循环结构的方法。 4.掌握利用while 语句实现循环结构的方法。 5.掌握MATLAB 函数的编写及调试方法。 二、 实验的设备及条件 计算机一台（带有MATLAB7.0以上的软件环境）。 M 文件的编写： 启动MATLAB 后，点击File|New|M-File ，启动MATLAB 的程序编辑及调试器（Editor/Debugger ），编辑以下程序，点击File|Save 保存程序，注意文件名最好用英文字符。点击Debug|Run 运行程序，在命令窗口查看运行结果，程序如有错误则改正 三、 实验容 1.编写求解方程02=++c bx ax 的根的函数（这个方程不一定为一元二次方程，因c b a 、、的不同取值而定），这里应根据c b a 、、的不同取值分别处理，有输入参数提示，当0~,0,0===c b a 时应提示“为恒不等式!”。并输入几组典型值加以检验。 （提示：提示输入使用input 函数） 2.输入一个百分制成绩，要求输出成绩等级A+、A 、B 、C 、D 、E 。其中100分为A+，90分～99分为A ，80分～89分为B ，70分～79分为C ，60分～69分为D ，60分以下为E 。 要求：（1）用switch 语句实现。 （2）输入百分制成绩后要判断该成绩的合理性，对不合理的成绩应输出出错信息。 （提示：注意单元矩阵的用法） 3.数论中一个有趣的题目：任意一个正整数，若为偶数，则用2除之，若为奇数，则与3相乘再加上1。重复此过程，最终得到的结果为1。如： 2→1 3→10→5→16→8→4→2→1 6→3→10→5→16→8→4→2→1 1．全景图到穹景图 这个程序我最初是用FreeImage写的，这两天改成了matlab，再不贴上来，我就要忘了。 看到一篇文章有这样的变换，挺有意思的，就拿来试了一下，文章点此。 全景图到穹顶图变换，通俗的说就是将全景图首尾相接做成一个圆环的样子。 先看下面这图： 下面的矩形就是我们要处理的全景图，上面的矩形是变换后的图像。下面图像的底边对应穹顶图的圆，顶边对应穹顶图的外圆，当然，反过来也是可以的。 程序流程： 1.定义穹顶图圆和外圆的半径，变换后的像素就填充在这个外半径的圆环中。 2.遍历穹顶图，当所处理当前像素位于圆环，则通过极坐标反变换去全景图中寻找相应位置的像素进行填充。 3.遍历完图像就行了。 用的技巧和图像旋转或放大缩小都是类似的。 处理结果： 原图： 结果： matlab代码如下： clear all; close all; clc; img=imread('pan.jpg'); imshow(img); [m,n]=size(img); r1=100; %环半径 r2=r1+m; %外环半径 imgn=zeros(2*r2,2*r2); [re_m,re_n]=size(imgn); for y=1:re_m for x=1:re_n dis_x=x-re_n/2; dis_y=y-re_m/2; l=sqrt(dis_x^2+dis_y^2); if l<=r2 && l>=r1 theta=0; if y>re_m/2 theta=atan2(dis_y,dis_x); end if y 动态规划方法的Matlab 实现与应用 动态规划(Dynamic Programming)是求解决策过程最优化的有效数学方法，它是根据“最优决策的任何截断仍是最优的”这最优性原理，通过将多阶段决策过程转化为一系列单段决策问题，然后从最后一段状态开始逆向递推到初始状态为止的一套最优化求解方法。 1．动态规划基本组成 (1) 阶段 整个问题的解决可分为若干个阶段依次进行，描述阶段的变量称为阶段变量，记为k (2) 状态 状态表示每个阶段开始所处的自然状况或客观条件，它描述了研究问题过程的状况。各阶段状态通常用状态变量描述，用k x 表示第k 阶段状态变量，n 个阶段决策过程有n+ 1个状态。 (3) 决策 从一确定的状态作出各种选择从而演变到下一阶段某一状态，这种选择手段称为决策。描述决策的变量称为决策变量，决策变量限制的取值范围称为允许决策集合。用()k k u x 表示第k 阶段处于状态k x 时的决策变量，它是k x 的函数。用()k k D x Dk(xk)表示k x 的允许决策的集合。 (4) 策略 每个阶段的决策按顺序组成的集合称为策略。由第k 阶段的状态k x 开始到终止状态的后部子过程的策略记为{}11(),(),,()k k k k n n u x u x u x ++ 。可供选择的策略的范围称为允许策略集合，允许策略集合中达到最优效果的策略称为最优策略。从初始状态* 11()x x =出发，过程按照最优策略和状态转移方程演变所经历的状态序列{ } **** 121,,,,n n x x x x + 称为最优轨线。 (5) 状态转移方程 如果第k 个阶段状态变量为k x ，作出的决策为k u ，那么第k+ 1阶段的状态变量1k x +也被完全确定。用状态转移方程表示这种演变规律，记为1(,)k k k x T x u +=。 (6) 指标函数 指标函数是系统执行某一策略所产生结果的数量表示，是衡量策略优劣的数量指标，它定义在全过程和所有后部子过程上，用()k k f x 表示。过程在某阶段j 的阶段指标函数是衡量该阶段决策优劣数量指标，取决于状态j x 和决策j u ，用(,)j j j v x u 表示。 2．动态规划基本方程 (){} 11()min ,,(),()k k k k k k k k k k f x g v x u f x u D x ++=∈???? Matlab 实现 (dynprog.m 文件) function [p_opt,fval]=dynprog (x,DecisFun,SubObjFun,TransFun,ObjFun) % x 是状态变量，一列代表一个阶段的所有状态； % M-函数DecisFun(k,x) 由阶段k 的状态变量x 求出相应的允许决策变量; % M-函数SubObjFun(k,x,u) 是阶段指标函数, % M-函数ObjFun(v,f) 是第k 阶段至最后阶段的总指标函数 % M-函数TransFun(k,x,u) 是状态转移函数, 其中x 是阶段k 的某状态变量, u 是相应的决策变量; %输出 p_opt 由4列构成，p_opt=[序号组;最优策略组;最优轨线组;指标函数值组]； %输出 fval 是一个列向量，各元素分别表示p_opt 各最优策略组对应始端状态x 的最优函数值。 Matlab 上机实验一、二 3.求下列联立方程的解???????=+-+-=-+=++-=--+4 1025695842475412743w z y x w z x w z y x w z y x >> a=[3 4 -7 -12;5 -7 4 2;1 0 8 -5;-6 5 -2 10]; >> b=[4;4;9;4]; >> c=a\b 4．设???? ??????------=81272956313841A ，??????????-----=793183262345B ，求C1=A*B’;C2=A’*B;C3=A.*B,并求上述所有方阵的逆阵。 >> A=[1 4 8 13;-3 6 -5 -9;2 -7 -12 -8]; >> B=[5 4 3 -2;6 -2 3 -8;-1 3 -9 7]; >> C1=A*B' >> C2=A'*B >> C3=A.*B >> inv(C1) >> inv(C2) >> inv(C3) 5．设 ?? ????++=)1(sin 35.0cos 2x x x y ，把x=0~2π间分为101点，画出以x 为横坐标，y 为纵坐标的曲线。 >> x=linspace(0,2*pi,101); >> y=cos(x)*(0.5+(1+x.^2)\3*sin(x)); >> plot(x,y,'r') 6．产生8×6阶的正态分布随机数矩阵R1, 求其各列的平均值和均方差。并求该矩阵全体数的平均值和均方差。 （mean var ） a=randn(8,6) mean(a) var(a) k=mean(a) k1=mean(k) i=ones(8,6) i1=i*k1 i2=a-i1 i3=i2.*i2 g=mean(i3) g2=mean(g) 正弦波的源程序： (一)，用到的函数 1,f2t函数 function x=f2t(X) global dt df t f T N %x=f2t(X) %x为时域的取样值矢量 %X为x的傅氏变换 %X与x长度相同并为2的整幂 %本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)]; x=ifft(X)/dt; end 2,t2f函数。 function X=t2f(x) global dt df N t f T %X=t2f(x) %x为时域的取样值矢量 %X为x的傅氏变换 %X与x长度相同,并为2的整幂。 %本函数需要一个全局变量dt(时域取样间隔) H=fft(x); X=[H(N/2+1:N),H(1:N/2)]*dt; end (二),主程序。 1，%（1）绘出正弦信号波形及频谱 global dt df t f N close all k=input('取样点数＝2^k, k取10左右'); if isempty(k), k=10; end f0=input('f0=取1(kz)左右'); if isempty(f0), f0=1; end N=2^k; dt=0.01; %ms df=1/(N*dt); %KHz T=N*dt; %截短时间 Bs=N*df/2; %系统带宽 f=[-Bs+df/2:df:Bs]; %频域横坐标 t=[-T/2+dt/2:dt:T/2]; %时域横坐标 s=sin(2*pi*f0*t); %输入的正弦信号 S=t2f(s); %S是s的傅氏变换 a=f2t(S); %a是S的傅氏反变换 a=real(a); as=abs(S); subplot(2,1,1) %输出的频谱 plot(f,as,'b'); grid axis([-2*f0,+2*f0,min(as),max(as)]) xlabel('f (KHz)') ylabel('|S(f)| (V/KHz)') %figure(2) subplot(2,1,2) plot(t,a,'black') %输出信号波形画图grid axis([-2/f0,+2/f0,-1.5,1.5]) xlabel('t(ms)') ylabel('a(t)(V)') gtext('频谱图') 最佳基带系统的源程序： (一)，用到的函数 f2t函数和t2f函数。代码>> (二),主程序 globaldt t f df N T close all clear Eb_N0 Pe k=input('取样点数＝2^k, k取13左右'); if isempty(k), k=13; end z=input('每个信号取样点数＝2^z, z 关闭文件用fclose函数，调用格式为：sta＝fclose(fid)说明：该函数关闭fid所表示的文件。其调用格式为：[A，COUNT]=fscanf(fid，format，size)说明：其中A用来存放读取的数据，COUNT返回所读取的数据元素个数，fid为文件句柄，format用来控制读取的数据格式，由%加上格式符组成，常见的格式符有：d（整型）、f（浮点型）、s（字符串型）、c（字符型）等，在%与格式符之间还可以插入附加格式说明符，如数据宽度说明等。 matlab fprintf.数据的格式化输出：fprintf(fid, format, variables)fprintf（fid，format，A）说明：fid为文件句柄，指定要写入数据的文件，format是用来控制所写数据格式的格式符，与fscanf函数相同，A是用来存放数据的矩阵。>> fid=fopen(""d:\char1.txt"",""w"");>> fid1=fopen(""d:\char1.txt"",""rt"");matlab读txt文件fid=fopen(""fx.txt"",""r"");%得到文件号[f,count]=fscanf(fid,""%f %f"",[12,90]);%把文件号1的数据读到f中。 matlab函数fgetl和fgets：按行读取格式文本函数Matlab提供了两个函数fgetl和fgets来从格式文本文件读取行，并存储到字符串向量中。这两个函数集几乎相同；不同之处是，fgets拷贝新行字符到字符向量，而fgetl则不。下面的M-file函数说明了fgetl的一个可能用法。此函数使用fgetl一次读取一整行。while f eof(fid) == 0 tline = fgetl(fid); %用Fourier变换求取信号的功率谱---周期图法 clf; Fs=1000; N=256;Nfft=256;%数据的长度和FFT所用的数据长度 n=0:N-1;t=n/Fs;%采用的时间序列 xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N); Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 subplot(2,1,1),plot(f,Pxx);%绘制功率谱曲线 xlabel('频率/Hz');ylabel('功率谱/dB'); title('周期图N=256');grid on; Fs=1000; N=1024;Nfft=1024;%数据的长度和FFT所用的数据长度 n=0:N-1;t=n/Fs;%采用的时间序列 xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N); Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列 subplot(2,1,2),plot(f,Pxx);%绘制功率谱曲线 xlabel('频率/Hz');ylabel('功率谱/dB'); title('周期图N=256');grid on; %用Fourier变换求取信号的功率谱---分段周期图法 %思想：把信号分为重叠或不重叠的小段，对每小段信号序列进行功率谱估计，然后取平均值作为整个序列的功率谱 clf;几种常见窗函数及其MATLAB程序实现
MATLAB编程作业
概念格matlab程序
随机走动-附matlab程序仿真
matlab编程实现求解最优解
MATLAB程序设计作业
matlab程序大全答案
matlab代码大全
CA码生成原理及matlab程序实现
matlab程序设计第三章课后习题答案
MATLAB简单程序大全
基于matlab程序实现人脸识别
Matlab程序代码
遗传算法的原理及MATLAB程序实现
实验二--MATLAB程序的设计(含实验报告)
MATLAB程序大全
基于Matlab的动态规划程序实现
Matlab编程与应用习题和一些参考答案
Matlab源程序代码
如何编写MATLAB程序才能实现对