个人收集的matlab代码截止到20110528
matlab代码
1-3%fft2Dbitmap_image.m%Simulation of Fourier transformation of bitmap imagesclearI=imread('triangle.bmp','bmp'); %Input bitmap imageI=I(:,:,1);figure(1) %displaying inputcolormap(gray(255));image(I)axis offFI=fft2(I);FI=fftshift(FI);max1=max(FI);max2=max(max1);scale=1.0/max2;FI=FI.*scale;figure(2) %Gray scale image of the absolute value of transform colormap(gray(255));image(10*(abs(256*FI)));axis off1-4%correlation.mclearI1=imread('smiley1.bmp','bmp'); %Input image 1 (reference image)I1=I1(:,:,1);figure(1) %displaying input image 1colormap(gray(255));image(I1)axis offFI1=fft2(I1);max1=max(FI1);max2=max(max1);scale=1.0/max2;FI1=FI1.*scale;I2=imread('smiley1.bmp','bmp'); %Input image 2 (image to be recognized) I2=I2(:,:,1);figure(2) %displaying input image 2colormap(gray(255));image(I2)axis offFI2=fft2(I2);max1=max(FI2);max2=max(max1);scale=1.0/max2;FI2=FI2.*scale;FPR=FI1.*conj(FI2);%calculating correlationPR=ifft2(FPR);PR=fftshift(PR);max1=max(PR);max2=max(max1);scale=1.0/max2;PR=PR.*scale;figure(3)%display of correlation in spatial domaincolormap(gray(255));image(abs(256*PR));axis offfftshift只是将fft2的结果移了下位，fft2的左下部分和右上部分对调，左上部分和右下部分对调，结果fft2的零频移到fft2得到矩阵的中心，这时可以看到中心一个亮点，要不然零频就分布在矩阵的四个角，中间一片黑2-1%Fresnel_diffraction.m%Simulation of Fresnel diffraction of a square aperture%Adapted from "Contemporary optical image processing with MATLAB®," %by T.-C. Poon and P. P. Banerjee, Elsevier 2001, pp. 64-65. clearL=1; %L : length of display areaN=256; %N : number of sampling pointsdx=L/(N-1); % dx : step size%Create square image, M by M square, rect(x/a), M=odd numberM=111;a=M/256R=zeros(256); %assign a matrix (256x256) of zerosr=ones(M); %assign a matrix (MxM) of onesn=(M-1)/2;R(128-n:128+n,128-n:128+n)=r;%End of creating input image%Axis Scalingfor k=1:256X(k)=1/255*(k-1)-L/2;Y(k)=1/255*(k-1)-L/2;%Kx=(2*pi*k)/((N-1)*dx)%in our case, N=256, dx=1/255Kx(k)=(2*pi*(k-1))/((N-1)*dx)-((2*pi*(256-1))/((N-1)*dx))/2;Ky(k)=(2*pi*(k-1))/((N-1)*dx)-((2*pi*(256-1))/((N-1)*dx))/2;end%Fourier transformation of RFR=(1/256)^2*fft2(R);FR=fftshift(FR);%Free space impulse response function% The constant factor exp(-jk0*z) is not calculated%sigma=ko/(2*z)=pi/(wavelength*z)%z=5cm,red light=0.6328*10^-4(cm)sigma=pi/((0.6328*10^-4)*5);for r=1:256,for c=1:256,%compute free-space impulse response with Gaussian apodization against aliasingh(r,c)=j*(sigma/pi)*exp(-4*200*(X(r).^2+Y(c).^2))*exp(-j*sigma*(X(r). ^2+Y(c).^2));endendH=(1/256)^2*fft2(h);H=fftshift(H);HR=FR.*H;H=(1/256)^2*fft2(h);H=fftshift(H);HR=FR.*H;hr=ifft2(HR);hr=(256^2)*hr;hr=fftshift(hr);%Image of the rectangle objectfigure(1)image(X,Y,255*R);colormap(gray(256));axis squarexlabel('cm')ylabel('cm')% plot of cross section of squarefigure(2)plot(X+dx/2,R(:,127))gridaxis([-0.5 0.5 -0.1 1.2])xlabel('cm')figure(3)plot(X+dx/2,abs(hr(:,127)))gridaxis([-0.5 0.5 0 max(max(abs(hr)))*1.1]) xlabel('cm')a =0.43362-2%Fresnel_zone_plate.m%Adapted from "Contemporary optical image processing with MATLAB®,"%by T.-C. Poon and P. P. Banerjee, Elsevier 2001, pp.177-178.%%display function is 1+sin(sigma*((x-x0)^2+(y-y0)^2)). All scales are arbitrary.%sigma=pi/(wavelength*z)clear;z0=input ('z0, distance from the point object to film, enter z0 (from 2 to 10)=');x0=input ('Inputting the location of the center of the FZP x0=y0,enter x0 (from -8 to 8) =');ROWS=256;COLS=256;colormap(gray(255))sigma=1/z0;y0=-x0;y=-12.8;for r=1:COLS,x=-12.8;for c=1:ROWS, %compute Fresnel zone platefFZP(r,c)=exp(j*sigma*(x-x0)*(x-x0)+j*sigma*(y-y0)*(y-y0));x=x+.1;endy=y+.1;end%normalizationmax1=max(fFZP);max2=max(max1);scale=1.0/max2;fFZP=fFZP.*scale;R=127*(1+imag(fFZP)); figure(1)image(R);axis square onaxis off。
matlab函数代码
MATLAB函数代码详解1. 任务名称：MATLAB函数代码MATLAB是一种功能强大的科学计算软件，它使用自己的语言和编程环境。
在MATLAB中，函数是用来封装可重用代码的基本工具之一。
本文将深入探讨MATLAB 函数代码并详细介绍其使用方法。
2. 简介MATLAB函数是一种可重复使用和模块化的代码块，它可以接受输入参数并返回一个或多个输出结果。
函数是用来实现特定功能的，可以在一个程序中多次调用，而不需要重复编写代码。
它提供了隔离和封装代码的能力，使得程序更加可维护和易于调试。
3. MATLAB函数的定义和结构MATLAB函数的定义以关键字”function”开头，后面跟着函数名和输入输出参数列表。
下面是一个简单的MATLAB函数的例子：function output = myFunction(input)% 函数体output = someOperation(input);end•function关键字用于声明函数的开始。
•函数名用于唯一标识一个函数。
•输入参数列表包含函数需要使用的输入数据，每个参数由逗号分隔。
•输出参数列表用于指定函数返回的结果，可以有一个或多个输出。
函数体是函数中的主要部分，它包含了实现特定功能的MATLAB代码。
函数体中可以包含各种语句、循环、条件判断等。
在函数体中，可以使用输入参数和本地变量进行计算，并通过输出参数返回结果。
函数体具有局部作用域，即只能在函数内部使用声明的变量。
此外，函数还可以具有嵌套函数，这些嵌套函数可以在主函数内部定义。
4. MATLAB函数的调用方法调用MATLAB函数需要使用函数名和输入参数。
下面是一个函数调用的示例：output = myFunction(input);实际上，函数调用是将输入参数传递给函数，并且函数执行相应的操作并返回结果。
调用一个函数可以在程序的任何地方进行，只要函数在调用之前已经定义。
通过函数调用，可以将函数拆分为更小的可管理的部分。
matlab程序代码及学习小结
1.信号傅里叶分解代码：（1）离散傅里叶变换N=256;dt=0.02;%数据的个数和采样间隔n=0:N-1; t=n*dt;%序号序列和时间序列x=sin(2*pi*t);%合成信号m=floor(N/2)+1;%分解a，b的最大序号值为分解的N/2个参数加参数a0a=zeros(1,m);b=zeros(1,m);%产生a，b两个为零的序列for k=0:m-1,for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N);b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);%MATLAB中的数组序号只能从1开始endc(k+1)=sqrt(a(k+1)^2+b(k+1)^2);endsubplot(211);plot(t,x);title('原始信号');xlabel('时间/s'); f=(0:m-1)/(N*dt);subplot(212);plot(f,c);title('Fourier变换’)；xlabel('频率/Hz');ylabel('振幅') （2）快速傅里叶变换fs=input('please input the fs:');%设定采样频率N=input('please input the N:');%设定数据长度(t,x)(先导入t,x的值)subplot(211);plot(t,x);%作正弦信号的时域波形axis([0,1,-1,1]);%设置坐标轴长度title('正弦信号时域波形');%进行FFT变换并做频谱图y=fft(x,N);%进行fft变换mag=abs(y);%求幅值f=(0:N-1)*fs/N;%横坐标频率的表达式为f=(0:M-1)*Fs/M; subplot(212);plot(f,mag);%做频谱图title('正弦信号幅频谱图');2.数据导入记事本数据导入代码方法一：x = load('E:\文件名.txt')plot(x(:, 1), x(:, 2))方法二：x=importdata(文件名)x1=x.data(行,列)%例如x1=x.data(1:100,2)(记事本中只有数据时，代码)x=importdata(文件名)x1=x(行,列)3.转速数据处理练习（1）转速处理代码离散变换a=importdata('speed1.txt')t=a(1:1000,1)x=a(1:1000,2)N=1000;dt=0.001;n=0:N-1; t=n*dt;%序号序列和时间序列m=floor(N/2)+1;%分解a，b的最大序号值为分解的N/2个参数加参数a0f=(0:m-1)/(N*dt);a=zeros(1,m);b=zeros(1,m);%产生a，b两个为零的序列for k=0:m-1,for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N);b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);%MATLAB中的数组序号只能从1开始endc(k+1)=sqrt(a(k+1)^2+b(k+1)^2);endsubplot(211);plot(t,x);title('原始信号');xlabel('时间/s');subplot(212);plot(f,c);axis([0,20,0,2]);title('Fourier变换’)；xlabel('频率/Hz');ylabel('振幅')频率1，4,8,17HzFFTa=importdata('speed1.txt')t=a(1:1000,1)x=a(1:1000,2)N=1000;fs=1000subplot(211);plot(t,x);%进行FFT变换并做频谱图y=fft(x,N);%进行fft变换mag=abs(y);%求幅值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(212);plot(f,mag);%做频谱图title('正弦信号幅频谱图');axis([0,20,0,400]);频率1，4,8,17Hz改进clc;clear;a=importdata('speed1.txt')t=a(1:1000,1)x=a(1:1000,2)N=1000;fs=1000subplot(311);plot(t,x);title('原始信号')%进行FFT变换并做频谱图y=fft(x,N);%进行fft变换mag=abs(y);%fft模值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(312);plot(f,mag);%做频谱图title('fft模值—频率');axis([0,20,0,400]);yy=mag/(N/2);yy(1)=mag(1)/N;subplot(313);plot(f,yy);axis([0,20,0,2]);title('幅—频');角度形式clc;clear;a=importdata('speed1.txt') t=a(1:1000,1)x=a(1:1000,2)xd=x*360/60N=1000;fs=1000subplot(411);plot(t,x);title('原始信号') subplot(412);plot(t,xd);title('原始信号')%进行FFT变换并做频谱图y=fft(xd,N);%进行fft变换mag=abs(y);%fft模值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(413);plot(f,mag);%做频谱图title('fft模值—频率');axis([0,20,0,2000]);yy=mag/(N/2);yy(1)=mag(1)/N;subplot(414);plot(f,yy);axis([0,20,0,10]);title('幅—频');弧度形式clc;clear;a=importdata('speed1.txt') t=a(1:1000,1)x=a(1:1000,2)xr=x*2*pi/60N=1000;fs=1000subplot(411);plot(t,x);title('原始信号')subplot(412);plot(t,xr);title('原始信号')%进行FFT变换并做频谱图y=fft(xr,N);%进行fft变换mag=abs(y);%fft模值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(413);plot(f,mag);%做频谱图title('fft模值—频率');axis([0,20,0,50]);yy=mag/(N/2);yy(1)=mag(1)/N;subplot(414);plot(f,yy);axis([0,20,0,0.1]);title('幅—频');。
matlab部分代码
% tolerance, in which case only the "principal" radius will
% be picked up.
% tic;
% [accum, circen, cirrad] = CircularHough_Grd(rawimg, [15 60]);
% toc;
% figure(1); imagesc(accum); axis image;
% title('Accumulation Array from Circular Hough Transform');
matlab1源程序
program1:(1)
function [accum, varargout] = CircularHough_Grd(img, radrange, varargin)
%Detect circular shapes in a grayscale image. Resolve their center
% with dramatically different maximum intensities, e.g.
% 10-bit bitmaps in stead of the assumedOT be used. A value of 4% to 10% of the maximum
% Mask from the search of local maxima in the accumulation
% array.
%
%%%%%%%%% EXAMPLE #0:
数学建模算法的matlab代码
二，hamiton回路算法提供一种求解最优哈密尔顿的算法---三边交换调整法，要求在运行jiaohuan3（三交换法）之前，给定邻接矩阵C和节点个数N，结果路径存放于R中。
bianquan.m文件给出了一个参数实例，可在命令窗口中输入bianquan,得到邻接矩阵C和节点个数N以及一个任意给出的路径R,，回车后再输入jiaohuan3,得到了最优解。
由于没有经过大量的实验，又是近似算法，对于网络比较复杂的情况，可以尝试多运行几次jiaohuan3,看是否能到进一步的优化结果。
%%%%%%bianquan.m%%%%%%%N=13;for i=1:Nfor j=1:NC(i,j)=inf;endendfor i=1:NC(i,i)=0;endC(1,2)=6.0;C(1,13)=12.9;C(2,3)=5.9;C(2,4)=10.3;C(3,4)=12.2;C(3,5)=17.6;C(4,13)=8.8;C(4,7)=7.4;C(4,5)=11.5;C(5,2)=17.6;C(5,6)=8.2;C(6,9)=14.9;C(6,7)=20.3;C(7,9)=19.0;C(7,8)=7.3;C(8,9)=8.1;C(8,13)=9.2;C(9,10)=10.3;C(10,11)=7.7;C(11,12)=7.2;C(12,13)=7.9;for i=1:Nfor j=1:Nif C(i,j) < infC(j,i)=C(i,j);endendendfor i=1:NC(i,i)=0;endR=[4 7 6 5 3 2 1 13 12 11 10 9 8];<pre name="code" class="plain">%%%%%%%%jiaohuan3.m%%%%%%%%%%n=0;for I=1:(N-2)for J=(I+1):(N-1)for K=(J+1):Nn=n+1;Z(n,:)=[I J K];endendendR=1:Nfor m=1:(N*(N-1)*(N-2)/6)I=Z(m,1);J=Z(m,2);K=Z(m,3); r=R;if J-I~=1&K-J~=1&K-I~=N-1 for q=1:(J-I)r(I+q)=R(J+1-q);endfor q=1:(K-J)r(J+q)=R(K+1-q);endendif J-I==1&K-J==1r(K)=R(J);r(J)=R(K);endif J-I==1&K-J~=1&K-I~=N-1 for q=1:(K-J)r(I+q)=R(I+1+q); endr(K)=R(J);endif K-J==1&J-I~=1&K~=Nfor q=1:(J-I)r(I+1+q)=R(I+q); endr(I+1)=R(K);endif I==1&J==2&K==Nfor q=1:(N-2)r(1+q)=R(2+q);endr(N)=R(2);endif I==1&J==(N-1)&K==Nfor q=1:(N-2)r(q)=R(1+q);endr(N-1)=R(1);endif J-I~=1&K-I==N-1for q=1:(J-1)r(q)=R(1+q);endr(J)=R(1);endif J==(N-1)&K==N&J-I~=1r(J+1)=R(N);for q=1:(N-J-1)r(J+1+q)=R(J+q);endendif cost_sum(r,C,N)<cost_sum(R,C,N)R=rendendfprintf('总长为%f\n',cost_sum(R,C,N))%%%%%%cost_sum.m%%%%%%%%functiony=cost_sum(x,C,N)y=0;for i=1:(N-1)y=y+C(x(i),x(i+1));endy=y+C(x(N),x(1));三，灰色预测代码<pre name="code" class="plain">clearclcX=[136 143 165 152 165 181 204 272 319 491 571 605 665 640 628];x1(1)=X(1);X1=[];for i=1:1:14x1(i+1)=x1(i)+X(i+1);X1=[X1,x1(i)];endX1=[X1,X1(14)+X(15)]for k=3:1:15p(k)=X(k)/X1(k-1);p1(k)=X1(k)/X1(k-1);endp,p1clear kZ=[];for k=2:1:15z(k)=0.5*X1(k)+0.5*X1(k-1);Z=[Z,z(k)];endZB=[-Z',ones(14,1)]Y=[];clear ifor i=2:1:15Y=[Y;X(i)];endYA=inv(B'*B)*B'*Yclear ky1=[];for k=1:1:15y(k)=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1); y1=[y1;y(k)];endy1clear kX2=[];for k=2:1:15x2(k)=y1(k)-y1(k-1);X2=[X2;x2(k)];endX2=[y1(1);X2]e=X'-X2m=abs(e)./X's=e'*en=sum(m)/13clear ksyms ky=(X(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1)Y1=[];for j=16:1:21y11=subs(y,k,j)-subs(y,k,j-1);Y1=[Y1;y11];endY1%程序中的变量定义:alpha是包含α、μ值的矩阵;%ago是预测后累加值矩阵;var是预测值矩阵;%error是残差矩阵; c是后验差比值function basicgrey(x,m) %定义函数basicgray(x)if nargin==1 %m为想预测数据的个数，默认为1 m=1;endclc; %清屏,以使计算结果独立显示if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换x=x';endn=length(x); %取输入数据的样本量x1(:,1)=cumsum(x); %计算累加值,并将值赋及矩阵be for i=2:n %对原始数列平行移位 Y(i-1,:)=x(i,:);endfor i=2:n %计算数据矩阵B的第一列数据z(i,1)=0.5*x1(i-1,:)+0.5*x1(i,:);endB=ones(n-1,2); %构造数据矩阵BB(:,1)=-z(2:n,1);alpha=inv(B'*B)*B'*Y; %计算参数α、μ矩阵for i=1:n+m %计算数据估计值的累加数列,如改n+1为n+m可预测后m个值ago(i,:)=(x1(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1, :)*(i-1))+alpha(2,:)/alpha(1,:);endvar(1,:)=ago(1,:);f or i=1:n+m-1 %可预测后m个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原,并计算出下m个预测值end[P,c,error]=lcheck(x,var); %进行后验差检验[rela]=relations([x';var(1:n)']); %关联度检验ago %显示输出预测值的累加数列alpha %显示输出参数α、μ数列var %显示输出预测值error %显示输出误差P %显示计算小残差概率 c %显示后验差的比值crela %显示关联度judge(P,c,rela) %评价函数显示这个模型是否合格<pre name="code" class="plain">function judge(P,c,rela) %评价指标并显示比较结果if rela>0.6'根据经验关联度检验结果为满意（关联度只是参考主要看后验差的结果）'else'根据经验关联度检验结果为不满意（关联度只是参考主要看后验差的结果）'endif P>0.95&c<0.5'后验差结果显示这个模型评价为“优”'else if P>0.8&c<0.5'后验差结果显示这个模型评价为“合格”'else if P>0.7&c<0.65'后验差结果显示这个模型评价为“勉强合格”' else'后验差结果显示这个模型评价为“不合格”' endendendfunction [P,c,error]=lcheck(x,var)%进行后验差检验n=length(x);for i=1:nerror(i,:)=abs(var(i,:)-x(i,:)); %计算绝对残差c=std(abs(error))/std(x); %调用统计工具箱的标准差函数计算后验差的比值cs0=0.6745*std(x);ek=abs(error-mean(error));pk=0;for i=1:nif ek(i,:)<s0pk=pk+1;endendP=pk/n; %计算小残差概率%附带的质料里有一部分讲了关联度function [rela]=relations(x)%以x(1,:)的参考序列求关联度[m,n]=size(x);for i=1:mfor j=n:-1:2x(i,j)=x(i,j)/x(i,1);endfor i=2:mx(i,:)=abs(x(i,:)-x(1,:)); %求序列差endc=x(2:m,:);Max=max(max(c)); %求两极差Min=min(min(c));p=0.5; %p称为分辨率,0<p<1,一般取p=0.5for i=1:m-1for j=1:nr(i,j)=(Min+p*Max)/(c(i,j)+p*Max); %计算关联系数endendfor i=1:m-1rela(i)=sum(r(i,:))/n; %求关联度end四，非线性拟合function f=example1(c,tdata)f=c(1)*(exp(-c(2)*tdata)-exp(-c(3)*tdata));<pre name="code" class="plain">function f=zhengtai(c,x) f=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c( 2)^2));x=1:1:12;y=[01310128212]';c0=[2 8];for i=1:1000c=lsqcurvefit(@zhengtai,c0,x,y);c0=c;endy1=(1./(sqrt(2.*3.14).*c(1))).*exp(-(x-c(1)).^2./(2.*c (2)^2));plot(x,y,'r-',x,y1);legend('实验数据','拟合曲线')x=[0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16]';y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4]';f=@(c,x)c(1)*(exp(-c(2)*x)-exp(-c(3)*x));c0=[114 0.1 2]';for i=1:50opt=optimset('TolFun',1e-3);[c R]=nlinfit(x,y,f,c0,opt)c0=c;hold onplot(x,c(1)*(exp(-c(2)*x)-exp(-c(3)*x)),'g')endt=[0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16];y=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4];c0=[1 1 1];for i=1:50 c=lsqcurvefit(@example1,c0,t,y);c0=c;endy1=c(1)*(exp(-c(2)*t)-exp(-c(3)*t));plot(t,y,' +',t,y1);legend('实验数据','拟合曲线')五，插值拟合相关知识在生产和科学实验中，自变量及因变量间的函数关系有时不能写出解析表达式，而只能得到函数在若干点的函数值或导数值，或者表达式过于复杂而需要较大的计算量。
MATLAB源代码
埃特金插值function f=Atken(x,y,x0)syms t;if(length(x)==length(y))n=length(x);elsedisp;return;endy1(1:n)=t;for(i=1:n-1)for(j=i+1:n)y1(j)=y(j)*(t-x(i))/(x(j)-x(i))+y(i)*(t-x(j))/(x(i)-x(j));endy=y1simplify(y1);endif(nargin==3)f=subs(y1(n),'t',x0);elsesimplify(y1(n));f=collect(y1(n));f=vpa(f,6);end特征多项式法function l=Chapoly(A)%特征多项式法求矩阵特征值%已知矩阵：A%求得的矩阵特征值：lsyms t;N=size(A);n=N(1,1);y=det(A-t*eye (n,n));l=solve(y);l=vpa(1,5); %function f = Chebyshev(y,k,x0)syms t;T(1:k+1) = t;T(1) = 1;T(2) = t;c(1:k+1) = 0.0;c(1)=int(subs(y,findsym(sym(y)),sym('t'))*T(1)/sqrt(1-t^2),t,-1,1)/pi; c(2)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(2)/sqrt(1-t^2),t,-1,1)/pi;f = c(1)+c(2)*t;for i=3:k+1T(i) = 2*t*T(i-1)-T(i-2);c(i) = 2*int(subs(y,findsym(sym(y)),sym('t'))*T(i)/sqrt(1-t^2),t,-1,1)/2;f = f + c(i)*T(i);f = vpa(f,6);if(i==k+1)if(nargin == 3)f = subs(f,'t',x0);elsef = vpa(f,6);endendend弦割法求解非线性方程function [x k t]=ChordsecantToEquation(f,x0,x1,eps)%弦割法求解非线性方程%[x k t]=ChordsecantToEquation(f,x0,x1,eps)%x:近似解%k:迭代次数%t:运算时间%f:原函数，定义为内联函数%x0,x1:初始值%eps:误差限%%应用举例：%f=inline('x^3+4*x^2-10');%x=ChordsecantToEquation(f,1,2,0.5e-6)%[x k]=ChordsecantToEquation(f,1,2,0.5e-6)%[x k t]=ChordsecantToEquation(f,1,2,0.5e-6)%函数的最后一个参数也可以不写，默认情况下，eps=0.5e-6%[x k t]=ChordsecantToEquation(f,1,2)if nargin==3eps=0.5e-6;endk=0;while 1x=x1-f(x1)*(x1-x0)./(f(x1)-f(x0));k=k+1;if abs(x-x1) < eps || k > 30break;endx0=x1;x1=x;endt=toc;if k >= 30disp('迭代次数太多。
Matlab程序代码
Matlab程序代码function [seg]=MRMRF(w,class_number,potential,maxIter)%MRMRF图像分割算法%w-待分解图像的多尺度序列%class_number-分类数%potential-potts模型势函数%maxIter-最大迭代次数%多分辨率表达的尺度数L=size(w,1);%分割结果的多尺度序列seg=cell(L,1);%使用ICM算法计算最高尺度的分割结果seg{L}=ICM(w{L},class_number,potential,maxIter);%计算其他尺度的分割结果for n=(L-1):-1:1%获取初始分割结果segn=myZoomOut(seg{n+1});%获取尺度n上的最终分割结果wn=w{n};segn=ICMn(wn,segn,class_number,potential,maxIter);%保存尺度n上的分割结果seg{n}=segn;endendfunction [seg]=ICMn(image,seg,class_number,potential,maxIter) %尺度n上已知初始结果的ICM分割[width,height,bands]=size(image);%将图像和初始分割结果分别转换为向量image=imstack2vectors(image);seg=imstack2vectors(seg);%ICM迭代iter=0;while(iter<maxIter)%估计特征场参数[mu,sigma]=GMM_parameter(image,seg,class_number);%计算特征场能量Ef=EnergyOfFeatureField(image,mu,sigma,class_number);%计算标记场能量El=EnergyOfLabelField(seg,potential,width,height,class_number);% 计算新的分割结果E=Ef+El;[tm,seg]=min(E,[],2);% 更新迭代条件iter=iter+1;endseg=reshape(seg,[width height]);endfunction [x]=myZoomOut(x)%放大x两倍[col,raw,B]=size(x);tm=zeros(2*col,2*raw,B);tm(1:2:2*col,1:2:2*raw,:)=x;tm(2:2:2*col,1:2:2*raw,:)=x;tm(1:2:2*col,2:2:2*raw,:)=x;tm(2:2:2*col,2:2:2*raw,:)=x;x=tm;end%高斯模型参数估计function [mu,sigma]=GMM_parameter(image,seg,class_number)% image-已经向量化的图像数据，每一行表示一个像素% seg-与image对应的图像标记，每一个对应一个像素标记% class_number-图像的分类数[n,d]=size(image);mu=zeros(class_number,d);sigma=zeros(d,d,class_number);for i=1:class_numberIm_i=image(seg==i,:);[sigma(:,:,i),mu(i,:)]=covmatrix(Im_i);endendfunction [C,m]=covmatrix(X)[k,n]=size(X);X=double(X);if k==1C=0;m=X;elsem=sum(X,1)/k;X=X-m(ones(k,1),:);C=(X'*X)/(k-1);m=m';endendfunction [E]=EnergyOfLabelField(seg,potential,width,height,class_number) % seg-向量化的像素标记矩阵% potential-Potts模型势函数取值% width,height-图像大小n=size(seg,1);seg=reshape(seg,[width height]);%计算领域nei8_image=neighbor8syn(seg);Nei8=imstack2vectors(nei8_image);%计算标记场能量E=zeros(n,class_number);for i=1:class_numberE(:,i)=sum(Nei8~=i,2);endE=E*potential;endfunction nei=neighbor8syn(Y)[s,t,P]=size(Y);Yul=zeros(s,t,P);Yul(2:s,2:t,:)=Y(1:s-1,1:t-1,:);Yu=zeros(s,t,P);Yu(2:s,:,:)=Y(1:s-1,:,:);Yur=zeros(s,t,P);Yur(2:s,1:t-1,:)=Y(1:s-1,2:t,:);Yr=zeros(s,t,P);Yr(:,1:t-1,:)=Y(:,2:t,:);Ydr=zeros(s,t,P);Ydr(1:s-1,1:t-1,:)=Y(2:s,2:t,:);Yd=zeros(s,t,P);Yd(1:s-1,:,:)=Y(2:s,:,:);Ydl=zeros(s,t,P);Ydl(1:s-1,2:t,:)=Y(2:s,1:t-1,:);Yl=zeros(s,t,P);Yl(:,2:t,:)=Y(:,1:t-1,:);nei=cat(3,Yul,Yu,Yur,Yr,Ydr,Yd,Ydl,Yl);endfunction [E]=EnergyOfFeatureField(image,mu,sigma,class_number) % mu-均值矩阵% sigma-方差矩阵n=size(image,1);E=zeros(n,class_number);% 计算能量for i=1:class_numbermu_i=mu(i,:);sigma_i=sigma(:,:,i);diff_i=image-repmat(mu_i,[n,1]);E(:,i)=sum(diff_i*inv(sigma_i).*diff_i,2)+log(det(sigma_i)); endend。
Matlab源程序代码
正弦波的源程序：(一)，用到的函数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;end2,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 Nclose allk=input('取样点数＝2^k, k取10摆布');if isempty(k), k=10; endf0=input('f0=取1(kz)摆布');if isempty(f0), f0=1; endN=2^k;dt=0.01; %msdf=1/(N*dt); %KHzT=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');gridaxis([-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') %输出信号波形画图gridaxis([-2/f0,+2/f0,-1.5,1.5])xlabel('t(ms)')ylabel('a(t)(V)')gtext('频谱图')最佳基带系统的源程序：(一)，用到的函数f2t函数和t2f函数。
MATLAB 代码
嵌入：M=256;%原图像长度N=32;%水印图像长度K=8;I=zeros(M,M);J=zeros(N,N);BLOCK=zeros(K,K);%显示水印图像subplot(1,8,2);J=imread('14','bmp');imshow(J);title('水印图像');%显示原图像subplot(1,3,2);I=imread('11','bmp');imshow(I);title('原始公开图像');%嵌入水印tem=1;for p=1:Nfor q=1:Nx=(p-1)*K+1;y=(q-1)*K+1;BLOCK=I(x:x+K-1,y:y+K-1); BLOCK=dct2(BLOCK);if J(p,q)==0a=-1;elsea=1;endBLOCK(2,1)=BLOCK(2,1)*(1+a*0.01); BLOCK=idct2(BLOCK);I(x:x+K-1,y:y+K-1)=BLOCK;endend%显示嵌入水印后的图像subplot(1,3,3);imshow(I);title('嵌入水印后的图像');imwrite(I,'embedded.bmp','bmp');基于小波变换的数字图像水印处理（MATLAB源代码）clear all; close all; clc;M=256;%原图像长度N=64; %水印长度[filename1,pathname]=uigetfile('*.*','select the image');image1=imread(num2str(filename1));subplot(2,2,1);imshow(image1); title('original image'); % orginal image for watermarkingimage1=double(image1);imagew=imread('dmg2.tif');subplot(2,2,2);imshow(imagew);title('original watermark');%original watermark%嵌入水印[ca,ch,cv,cd] = dwt2(image1,'db1');[cas,chs,cvs,cds] = dwt2(ca,'db1');for i=1:Nfor j=1:Nif imagew(i,j)==0a=-1;elsea=1;endCa(i,j)=cas(i,j)*(1+a*0.03);endendIM= idwt2(Ca,chs,cvs,cds,'db1') ;markedimage=double(idwt2(IM,ch,cv,cd,'db1'));%显示嵌入后水印图像subplot(2,2,3);colormap(gray(256));image(markedimage);title('mark ed image');imwrite(markedimage,gray(256),'watermarked.bmp','bmp');%提取水印image1=imread(num2str(filename1));image1=double(image1);imaged=imread('watermarked.bmp');[ca,ch,cv,cd] = dwt2(image1,'db1');[cas,chs,cvs,cds]=dwt2(ca,'db1');[caa,chh,cvv,cdd]=dwt2(imaged,'db1');[caas,chhs,cvvs,cdds]=dwt2(caa,'db1');for p=1:Nfor q=1:Na=caas(p,q)/cas(p,q)-1;if a<0W(p,q)=0;elseW(p,q)=255;endendend%显示提取的水印subplot(2,2,4);colormap(gray(256));image(W);title('从含水印图像中提取的水印'); imwrite(W,gray(256),'watermark.bmp','bmp');MAIN.m 主程序%-------------------水印嵌入------------------------------------------------while 1clear;c=0.3;a=imread('lina.jpg');%原图像b=imread('changsha.bmp')*255;%二值水印图像[m1,n1]=size(a);[m2,n2]=size(b);e0=(sum(sum(a.^2)))/(m1*n1);e0=c*e0;%计算基准能量[t,tkey]=sdwt(double(a),'db2',m2,n2,e0);%树状小波分解[t,tkey]=embed(t,tkey,b);%嵌入水印aw=sidwt(t,'db2');%重组figure(1);subplot(1,2,1);imshow(uint8(a));title('原图');subplot(1,2,2);imshow(uint8(aw));title('嵌入图');imwrite(uint8(aw),'watermark.jpg');% csvwrite('key.txt',reshape(tkey,m2,n2));v1=m1*m1*255*255;v2=sum(sum((double(a)-aw).^2));snr=10*log10(v1/v2);% 峰值信噪比snr。
（完整word版）matlab经典代码大全
（完整word版）matlab经典代码大全哈哈哈MATLAB显示正炫余炫图：plot(x,y1,'* r',x,y2,'o b')定义【0,2π】；t=0:pi/10:2*pi;定义函数文件：function [返回变量列表]=函数名（输入变量列表）顺序结构：选择结构1）if-else-end语句其格式为：if 逻辑表达式程序模块1；else程序模块2；End图片读取：%选择图片路径[filename, pathname] = ...uigetfile({'*.jpg';'*.bmp';'*.gif'},'选择图片');%合成路径+文件名str=[pathname,filename];%为什么pathname和filename要前面出现的位置相反才能运行呢%读取图片im=imread(str);%使用图片axes(handles.axes1);%显示图片imshow(im);边缘检测：global imstr=get(hObject,'string');axes (handles.axes1);switch strcase ' 原图'imshow(im);case 'sobel'BW = edge(rgb2gray(im),'sobel');imshow(BW);case 'prewitt'BW = edge(rgb2gray(im),'prewitt');imshow(BW);case 'canny'BW = edge(rgb2gray(im),'canny');imshow(BW);Canny算子边缘定位精确性和抗噪声能力效果较好，是一个折中方案end;开闭运算：se=[1,1,1;1,1,1;1,1,1;1,1,1]; %Structuring ElementI=rgb2gray(im);imshow(I,[]);title('Original Image');I=double(I);[im_height,im_width]=size(I);[se_height,se_width]=size(se);halfheight=floor(se_height/2);halfwidth=floor(se_width/2);[se_origin]=floor((size(se)+1)/2);image_dilation=padarray(I,se_origin,0,'both'); %Image to be used for dilationimage_erosion=padarray(I,se_origin,256,'both'); %Image to be used for erosion %%%%%%%%%%%%%%%%%% %%% Dilation %%%%%%%%%%%%%%%%%%%%%for k=se_origin(1)+1:im_height+se_origin(1)for kk=se_origin(2)+1:im_width+se_origin(2)dilated_image(k-se_origin(1),kk-se_origin(2))=max(max(se+image_dilation(k-se_origin(1):k+halfh eight-1,kk-se_origin(2):kk+halfwidth-1)));endendfigure;imshow(dilated_image,[]);title('Image after Dilation'); %%%%%%%%%%%%%%%%%%%% Erosion %%%%%%%%%%%%%%%%%%%%se=se';for k=se_origin(2)+1:im_height+se_origin(2)for kk=se_origin(1)+1:im_width+se_origin(1)eroded_image(k-se_origin(2),kk-se_origin(1))=min(min(image_erosion(k-se_origin(2):k+halfwidth -1,kk-se_origin(1):kk+halfheight-1)-se));endendfigure;imshow(eroded_image,[]);title('Image after Erosion'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%% Opening(Erosion first, then Dilation) %%% %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%se=se';image_dilation2=eroded_image; %Image to be used for dilationfor k=se_origin(1)+1:im_height-se_origin(1)for kk=se_origin(2)+1:im_width-se_origin(2)opening_image(k-se_origin(1),kk-se_origin(2))=max(max(se+image_dilation2(k-se_origin(1):k+hal fheight-1,kk-se_origin(2):kk+halfwidth-1)));endendfigure;imshow(opening_image,[]);title('OpeningImage'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% Closing(Dilation first, then Erosion) %%% %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%se=se';image_erosion2=dilated_image; %Image to be used for erosionfor k=se_origin(2)+1:im_height-se_origin(2)for kk=se_origin(1)+1:im_width-se_origin(1)closing_image(k-se_origin(2),kk-se_origin(1))=min(min(image_erosion2(k-se_origin(2):k+halfwidt h-1,kk-se_origin(1):kk+halfheight-1)-se));endendfigure;imshow(closing_image,[]);title('Closing Image');Warning: Image is too big to fit on screen; displaying at 31% scale.> In truesize>Resize1 at 308In truesize at 44In imshow at 161图像的直方图归一化：I=imread(‘red.bmp’);%读入图像figure;%打开新窗口[M,N]=size(I);%计算图像大小[counts,x]=imhist(I,32);%计算有32个小区间的灰度直方图counts=counts/M/N;%计算归一化灰度直方图各区间的值stem(x,counts);%绘制归一化直方图图像平移：I=imread('shuichi.jpg');se=translate(strel(1),[180 190]);B=imdilate(I,se);figure;subplot(1,2,1),subimage(I);title('原图像');subplot(1,2,2),subimage(B);title('平移后图像');图像的转置；A=imread('nir.bmp');tform=maketform('affine',[0 1 0;1 0 0;0 0 1]);B=imtransform(A,tform,'nearest');figure;imshow(A);figure;imshow(B);imwrite(B,'nir转置后图像.bmp');图像滤波：B = imfilter(A,H,option1,option2,...)或写作g = imfilter(f, w, filtering_mode, boundary_options, size_options)其中，f为输入图像，w为滤波掩模，g为滤波后图像。
matlab函数代码
matlab函数代码Matlab是一种高级的数学计算软件，它可以用于各种科学和工程领域的计算和数据分析。
Matlab的强大之处在于它具有丰富的函数库和易于使用的编程语言，使得用户可以快速地编写复杂的数学计算和数据分析程序。
下面是一些常用的Matlab函数代码。
1. linspace函数linspace函数用于生成一组等间隔的数值。
它的语法如下：x = linspace(start, stop, n)其中start和stop是生成数值的起始值和结束值，n是生成数值的个数。
例如，要生成从0到1之间的10个等间隔的数值，可以使用以下代码：x = linspace(0, 1, 10)2. plot函数plot函数用于绘制二维图形。
它的语法如下：plot(x, y)其中x和y是要绘制的数据点的坐标。
例如，要绘制y = sin(x)的图形，可以使用以下代码：x = linspace(0, 2*pi, 100);y = sin(x);plot(x, y)3. rand函数rand函数用于生成随机数。
它的语法如下：x = rand(n)其中n是要生成的随机数的个数。
例如，要生成10个随机数，可以使用以下代码：x = rand(10)4. mean函数mean函数用于计算一组数据的平均值。
它的语法如下：m = mean(x)其中x是要计算平均值的数据。
例如，要计算向量x的平均值，可以使用以下代码：x = [1, 2, 3, 4, 5];m = mean(x)5. std函数std函数用于计算一组数据的标准差。
它的语法如下：s = std(x)其中x是要计算标准差的数据。
例如，要计算向量x的标准差，可以使用以下代码：x = [1, 2, 3, 4, 5];s = std(x)6. eig函数eig函数用于计算矩阵的特征值和特征向量。
它的语法如下：[V, D] = eig(A)其中A是要计算特征值和特征向量的矩阵，V是特征向量的矩阵，D 是特征值的对角矩阵。
matlab源代码
例错误！文档中没有指定样式的文字。
-1%周期信号（方波）的展开,fb_jinshi.mclose all;clear all;N=100; %取展开式的项数为2N＋1项T=1;fs=1/T;N_sample=128; %为了画出波形，设置每个周期的采样点数dt = T/N_sample;t=0:dt:10*T-dt;n=-N:N;Fn = sinc(n/2).*exp(-j*n*pi/2);Fn(N+1)=0;ft = zeros(1,length(t));for m=-N:Nft = ft + Fn(m+N+1)*exp(j*2*pi*m*fs*t);endplot(t,ft)例错误！文档中没有指定样式的文字。
-4利用FFT计算信号的频谱并与信号的真实频谱的抽样比较。
脚本文件T2F.m定义了函数T2F，计算信号的傅立叶变换。
function [f,sf]= T2F(t,st)%This is a function using the FFT function to calculate a signal's Fourier %Translation%Input is the time and the signal vectors,the length of time must greater %than 2%Output is the frequency and the signal spectrumdt = t(2)-t(1);T=t(end);df = 1/T;N = length(st);f=-N/2*df:df:N/2*df-df;sf = fft(st);sf = T/N*fftshift(sf);脚本文件F2T.m定义了函数F2T，计算信号的反傅立叶变换。
function [t st]=F2T(f,sf)%This function calculate the time signal using ifft function for the input %signal's spectrumdf = f(2)-f(1);Fmx = ( f(end)-f(1) +df);dt = 1/Fmx;N = length(sf);T = dt*N;%t=-T/2:dt:T/2-dt;t = 0:dt:T-dt;sff = fftshift(sf);st = Fmx*ifft(sff);另写脚本文件fb_spec.m如下：%方波的傅氏变换, fb_spec.mclear all;close all;T=1;N_sample = 128;dt=T/N_sample;t=0:dt:T-dt;st=[ones(1,N_sample/2), -ones(1,N_sample/2)]; %方波一个周期subplot(211);plot(t,st);axis([0 1 -2 2]);xlabel('t'); ylabel('s(t)');subplot(212);[f sf]=T2F(t,st); %方波频谱plot(f,abs(sf)); hold on;axis([-10 10 0 1]);xlabel('f');ylabel('|S(f)|');%根据傅氏变换计算得到的信号频谱相应位置的抽样值sff= T^2*j*pi*f*0.5.*exp(-j*2*pi*f*T).*sinc(f*T*0.5).*sinc(f*T*0.5);plot(f,abs(sff),'r-')例错误！文档中没有指定样式的文字。
matlab相关代码
matlab相关代码Matlab是一种高级技术计算语言和交互式环境，广泛应用于科学、工程和金融等领域。
它具有强大的数学和图形处理能力，可以快速处理大量数据和进行复杂的算法分析。
以下是一些常用的Matlab相关代码：1. 矩阵操作Matlab中的矩阵操作非常方便，可以使用以下代码实现：a = [1 2 3; 4 5 6; 7 8 9]; % 定义一个3x3的矩阵b = [4 5 6; 7 8 9; 10 11 12]; % 定义另一个3x3的矩阵c = a + b; % 矩阵加法d = a - b; % 矩阵减法e = a * b; % 矩阵乘法f = a .* b; % 矩阵对应元素相乘g = a' % 矩阵转置2. 绘图Matlab中的绘图功能非常强大，可以使用以下代码实现：x = 0:0.1:2*pi; % 定义x轴的范围y = sin(x); % 定义y轴的值plot(x,y); % 绘制sin函数图像xlabel('x'); % 设置x轴标签ylabel('y'); % 设置y轴标签title('sin函数图像'); % 设置图像标题grid on; % 显示网格线3. 数据处理Matlab中的数据处理功能非常强大，可以使用以下代码实现：data = [1 2 3 4 5 6 7 8 9]; % 定义一个数据数组mean_data = mean(data); % 计算平均值std_data = std(data); % 计算标准差max_data = max(data); % 计算最大值min_data = min(data); % 计算最小值median_data = median(data); % 计算中位数4. 机器学习Matlab中的机器学习功能非常强大，可以使用以下代码实现：load fisheriris; % 加载鸢尾花数据集X = meas(:,1:2); % 取前两个特征Y = species; % 取标签svmModel = fitcsvm(X,Y); % 训练SVM模型svmModel.ClassNames % 显示类别名称svmModel.Beta % 显示模型系数总之，Matlab是一种非常强大的工具，可以帮助科学家、工程师和金融分析师等快速处理数据和进行复杂的算法分析。
实用matlab代码
实⽤matlab代码1.图像反转MATLAB程序实现如下：I=imread('xian.bmp');J=double(I);J=-J+(256-1); %图像反转线性变换H=uint8(J);subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(H);2.灰度线性变换MATLAB程序实现如下：I=imread('xian.bmp');subplot(2,2,1),imshow(I);title('原始图像');axis([50,250,50,200]);axis on; %显⽰坐标系I1=rgb2gray(I);subplot(2,2,2),imshow(I1);title('灰度图像');axis([50,250,50,200]);axis on; %显⽰坐标系J=imadjust(I1,[0.1 0.5],[]); %局部拉伸，把[0.1 0.5]内的灰度拉伸为[0 1] subplot(2,2,3),imshow(J); title('线性变换图像[0.1 0.5]');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系K=imadjust(I1,[0.3 0.7],[]); %局部拉伸，把[0.3 0.7]内的灰度拉伸为[0 1] subplot(2,2,4),imshow(K); title('线性变换图像[0.3 0.7]');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系3.⾮线性变换MATLAB程序实现如下：I=imread('xian.bmp');I1=rgb2gray(I);subplot(1,2,1),imshow(I1);title('灰度图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系J=double(I1);J=40*(log(J+1));H=uint8(J);subplot(1,2,2),imshow(H);title('对数变换图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系4.直⽅图均衡化MATLAB程序实现如下：I=imread('xian.bmp');I=rgb2gray(I);figure;subplot(2,2,1);imshow(I);subplot(2,2,2);imhist(I);I1=histeq(I);figure;subplot(2,2,1);imshow(I1);subplot(2,2,2);imhist(I1);5.线性平滑滤波器⽤MATLAB实现领域平均法抑制噪声程序：I=imread('xian.bmp');subplot(231)imshow(I)title('原始图像')I=rgb2gray(I);I1=imnoise(I,'salt & pepper',0.02);subplot(232)imshow(I1)title('添加椒盐噪声的图像')k1=filter2(fspecial('average',3),I1)/255; %进⾏3*3模板平滑滤波k2=filter2(fspecial('average',5),I1)/255; %进⾏5*5模板平滑滤波k3=filter2(fspecial('average',7),I1)/255; %进⾏7*7模板平滑滤波k4=filter2(fspecial('average',9),I1)/255; %进⾏9*9模板平滑滤波subplot(233),imshow(k1);title('3*3模板平滑滤波');subplot(234),imshow(k2);title('5*5模板平滑滤波');subplot(235),imshow(k3);title('7*7模板平滑滤波');subplot(236),imshow(k4);title('9*9模板平滑滤波');6.中值滤波器⽤MATLAB实现中值滤波程序如下：I=imread('xian.bmp');I=rgb2gray(I);J=imnoise(I,'salt&pepper',0.02);subplot(231),imshow(I);title('原图像');subplot(232),imshow(J);title('添加椒盐噪声图像');k1=medfilt2(J); %进⾏3*3模板中值滤波k2=medfilt2(J,[5,5]); %进⾏5*5模板中值滤波k3=medfilt2(J,[7,7]); %进⾏7*7模板中值滤波k4=medfilt2(J,[9,9]); %进⾏9*9模板中值滤波subplot(233),imshow(k1);title('3*3模板中值滤波'); subplot(234),imshow(k2);title('5*5模板中值滤波');subplot(235),imshow(k3);title('7*7模板中值滤波');subplot(236),imshow(k4);title('9*9模板中值滤波');7.⽤Sobel算⼦和拉普拉斯对图像锐化：I=imread('xian.bmp');subplot(2,2,1),imshow(I);title('原始图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系I1=im2bw(I);subplot(2,2,2),imshow(I1);title('⼆值图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系H=fspecial('sobel'); %选择sobel算⼦J=filter2(H,I1); %卷积运算subplot(2,2,3),imshow(J);title('sobel算⼦锐化图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系h=[0 1 0,1 -4 1,0 1 0]; %拉普拉斯算⼦J1=conv2(I1,h,'same'); %卷积运算subplot(2,2,4),imshow(J1); title('拉普拉斯算⼦锐化图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系8.梯度算⼦检测边缘⽤MATLAB实现如下：I=imread('xian.bmp');subplot(2,3,1);imshow(I);title('原始图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系I1=im2bw(I); subplot(2,3,2);imshow(I1);title('⼆值图像');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系I2=edge(I1,'roberts'); figure;subplot(2,3,3);imshow(I2);title('roberts算⼦分割结果');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系I3=edge(I1,'sobel'); subplot(2,3,4);imshow(I3);title('sobel算⼦分割结果');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系I4=edge(I1,'Prewitt'); subplot(2,3,5);imshow(I4);title('Prewitt算⼦分割结果');axis([50,250,50,200]);grid on; %显⽰⽹格线axis on; %显⽰坐标系9.LOG算⼦检测边缘⽤MATLAB程序实现如下：I=imread('xian.bmp');subplot(2,2,1);imshow(I);。
matlab代码
n=7:12;y=1./abs(n-6);plo t(n,y,'r*','Mar kersi ze',20)gr id on//画图x=l inspa ce(0,2*pi,100);y=si n(x);plot(x,y);加注解x=l inspa ce(0,2*pi,100);>> y=sin(x);>> plo t(x,y);>> x=li nspac e(0,2*pi,100);>> y=sin(x);>> plot(x,y,'r',x,y,'*');>> leg end('y=sin(x)')>> t ext(1.4,0.91,'波峰')x label('Inp ut Va lue');>>ylabe l('Fu nctio n Val ue');>> t itle('一个正弦函数');//a xis命令控制坐标轴x=l inspa ce(0,2*pi,100);y=si n(x);z=co s(x);plot(x,y,x,z);axis on>> axi s('sq uare')>>axis('equa l')//plo t函数t=(0:p i/100:pi);y1=s in(t)*[1,-1];y2=sin(t).*sin(9*t);t3=pi*(0:9)/9;y3=si n(t3).*sin(9*t3);pl ot(t,y1,'r:',t,y2,'b',t3,y3,'b o')a xis([0,pi,-1,1])//hold on 命令不会取消之前画的曲线>> hold on>> y4=-1:0.2:1;>> pl ot(1.5,y4,'*')//同一窗口用不同的坐标系绘制几组数据subp lot(2,2,1);plot(x,si n(x));>>subpl ot(2,2,2);plot(x,cos(x));>> s ubplo t(2,2,3);p lot(x,sinh(x));>> s ubplo t(2,2,4);p lot(x,cosh(x));//f igure函数x=lins pace(0,2*p i,30);>>figur e>>plot(x,sin(x));>> x label('Inp ut Va lue');>>ylabe l('Fu nctio n Val ue');>> t itle('一个正弦函数');>> l egend('y=s in(x)')>> figu re>> plot(x,si nh(x));>> clos e//直方图close all>> x=1:10;>> y=rand(size(x));>> b ar(x,y);>> x=1:10x=lin space(0,2,20)*p iy=s in(x)e=st d(y)*ones(size(x));error bar(x,y,e)//柄图x=l inspa ce(0,10,50);>> y=si n(x).*exp(-x/3);>>stem(x,y)//阶梯图> x=lin space(0,10,50);y=si n(x).*exp(-x/3);>>stair s(x,y);//饼图x=[11.4,23.5,35.4,15.6];>> ex plode=zero s(siz e(x));>>[c,of fset]=min(x);>> exp lode(offse t)=1;>> p ie(x,explo de);//频数累计柱状图x=ra ndn(5000,1);>> hist(x,20);//矢量图//羽状图thet a=90:-10:0;>>r=one s(siz e(the ta));>> [u,v]=pol2c art(t heta*pi/180,r*10);>> fea ther(u,v)>> ax is eq ual//极坐标图th eta=l inspa ce(0,2*pi);>>r=cos(4*th eta);>> p olar(theta,r);//精确绘图fp lot('sin(1./x)',[0.02,0.2]); 0.02,0.2是绘图范围三维绘图三维线性图形>> t=(0:0.02:2)*pi;>> x=sin(t);>> y=c os(t);>>z=cos(2*t);>>plot3(x,y,z,'b-',x,y,z,'b d');>> vi ew([-82,58]);>> box on;>> le gend('链','宝石');//三维曲面x=-4:0.2:4;>>y=x;>> [X,Y]=m eshgr id(x,y);>> Z=X.^2/9+Y.^2/9;>> sub plot(1,2,1);>> mesh(X,Y,Z)>> titl e('椭圆抛物面网线图')>> sub plot(1,2,2);>> surf(X,Y,Z)>> tit le('椭圆抛物面网面图')//设置视角vi ew(az,el):az:方位角，el:仰视角v iew([x,y,z]):由坐标指定观察点，即确定视角。
(完整word版)MatLab代码大全
第2章图像获取2.3.2 二维连续傅里叶变换例2.2figure(1); %建立图形窗口1[u,v] = meshgrid(-1:0.01:1); %生成二维频域网格F1 = abs(sinc(u.*pi));F2 = abs(sinc(v.*pi));F=F1.*F2; %计算幅度频谱F=|F(u,v)|surf(u,v,F); %显示幅度频谱，如图2.3(b)shading interp; %平滑三维曲面上的小格axis off; %关闭坐标系figure(2); %建立图形窗口2F1=histeq(F); %扩展F的对比度以增强视觉效果imshow(F1); %用图像来显示幅度频谱，如图2.3(c)第3章图像变换3.4.4 二维FFT的MATLAB实现例3.2 简单图像及其傅里叶变换MATLAB程序：%建立简单图像d并显示之d = zeros(32,32); %图像大小32⨯32d(13:20,13:20) = 1; %中心白色方块大小为8⨯8figure(1); %建立图形窗口1imshow(d,'notruesize')；%显示图像d如图3.5(a)所示%计算傅里叶变换并显示之D = fft2(d); %计算图像d的傅里叶变换，fft2(d) = fft(fft(d).').'figure(2); %建立图形窗口2imshow(abs(D),[-1 5],'notruesize'); %显示图像d的傅里叶变换谱如3.5(b)所示例3.3 MATLAB图像及其傅里叶变换谱MATLAB程序：figure(1);load imdemos saturn2; %装入MA TLAB图像saturn2imshow(saturn2); %显示图像saturn2如图3.6(a)所示figure(2);S= fftshift(fft2(saturn2)); %计算傅里叶变换并移位imshow(log(abs(S)),[ ]); %显示傅里叶变换谱如3.6(b)所示例3.4 真彩图像及其傅里叶变换谱MATLAB程序：figure(1);A=imread('image1.jpg'); %装入真彩图像，见图1.1(b)B=rgb2gray(A); %将真彩图像转换为灰度图像imshow(B); %显示灰度图像如图3.7(a)所示C=fftshift(fft2(B)); %计算傅里叶变换并移位figure(2);imshow(log(abs(C)),[ ]); %显示傅里叶变换谱如3.7(b)所示3.5.4 离散余弦变换的MATLAB实现例3.5 计算并显示真彩图像余弦变换的MATLAB程序如下：RGB=imread('image2.jpg'); %装入真彩图像figure(1);imshow(RGB); %显示彩色图像GRAY=rgb2gray(RGB); %将真彩图像转换为灰度图像figure(2);imshow(GRAY); %显示灰度图像如图3.10(a)所示DCT=dct2(GRAY); %进行余弦变换figure(3);imshow(log(abs(DCT)),[ ]); %显示余弦变换如图3.10(b)所示。
matlab相关代码
matlab相关代码Matlab代码实现图像的二值化图像二值化是图像处理中的一种基本操作，它将灰度图像转换为二值图像，使得图像中的像素只有两种取值，即黑色和白色。
在Matlab中，可以使用im2bw函数实现图像的二值化。
im2bw函数的语法格式为：BW = im2bw(I, level)其中，I为输入的灰度图像，level为二值化的阈值。
函数返回值BW为二值图像。
下面是一个简单的Matlab代码实现图像的二值化：% 读取灰度图像I = imread('lena_gray.jpg');% 显示原始图像subplot(1,2,1);imshow(I);title('原始图像');% 对图像进行二值化level = graythresh(I);BW = im2bw(I, level);% 显示二值化后的图像subplot(1,2,2);imshow(BW);title('二值化后的图像');上述代码中，首先使用imread函数读取灰度图像lena_gray.jpg，并使用subplot函数将原始图像和二值化后的图像显示在同一窗口中。
然后，使用graythresh函数计算二值化的阈值level，并使用im2bw函数将灰度图像I转换为二值图像BW。
最后，使用imshow 函数显示二值化后的图像。
需要注意的是，im2bw函数默认将灰度值大于等于阈值的像素设置为白色，将灰度值小于阈值的像素设置为黑色。
如果需要反转二值图像的颜色，可以使用imcomplement函数。
Matlab提供了丰富的图像处理函数，可以方便地实现图像的二值化、滤波、边缘检测等操作。
熟练掌握这些函数，可以大大提高图像处理的效率和精度。
matlab 内置的c代码
MATLAB 是一种高级编程语言和交互式环境，主要用于算法开发、数据可视化、数据分析以及数值计算。
它也有内置的 C/C++ 函数，可以在 MATLAB 中直接调用。
以下是一些 MATLAB 中内置的 C/C++ 函数的例子：1. `mexFunction`：这是创建 MATLAB 可执行文件（.mex 文件）所必需的入口函数。
它允许您在 MATLAB 中直接调用 C/C++ 代码。
```c#include "mex.h"void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {// 在这里写你的 C/C++ 代码}```2. `mxCreateDoubleMatrix`：这个函数用于创建一个双精度浮点数矩阵。
```cmxArray *mxCreateDoubleMatrix(int m, int n, mxComplexityComplexFlag complexFlag);```3. `mxGetPr` 和 `mxSetPr`：这两个函数用于获取和设置矩阵的原始指针。
```cdouble *mxGetPr(const mxArray *a);void mxSetPr(mxArray *a, const double *pr);```4. `mexCallMATLAB`：这个函数用于在 C/C++ 代码中调用 MATLAB 函数。
```cvoid mexCallMATLAB(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[], const char *name);```以上都是 MATLAB 提供的 C/C++ API，你可以在 MATLAB 的官方文档中找到更多关于这些 API 的详细信息。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Rbf神经网络预测：clcclose allclear%ld = 100; % the number of learning data %x = rand(2,200); % generate data(0--1)x=[2.1000 0.20003.1000 0.20004.1000 0.20006.1000 0.20007.1000 0.20002.1000 1.00003.1000 1.00004.1000 1.00006.1000 1.00007.1000 1.00002.1000 1.80003.1000 1.80004.1000 1.80006.1000 1.80007.1000 1.80002.1000 2.60003.1000 2.60004.1000 2.60006.1000 2.60007.1000 2.6000]%(x -0.5)*1.5*2 % define x to -1.5--0.5; [xn,xs]=mapminmax(x);xn=xn';%x1 = x(1,:);%x2 = x(2,:);f = [21.300022.200023.400021.300021.300017.500019.200025.400017.600017.500017.000023.700029.000022.100017.000025.100028.400029.200027.600025.1000]%10 + x1.^2.*cos(2*pi.*x1) + x2.^2.*cos(2*pi.*x2);f=f';%subplot(2,1,1)%plot(x,f)%hold on%subplot(2,1,2)plot(xn,f)title('initial data')% generate neural networknet =newrb(xn,f)%x1=1;%x2=1;%a=10 + x1.^2.*cos(2*pi.*x1) + x2.^2.*cos(2*pi.*x2)h=sim(net,xn)%subplot(2,1,1)plot(x,h)title('forecast function')%hold on、////////////////////////////////////////////////////////////////////////////////////////////////////////clearclcclose allx=0:20:500;y=[3200 3400 3600 3800 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015 4015]; f=polyfit(x,y,2);x0=0:20:500;y0=polyval(f,x0);plot(x,y,'*',x0,y0);xlabel('横坐标x')ylabel('纵坐标y')title('描点绘制曲线')legend('描点坐标','绘制的曲线')grid offf1=5;f2=100;T=10;f0=(f2+f1)/2;f3=f2-f1;t=0:0.1:10;s=sin(2*pi*(f0+(f3*t)/(2*T)).*t); plot(t,s)function K=f1(Nx,i,j,Yj,Y1,Y2,Y3) K1=0;K2=1.51;K3=0.93;K4=2.91;if j<Nx(i)&Nx(i)>0K=K1;elseif Yj(i)<Y1K=K2;endif Yj(i)<Y2&Yj(i)>Y1K=K3;endif Yj(i)>Y2K=K4;endEndfunction S=f2(Nx,i,j,Xw,Xj,Yw,Yj,dl,K,MN)R=0.01;h1=2074;Tf1=50.00000;if Nx(i)==jdt=abs(R-(Xj(j)^2+Yj(i)^2)^0.5);dV=(Xw(j)-Xw(j-1))*(Yw(i)-Yw(i-1));Sc=pi/2/MN*R*h1*K*Tf1/(K+h1*dt)/dV;Sp=-pi/2/MN*R*h1*K/(K+h1*dt)/dV;elseSc=0;Sp=0;endS=[Sc;Sp];function Heat4%边界不规则R=0.01; X=0.125; Y1=0.03;Y_1=-0.03;Y2=0.05;Y3=0.07;MN=25;N1=115;M1=20;M_1=20;M2=20;M3=20;K1=0;K2=1.51;K3=0.93;K4=2.91;h1=2309;Tf1=50.00000;h2=60;Tf2=18.00000; dl=pi/2/MN*R; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for i=1:MN+1Xw0(i)=R*sin((i-1)*pi/2/MN);endfor i=1:2*MN+1Yw0(i)=-R*cos((i-1)*pi/2/MN);endXw1=R:(X-R)/N1:X;Xw=[Xw0,Xw1(2:N1+1)];Yw_1=Y_1:(-R-Y_1)/M_1:-R;Yw1=R:(Y1-R)/M1:Y1;Yw2=Y1:(Y2-Y1)/M2:Y2;Yw3=Y2:(Y3-Y2)/M3:Y3;Yw=[Yw_1,Yw0(2:2*MN),Yw1,Yw2(2:M2+1),Yw3(2:M3+1)]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%M=length(Yw)-1;N=length(Xw)-1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Xj(1)=Xw(1);for j=2:N+1Xj(j)=(Xw(j)+Xw(j-1))/2;endXj(N+2)=Xw(N+1);Yj(1)=Yw(1);for i=2:M+1Yj(i)=(Yw(i)+Yw(i-1))/2;endYj(M+2)=Yw(M+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:M+2if i<=M_1+1Nx(i)=0;endif i>M_1+1&i<=M_1+MN+1Nx(i)=i-M_1;endif i>M_1+MN+1&i<M_1+2*MN+2Nx(i)=M_1+2*MN+3-i;endif i>M_1+2*MN+1Nx(i)=0;endend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=2:N+2iii=1;for i=1:M+2if Nx(i)==jNy(j,iii)=i;iii=iii+1;endendend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T=ones(M+2,N+2)*20.00000;for ii=1:3000iifor i=2:M+1for j=2:N+1Kp=f1(Nx,i,j,Yj,Y1,Y2,Y3);Kw=f1(Nx,i,j-1,Yj,Y1,Y2,Y3);Ke=f1(Nx,i,j+1,Yj,Y1,Y2,Y3);Ks=f1(Nx,i-1,j,Yj,Y1,Y2,Y3);Kn=f1(Nx,i+1,j,Yj,Y1,Y2,Y3);S=f2(Nx,i,j,Xw,Xj,Yw,Yj,dl,Kp,MN);Sc=S(1);Sp=S(2);if (Kp+Kw)==0Aw(j)=0;elseAw(j)=(Yw(i)-Yw(i-1))*Kw*Kp/(Kw*(Xj(j)-Xj(j-1))/2+Kp*(Xj(j)-Xj(j-1))/2) ;endif (Kp+Ke)==0Ae(j)=0;elseAe(j)=(Yw(i)-Yw(i-1))*Ke*Kp/(Ke*(Xj(j+1)-Xj(j))/2+Kp*(Xj(j+1)-Xj(j))/2) ;endif (Kp+Ks)==0As(j)=0;elseAs(j)=(Xw(j)-Xw(j-1))*Ks*Kp/(Ks*(Yj(i)-Yj(i-1))/2+Kp*(Yj(i)-Yj(i-1))/2);endif (Kp+Kn)==0An(j)=0;elseAn(j)=(Xw(j)-Xw(j-1))*Kn*Kp/(Kp*(Yj(i+1)-Yj(i))/2+Kn*(Yj(i+1)-Yj(i))/2);endAp(j)=-Aw(j)-Ae(j)-As(j)-An(j)+Sp*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1));B(j)=-Sc*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1))-As(j)*T(i-1,j)-An(j)*T(i+1,j);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%上下边界条件if i==2Ap(j)=-Aw(j)-Ae(j)-An(j)+Sp*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1));B(j)=-Sc*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1))-An(j)*T(i+1,j);endif i==M+1KK=K4;%2.7;Ap(j)=Ap(j)+An(j)*KK/((Yj(M+2)-Yj(M+1))*h2+KK);B(j)=-Sc*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1))-As(j)*T(i-1,j)-An(j)*Tf2*h2*(Yj(M+2)-Yj(M+1))/(( Yj(M+2)-Yj(M+1))*h2+KK);endend%%%%%%%%%%%%%%%%%%%%%%%%左右边界条件Ap(2)=Ap(2)+Aw(2);Ap(N+1)=Ap(N+1)+Ae(N+1);if Nx(i)==0A=diag(Ap(2:N+1))+diag(Aw(3:N+1),-1)+diag(Ae(2:N),1);T1=A\B(2:N+1)';T(i,2:N+1)=T1';elseA=diag(Ap(Nx(i):N+1))+diag(Aw(Nx(i)+1:N+1),-1)+diag(Ae(Nx(i):N),1);T2=A\B(Nx(i):N+1)';T(i,Nx(i):N+1)=T2';T(i,2:Nx(i)-1)=Tf1;%20;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=2:N+1for i=2:M+1Kp=f1(Nx,i,j,Yj,Y1,Y2,Y3);Kw=f1(Nx,i,j-1,Yj,Y1,Y2,Y3);Ke=f1(Nx,i,j+1,Yj,Y1,Y2,Y3);Ks=f1(Nx,i-1,j,Yj,Y1,Y2,Y3);Kn=f1(Nx,i+1,j,Yj,Y1,Y2,Y3);S=f2(Nx,i,j,Xw,Xj,Yw,Yj,dl,Kp,MN);Sc=S(1);Sp=S(2);if (Kp+Kw)==0Aw(i)=0;elseAw(i)=(Yw(i)-Yw(i-1))*Kw*Kp/(Kw*(Xj(j)-Xj(j-1))/2+Kp*(Xj(j)-Xj(j-1))/2) ;endif (Kp+Ke)==0Ae(i)=0;elseAe(i)=(Yw(i)-Yw(i-1))*Ke*Kp/(Ke*(Xj(j+1)-Xj(j))/2+Kp*(Xj(j+1)-Xj(j))/2) ;endif (Kp+Ks)==0As(i)=0;elseAs(i)=(Xw(j)-Xw(j-1))*Ks*Kp/(Ks*(Yj(i)-Yj(i-1))/2+Kp*(Yj(i)-Yj(i-1))/2);endif (Kp+Kn)==0An(i)=0;elseAn(i)=(Xw(j)-Xw(j-1))*Kn*Kp/(Kp*(Yj(i+1)-Yj(i))/2+Kn*(Yj(i+1)-Yj(i))/2);endAp(i)=-Aw(i)-Ae(i)-As(i)-An(i)+Sp*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1));B(i)=-Sc*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1))-Aw(i)*T(i,j-1)-Ae(i)*T(i,j+1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%左右边界条件Ap(i)=Ap(i)+Aw(i);%-Ae(i)-As(i)-An(i)+Sp*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1));B(i)=-Sc*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1))-Ae(i)*T(i,j+1);endif j==N+1Ap(i)=Ap(i)+Ae(i);%-Aw(i)-As(i)-An(i)+Sp*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1));B(i)=-Sc*(Yw(i)-Yw(i-1))*(Xw(j)-Xw(j-1))-Aw(i)*T(i,j-1);endend%%%%%%%%%%%%%%%%%%%%%%%%上下边界条件Ap(2)=Ap(2)+As(2);KK=K4;%2.7;Ap(M+1)=Ap(M+1)+An(M+1)*KK/((Yj(M+2)-Yj(M+1))*h2+KK);B(M+1)=B(M+1)-An(M+1)*Tf2*h2*(Yj(M+2)-Yj(M+1))/((Yj(M+2)-Yj(M+1))*h2+KK);if j<=length(Ny)A1=diag(Ap(2:Ny(j,1)))+diag(As(3:Ny(j,1)),-1)+diag(An(2:Ny(j,1)-1),1);T3=A1\B(2:Ny(j,1))';A2=diag(Ap(Ny(j,2):M+1))+diag(As(Ny(j,2)+1:M+1),-1)+diag(An(Ny(j,2):M),1);T4=A2\B(Ny(j,2):M+1)';T(2:Ny(j,1),j)=T3;T(Ny(j,2):M+1,j)=T4;elseA3=diag(Ap(2:M+1))+diag(As(3:M+1),-1)+diag(An(2:M),1);T5=A3\B(2:M+1)';T(2:M+1,j)=T5;endend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%A=0for i=2:N+1A=A+T(M+1,i)endt=A/Nend%迭代循环%K=2.7%for i=1:M+2% S=f2(Nx,i,j,Xw,Xj,Yw,Yj,dl,K,MN);% KK(i,j)=S(1);% KK1(i,j)=S(2);% end%end%KK1surf(Xj(2:N+1),Yj(2:M+1),T(2:M+1,2:N+1))shading interpview(0, 90);figure,contour(Xj(2:N+1),Yj(2:M+1),T(2:M+1,2:N+1));mydir=uigetdir('选择一个目录');if mydir(end)~='\'mydir=[mydir,'\'];end;files=dir([mydir,'*.bmp']);y=length(files);mkdir('C:\Documents and Settings\Administrator\桌面\','huidutuzhuose'); for x=1:yI=imread([mydir,files(x).name]);I1=rgb2gray(I)B=I1;[m,n]=size(I1);[L,num]=bwlabel(I1);stats=regionprops(L,'Area');allArea=[stats.Area];M=max(allArea)allArea=allArea/(M/60);B=double(B);%set('B',[1 1 1]);for l=1:numfor a=1:mfor b=1:nif L(a,b)==lB(a,b)=B(a,b)-1+allArea(l);else B(a,b)=B(a,b);endendendendfigure,imagesc(B)daspect([1 1 1]);%将图片还原为原来大小colorbarF=getframe(gcf);imwrite(F.cdata,['C:\Documents and Settings\Administrator\桌面\huidutuzhuose\',files(x).name]); %close all;endfunction [ f ] = dfft( x,n )%dfft :计算序列的n点fft并绘制图形；% x ：输入序列向量；% n ：序列fft的点数；% f : x的fft；f=fft(x,n);subplot(2,1,1)stem(0:length(x)-1,x,'.')title('信号原型')axis([0 length(x)-1 min(x) max(x)]) subplot(2,1,2)stem(0:n-1,abs(f),'.')title('上图序列FFT后的图象')end调用：a=1:5dfft(a,5)自组织竞争神经网络：load abs.txt;data = abs;data = data(:,2:115);p1 = data(1:40,:);t1 = data(41:60,:);p1 = p1';t1 = t1';% Data Normalization[pn,ps] = mapminmax(p1); [tn,ts] = mapminmax(t1);Q = minmax(pn);% Generate Neural Network net = newc(Q,2,0.01)net = init(net)net.trainParam.epochs = 300;% Trainingnet = train(net,pn)% Textinga = sim(net,pn)建模2012：%% 通道3温度和风速clcclearclose allx=[0.3 0.9 1.5 2.1 2.7];y=[2.4 5 7.2];Z=[27 29 29 30 2930 29 31 32 3027 31 31 52 31];ZN=[0.4 0.6 0.7 0.8 0.90.4 0.5 0.6 0.7 0.60.4 0.6 0.6 0.6 0.5];%subplot(2,1,1);mesh(x,y,z);title('通道3温度分布'); %% 温度分布xi=0.3:0.01:2.7;yi=2:0.01:7.5;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=interp2(x,y,Z,XI,YI,'cubic');figuremesh(XI,YI,ZI);title('热通道3温度分布');xlabel('高度z')ylabel('离空调的距离y')%% 风速分布xi=0.3:0.01:2.7;yi=2:...0.01:7.5;[XI,YI]=meshgrid(xi,yi);ZI=griddata(x,y,ZN,XI,YI,'cubic');figuremesh(XI,YI,ZI);title('热通道3风速分布');xlabel('高度z')ylabel('离空调的距离y')clcclearclose alla=[1.1 2.1 0.2 21.3 21.3 21.3 21.3 1.1 3.1 0.2 22.2 22.2 22.2 22.21.1 4.1 0.2 23.4 23.6 23.4 23.41.1 6.1 0.2 21.3 21.3 21.3 21.31.1 7.1 0.2 21.3 21.3 21.3 21.31.12.1 1 17.5 17.5 17.5 17.51.1 3.1 1 19.2 19.2 19.2 19.21.1 4.1 1 25.4 26.4 25.4 25.41.1 6.1 1 17.6 17.6 17.6 17.61.1 7.1 1 17.5 17.5 17.5 17.51.12.1 1.8 17.0 17.0 17.0 17.01.1 3.1 1.8 23.7 23.8 24.3 23.71.1 4.1 1.8 29.0 31.5 29.0 29.01.1 6.1 1.8 22.1 22.1 22.1 22.11.1 7.1 1.8 17.0 17.0 17.0 17.01.12.1 2.6 25.1 25.1 25.1 25.11.1 3.12.6 28.4 28.9 29.6 28.41.1 4.12.6 29.2 31.1 29.7 29.21.1 6.12.6 27.6 27.6 27.7 27.61.1 7.12.6 25.1 25.1 25.1 25.12.4 2.1 0.2 16.0 16.0 16.0 16.02.43.1 0.2 16.0 16.0 16.0 16.02.4 4.1 0.2 25.3 25.9 25.3 25.32.4 6.1 0.2 16.0 16.0 16.0 16.02.4 7.1 0.2 16.0 16.0 16.0 16.02.4 2.1 1 16.0 16.0 16.0 16.02.43.1 1 16.0 16.0 16.0 16.02.4 4.1 1 28.7 31.0 28.7 28.72.4 6.1 1 16.0 16.0 16.0 16.02.4 7.1 1 16.0 16.0 16.0 16.02.4 2.1 1.8 16.5 16.5 16.5 16.52.43.1 1.8 20.4 20.4 20.7 20.42.4 4.1 1.8 32.1 35.9 32.4 32.22.4 6.1 1.8 20.4 20.4 20.4 20.42.4 7.1 1.8 16.5 16.5 16.5 16.52.4 2.1 2.6 24.7 24.8 24.8 24.82.43.1 2.6 27.5 27.6 28.5 27.52.4 4.1 2.6 29.3 30.7 30.3 29.32.4 6.1 2.6 27.5 27.5 27.7 27.52.4 7.1 2.6 24.7 24.8 24.8 24.84.1 3.1 0.2 16.0 16.0 16.0 16.0 4.1 4.1 0.2 25.3 25.3 25.9 25.3 4.16.1 0.2 16.0 16.0 16.0 16.0 4.17.1 0.2 16.0 16.0 16.0 16.0 4.1 2.1 1 16.0 16.0 16.0 16.0 4.1 3.1 1 16.0 16.0 16.0 16.0 4.1 4.1 1 28.7 28.7 31.0 28.7 4.1 6.1 1 16.0 16.0 16.0 16.0 4.1 7.1 1 16.0 16.0 16.0 16.0 4.1 2.1 1.8 16.4 16.4 16.4 16.4 4.1 3.1 1.8 19.4 19.3 19.6 19.4 4.1 4.1 1.8 32.0 32.0 35.9 32.0 4.1 6.1 1.8 19.4 19.3 19.4 19.3 4.1 7.1 1.8 16.4 16.4 16.4 16.4 4.1 2.1 2.6 24.6 24.6 24.7 24.6 4.1 3.1 2.6 27.4 27.4 28.4 27.5 4.1 4.1 2.6 29.2 29.2 32.5 29.3 4.1 6.1 2.6 27.4 27.4 27.5 27.44.1 7.1 2.6 24.6 24.6 24.6 24.65.1 2.1 0.2 16.0 16.0 16.0 16.0 5.1 3.1 0.2 16.0 16.0 16.0 16.0 5.1 4.1 0.2 25.3 25.3 25.3 25.9 5.1 6.1 0.2 16.0 16.0 16.0 16.0 5.17.1 0.2 16.0 16.0 16.0 16.0 5.1 2.1 1 16.0 16.0 16.0 16.0 5.1 3.1 1 16.0 16.0 16.0 16.0 5.1 4.1 1 28.7 28.7 28.7 32.5 5.1 6.1 1 16.0 16.0 16.0 16.0 5.1 7.1 1 16.0 16.0 16.0 16.0 5.1 2.1 1.8 16.1 16.1 16.1 16.1 5.1 3.1 1.8 16.5 16.5 16.6 16.6 5.1 4.1 1.8 31.2 31.2 31.2 34.7 5.1 6.1 1.8 16.5 16.5 16.5 16.5 5.1 7.1 1.8 16.1 16.1 16.1 16.1 5.1 2.1 2.6 24.6 24.6 24.7 24.7 5.1 3.1 2.6 27.4 27.4 27.9 28.1 5.1 4.1 2.6 29.2 29.2 29.5 31.4 5.1 6.1 2.6 27.4 27.4 27.5 27.5 5.1 7.1 2.6 24.6 24.6 24.6 24.6 7.2 2.1 0.2 16.0 16.0 16.0 16.0 7.2 3.1 0.2 16.0 16.0 16.0 16.0 7.2 4.1 0.2 25.3 25.3 25.3 25.3 7.2 6.1 0.2 16.0 16.0 16.0 16.07.2 2.1 1 16.0 16.0 16.0 16.07.2 3.1 1 16.0 16.0 16.0 16.07.2 4.1 1 28.7 28.7 28.7 28.77.2 6.1 1 16.0 16.0 16.0 16.07.2 7.1 1 16.0 16.0 16.0 16.07.2 2.1 1.8 16.5 16.5 16.5 16.57.2 3.1 1.8 20.7 20.6 20.7 20.97.2 4.1 1.8 32.2 32.2 32.2 32.67.2 6.1 1.8 20.7 20.6 20.6 20.67.2 7.1 1.8 16.5 16.5 16.5 16.57.2 2.1 2.6 24.8 24.8 24.8 24.87.2 3.1 2.6 27.5 27.5 27.7 28.47.2 4.1 2.6 29.3 29.3 29.4 30.57.2 6.1 2.6 27.5 27.5 27.6 27.67.2 7.1 2.6 24.8 24.8 24.8 24.88 2.1 0.2 21.1 21.1 21.1 21.18 3.1 0.2 21.2 21.2 21.2 21.28 4.1 0.2 23.5 23.5 23.5 23.58 6.1 0.2 21.2 21.2 21.2 21.28 7.1 0.2 21.1 21.1 21.1 21.18 2.1 1 17.2 17.2 17.2 17.28 3.1 1 17.3 17.3 17.3 17.38 4.1 1 25.7 25.7 25.7 25.78 6.1 1 17.3 17.3 17.3 17.38 7.1 1 17.2 17.2 17.2 17.28 2.1 1.8 17.1 17.1 17.1 17.18 3.1 1.8 22.6 22.5 22.7 23.08 4.1 1.8 29.5 29.5 29.5 29.58 6.1 1.8 22.6 22.5 22.6 22.68 7.1 1.8 17.1 17.1 17.1 17.18 2.1 2.6 25.1 25.2 25.2 25.28 3.1 2.6 27.6 27.6 27.8 28.58 4.1 2.6 29.2 29.2 29.3 29.98 6.1 2.6 27.6 27.6 27.6 27.78 7.1 2.6 25.1 25.2 25.2 25.2]y=a(:,1);z=a(:,2);x=a(:,3);scatter3(z,y,x,'DisplayName','z,y,x');figure(gcf)%探测点的分布grid offxlabel('y')ylabel('x')zlabel('z')%hold on%boxsurface([0.5 1 1],[1 2 2])hold onboxplot3(0,0,0, 10.5,9.2,3.2)%% 第一个机柜群hold onboxplot3(2.5,1.2,0, 6.4,0.8,2)%% 第二个机柜群hold onboxplot3(2.5,3.2,0, 6.4,0.8,2) %% 第三个机柜群hold onboxplot3(2.5,5.2,0, 6.4,0.8,2) %% 第四个机柜群hold onboxplot3(2.5,7.2,0, 6.4,0.8,2)%% 空调1hold onboxplot3(0,1.7,0, 0.9,1.8,2)%% 空调2hold onboxplot3(0,5.7,0, 0.9,1.8,2)%% 通风口1hold onboxplot3(0.2,1.9,2, 0.5,1.4,.01) %% 通风口1hold onboxplot3(0.2,5.9,2, 0.5,1.4,.01)if 0%% 所有机柜均为0.5时，高度为0.2的温度分布figurem=a(1:5,4);n=a(21:25,4);o=a(41:45,4);p=a(61:65,4);q=a(81:85,4);r=a(101:105,4);x=[1.1 2.4 4.1 5.1 7.2 8];y=[2.1 3.1 4.1 6.1 7.1];Z=[m,n,o,p,q,r];xi=1:0.01:9;yi=0:0.02:9;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(x,y,Z,XI,YI,'v4');mesh(XI,YI,ZI);title('所有机柜工作量为0.5，高度为0.2的温度分布');%% 所有机柜均为0.5时，高度为1米处的温度分布figurex=[1.1 2.4 4.1 5.1 7.2 8];y=[2.1 3.1 4.1 6.1 7.1];Z=[m,n,o,p,q,r];m=a(6:10,4);n=a(26:30,4);o=a(46:50,4);p=a(66:70,4);q=a(86:90,4);r=a(106:110,4);%subplot(2,1,1);mesh(x,y,z);title('通道3温度分布');xi=1:0.01:9;yi=0:0.02:9;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(x,y,Z,XI,YI,'v4');mesh(XI,YI,ZI);title('所有机柜工作量为0.5，高度为1米处的温度分布');%% 所有机柜均为0.5时，高度为1.8米处的温度分布figurex=[1.1 2.4 4.1 5.1 7.2 8];y=[2.1 3.1 4.1 6.1 7.1];Z=[m,n,o,p,q,r];m=a(11:15,4);n=a(31:35,4);o=a(51:55,4);p=a(71:75,4);q=a(91:95,4);r=a(111:115,4);xi=1:0.01:9;yi=0:0.02:9;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(x,y,Z,XI,YI,'v4');mesh(XI,YI,ZI);title('所有机柜工作量为0.5，高度为1.8米处的温度分布'); %% 所有机柜均为0.5时，高度为2.6米处的温度分布figurex=[1.1 2.4 4.1 5.1 7.2 8];y=[2.1 3.1 4.1 6.1 7.1];Z=[m,n,o,p,q,r];m=a(16:20,4);n=a(36:40,4);o=a(56:60,4);p=a(76:80,4);q=a(96:100,4);r=a(116:120,4);%subplot(2,1,1);mesh(x,y,z);title('通道3温度分布');xi=1:0.01:9;yi=0:0.02:9;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(x,y,Z,XI,YI,'v4');mesh(XI,YI,ZI);title('所有机柜工作量为0.5，高度为2.6米处的温度分布');%% 所有机柜均为0.5时,通道5垂直截面温度分布figures=a(4:8,4);t=a(9:13,4);u=a(14:18,4);v=a(19:23,4);Z=[s,t,u,v]x=[2.1 3.1 4.1 6.1 7.1];y=[0.2 1 1.8 2.6];% subplot(2,1,1);mesh(x,y,z);title('所有机柜均为0.5时，通道4垂直截面温度分布');xi=1:0.001:3;yi=2:0.02:7.5;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(y,x,Z,XI,YI,'cubic');mesh(XI,YI,ZI);title('所有机柜工作量均为0.5时，通道5垂直截面温度分布');%% 所有机柜均为0.5时, 通道3右垂直截面温度分布figures=a(41:45,4);t=a(46:50,4);u=a(51:55,4);v=a(56:60,4);Z=[s,t,u,v]x=[2.1 3.1 4.1 6.1 7.1];y=[0.2 1 1.8 2.6];% subplot(2,1,1);mesh(x,y,z);title('所有机柜均为0.5时，通道4垂直截面温度分布');xi=1:0.001:3;yi=2:0.02:7.5;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(y,x,Z,XI,YI,'cubic');mesh(XI,YI,ZI);title('所有机柜工作量均为0.5时，通道3右垂直截面温度分布'); %% 所有机柜均为0.5时, 通道3左垂直截面温度分布figures=a(61:65,4);t=a(66:70,4);u=a(71:75,4);v=a(76:80,4);Z=[s,t,u,v]x=[2.1 3.1 4.1 6.1 7.1];y=[0.2 1 1.8 2.6];% subplot(2,1,1);mesh(x,y,z);title('所有机柜均为0.5时，通道4垂直截面温度分布');xi=1:0.001:3;yi=2:0.02:7.5;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(y,x,Z,XI,YI,'cubic');mesh(XI,YI,ZI);title('所有机柜工作量均为0.5时，通道3左垂直截面温度分布'); %% 所有机柜均为0.5时, 通道2垂直截面温度分布figures=a(81:85,4);t=a(86:90,4);u=a(91:95,4);v=a(96:100,4);Z=[s,t,u,v]x=[2.1 3.1 4.1 6.1 7.1];y=[0.2 1 1.8 2.6];% subplot(2,1,1);mesh(x,y,z);title('所有机柜均为0.5时，通道4垂直截面温度分布');xi=1:0.001:3;yi=2:0.02:7.5;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(y,x,Z,XI,YI,'cubic');mesh(XI,YI,ZI);title('所有机柜工作量均为0.5时，通道2垂直截面温度分布'); %% 所有机柜均为0.5时, 通道1垂直截面温度分布figures=a(101:105,4);t=a(106:110,4);u=a(111:115,4);v=a(116:120,4);Z=[s,t,u,v]x=[2.1 3.1 4.1 6.1 7.1];y=[0.2 1 1.8 2.6];% subplot(2,1,1);mesh(x,y,z);title('所有机柜均为0.5时，通道4垂直截面温度分布');xi=1:0.001:3;yi=2:0.02:7.5;[XI,YI]=meshgrid(xi,yi);%ZI=INTERP2(x,y,Z,XI,YI,'cubic');ZI=griddata(y,x,Z,XI,YI,'cubic');mesh(XI,YI,ZI);title('所有机柜工作量均为0.5时，通道1垂直截面温度分布'); endclcclearclf[X,Y]=meshgrid(-pi/2:0.1:pi/2,-pi:0.2:pi);Z=abs(sin(Y)*cos(X));[DX,DY]=gradient(Z,0.1,0.2);%contour(Z);hold on;quiver(DX,DY);hold off;function boxplot3(x0,y0,z0,Lx,Ly,Lz)%(x0,y0,z0)是第一个顶点的位置; (Lx,Ly,Lz)是长方体的长宽高.x=[x0 x0 x0 x0 x0+Lx x0+Lx x0+Lx x0+Lx];y=[y0 y0 y0+Ly y0+Ly y0 y0 y0+Ly y0+Ly];z=[z0 z0+Lz z0+Lz z0 z0 z0+Lz z0+Lz z0];index=zeros(6,5);index(1,:)=[1 2 3 4 1];index(2,:)=[5 6 7 8 5];index(3,:)=[1 2 6 5 1];index(4,:)=[4 3 7 8 4];index(5,:)=[2 6 7 3 2];index(6,:)=[1 5 8 4 1];for k=1:6plot3(x(index(k,:)),y(index(k,:)),z(index(k,:)))hold onendboxplot3(0,0,0,5,5,5)function boxsurface(p0,l)[x,y,z]=planarsurface(p0,p0+[0 0 l(3)],p0+[0 l(2) 0]);surf(x,y,z)hold on[x,y,z]=planarsurface(p0+[l(1) 0 0],p0+[l(1) 0 l(3)],p0+[l(1) l(2) 0]);surf(x,y,z)[x,y,z]=planarsurface(p0,p0+[0 0 l(3)],p0+[l(1) 0 0]);surf(x,y,z)[x,y,z]=planarsurface(p0+[0 l(2) 0],p0+[0 l(2) l(3)],p0+[l(1) l(2) 0]);surf(x,y,z)[x,y,z]=planarsurface(p0,p0+[l(1) 0 0],p0+[0 l(2) 0]);surf(x,y,z)[x,y,z]=planarsurface(p0+[0 0 l(3)],p0+[l(1) 0 l(3)],p0+[0 l(2) l(3)]);surf(x,y,z)axis equalaxis offendfunction [xx,yy,zz,l]=planarsurface(p0,p1,p2)v=p1-p0;w=p2-p0;s=0:0.2:1;l=length(s);[s,t]=meshgrid(s,s);xx=p0(1)+s*v(1)+t*w(1);yy=p0(2)+s*v(2)+t*w(2);zz=p0(3)+s*v(3)+t*w(3);调用：p0=[0 0 0];I=[5 5 5];boxsurface(p0,I)clear allclcclose allim=[];%%for i1=1:26;for i2=1:26;out=strcat(num2str(i1),'x',num2str(i2),'y.txt');%文件名 I=load(out);out1(i1,i2)=max(abs(I(:,2)));endend%%im=reshape(out1,26,26);imagesc(im);光学实验：药物图象。