日志
matlab source
|
function [x,n]=jacobi(A,b,x0,eps,varargin)
if nargin==3, eps= 1.0e-6; M = 200;
elseif nargin<3, error; return
elseif nargin ==5, M = varargin{1};
endD=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
B=D\(L+U); f=D\b; x=B*x0+f;
n=1; %迭代次数while norm(x-x0)>=eps
x0=x; x=B*x0+f; n=n+1;
if(n>=M)
disp('Warning: 迭代次数太多，可能不收敛！'); return;
end
end
function FilesAll(iniPath)
%iniPath='D:\My Documents\MATLAB';%初始化路径
if (iniPath(end)~='\'), iniPath=[iniPath,'\']; end
files=dir(fullfile(iniPath,'*.*'));%过滤
n=size(files,1);%文件个数
for i=1:n
tName=files(i,1).name;
if (length(tName)>2) %排除'.'和'..'
if (~files(i,1).isdir) %非目录
filename=fullfile(iniPath,tName);%文件名及路径
disp(filename);
else %是目录
iniPaths=[iniPath,tName]; FilesAll(iniPaths);
end
end
end
end
% 一维极值点
a = cumsum(randn(1000,1)); peaks = localMaximum(a,[100]);
figure; plot(a); hold on; plot(peaks,a(peaks),'ro');%二维矩阵极值点
[x y] = meshgrid(-4:0.1:4,-4:0.1:4); a = sinc(x)+sinc(y);
lMaxInd = localmaximum(a,[5 5]); %极大值点
lMinInd = localmaximum(-a,[5 5]);%极小值点
surf(x,y,a); shading interp; hold on;
plot3(x(lMaxInd),y(lMaxInd),a(lMaxInd),'k*','markersize',10,'linewidth',1.5);
plot3(x(lMinInd),y(lMinInd),a(lMinInd),'g*','markersize',10','linewidth',1.5);
% legend('sinc(x)+sinc(y)','peaks','valleys','location','best')localmaximum.m函数：http://www.mathworks.com/matlabcentral/fileexchange/14498-local-maxima-minima
%方法一
x=0:0.01:20; y=2*sin(x/2)+cos(2*x)/2;
% indmax=find(diff(sign(diff(y)))<0)+1;%极大值点
% indmin=find(diff(sign(diff(y)))>0)+1;%极小值点
ind=find(diff(sign(diff(y)))~=0)+1;%极值点
plot(x,y,x(ind),y(ind),'ro')%方法二
x=0:0.01:20; y=2*sin(x/2)+cos(2*x)/2;
indmax=intersect(find(diff(y)>0)+1,find(diff(y)<0));%极大值点
indmin=intersect(find(diff(y)<0)+1,find(diff(y)>0));%极小值点
plot(x,y,x(indmax),y(indmax),'ro',x(indmin),y(indmin),'go')%方法三
x=0:0.01:20; y=2*sin(x/2)+cos(2*x)/2;
yf=y(2:end-1)-y(1:end-2); %前向差分
yb=y(2:end-1)-y(3:end); %后向差分
I=find(yf.*yb>=0)+1; %寻找极值点(不考虑首尾两点)
plot(x,y,'b',x(I),y(I),'r*');
【转】matlab求图像的局部峰值点集 2010-10-21 23:28
经过辛苦寻找，目前发现三种方法来做（其中alpha控制强度）：
I = imread('lena-gray.bmp'); BW = nlfilter(I, [3 3], @(x) all(x(5) > x([1:4 6:9])+ alpha ) ); imshow(BW); %方法1
I = imread('lena-gray.bmp'); BW = imregionalmax(I); imshow(BW); %方法2CC = bwconncomp(BW);
for i = 1:CC.NumObjects,
index = CC.PixelIdxList{i};
if (numel(index) > 1), BW(index) = false; end
end
I = imread('lena-gray.bmp'); BW = I > imdilate(I, [1 1 1; 1 0 1; 1 1 1]) +alpha ; imshow(BW); %方法3
自定义二维隐函数绘图
Matlab中提供的ez开头的函数，基本都支持二维隐函数绘图如果细心的网友，肯定注意到了，其实所谓的隐函数绘图，其实器本质是使用了一个contour()函数下面是我们自己使用contour函数改造的二维隐函数绘图函数function implot(fun,rangexy,ngrid)
%二维隐函数绘图
%
%输入参数说明
% -fun 函数句柄，可以是匿名、inline和M函数
% -rangexy=[xmin xmax ymin ymax] 绘图范围，默认[-2*pi 2*pi]
% -ngrid 绘图时计算的点数，初值是20，然后逐步加细，默认50
%
% Example
% 绘制y^3+exp(y)-tanh(x)=0的图形
%>> f=inline('y^3+exp(y)-tanh(x)','x','y')
% >>implot(f,[-3 3 -2 1])
%
%原理说明：
%其实该函数就是调用了contour()函数，绘制隐函数在xoy平面上的等高线，就得到了二维隐函数的图像
%
%by dynamic
%all rights reserved by%if nargin == 1 ;% 使用默认的rangxy和ngrid
rangexy=[-2*pi,2*pi,-2*pi,2*pi];
ngrid=50;
endif nargin == 2; % 使用默认ngrid
ngrid=50;
end
%生成2D网格
xm=linspace(rangexy(1),rangexy(2),ngrid);
ym=linspace(rangexy(3),rangexy(4),ngrid);
[x,y]=meshgrid(xm,ym);
fvector=vectorize(fun);% 将目标函数矢量化
fvalues=feval(fvector,x,y); %计算函数值
%fvalues=fvector(x,y); % can also calculate directly from the vectorized inline function
contour(x,y,fvalues,[0,0],'b-');% 绘制z=0的等高线，即隐函数的图形
xlabel('x');ylabel('y');
grid
end
自定义三维隐函数绘制函数全世界人都知道Matlab那几千个库函数中，可恨的是，却偏偏不提供个三维隐函数的图像绘制的函数，真是郁闷死了
好，下面我们看一个三维隐函数绘制的应用，它灵活的应用了Matlab的isosurface()(等值面函数)，绘制出来的图效果还比较好function implot3(fun,rangexyz,ngrid,varargin)
%三维隐函数绘图
%
%输入参数说明
% -fun 函数句柄，可以是匿名、inline和M函数
% -rangexy=[xmin xmax ymin ymax,zmin,zmax] 绘图范围
% -ngrid 绘图时计算的点数
%
%Example
%fun=@(x,y,z)(x+y+z).*(x.*y+x.*z+y.*z)-10*x.*y.*z;
%rangexyz=[1 10 1 10 1 10];ngrid=50;
%implot3(fun,rangexyz,ngrid)
%
%rewrite by dynamic
%all rights reserved byx=linspace(rangexyz(1),rangexyz(2),ngrid);
y=linspace(rangexyz(3),rangexyz(4),ngrid);
z=linspace(rangexyz(5),rangexyz(6),ngrid);
[xx,yy,zz]=meshgrid(x,y,z);
%fvector=vectorize(fun);% 将目标函数矢量化
f=feval(fun,xx,yy,zz,varargin{:});
p=patch(isosurface(xx,yy,zz,f,0),varargin{:});
set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1]) ; view(3) ; camlight; lighting phong
end
【转】Matlab画三维隐函数曲面 2010-10-13 17:34
1.隐函数曲面网格图函数implicitmesh
function h=implicitmesh(f,xlimit,ylimit,zlimit,gd)
%implicitmesh(f,span,gd):画隐函数曲面f(x,y,z)=0的网格图,
% 各坐标范围均限定在span=[lb,ub], 网格数为gd，默认为25
%implicitmesh(f,xspan,yspan,zspan,gd):画隐函数曲面f(x,y,z)=0, 各坐标范围分别限定在xspan,yspan,zspan
%h=implicitmesh(...)：画隐函数曲面并输出句柄
%例一：implicitmesh(inline('x.*y+z.^2'),[-5 5])%注意*\^一定要设成点运算
%例二：f=@(x,y,z)x.^2+y.^2+0*z-1;%注意如果f中不含某个变量一定要加上诸如0*y的项。
%implicitmesh(f,[-1 1],10)
%例三：
%f=@(x,y,z)(x.^2 + (9/4)*y.^2 + z.^2 - 1).^3 - x.^2.*z.^3 - (9/80)*y.^2.*z.^3;
%g=@(x,y,z)(sqrt(x.^2+y.^2)-2).^2+z.^2-.09;
%implicitmesh(f,[-1.5 1.5],[-.8 .8],[-1.5 1.5],50);
%hold on%可以添加图形
%h=implicitmesh(g,[-2.3,2.3]);
%colormap hsv;set(h,'facecolor','none');%可以设置各种效果
%axis off;axis equal;
if nargin==2, ylimit=xlimit;zlimit=xlimit;gd=25;
elseif nargin==3, gd=ylimit;ylimit=xlimit;zlimit=xlimit;
elseif nargin==4, gd=25;
elseif nargin==5
else error('Error in input arguments')
end
x=linspace(xlimit(1),xlimit(2),gd); y=linspace(ylimit(1),ylimit(2),gd); z=linspace(zlimit(1),zlimit(2),gd);
[x,y,z]=meshgrid(x,y,z);val=f(x,y,z); [f,v]=isosurface(x,y,z,val,0);
if isempty(f), warning('There is no graph in the range.'); p=[];
else newplot; p=patch('Faces',f,'Vertices',v,'CData',v(:,3),'facecolor','w','EdgeColor','flat');
isonormals(x,y,z,val,p);view(3);grid on
end
if nargout==0, else h=p; end2.隐函数曲面表面图函数implicitsurf
function h=implicitsurf(f,xlimit,ylimit,zlimit,gd)
%implicitsurf(f,span,gd):画隐函数曲面f(x,y,z)=0的网格图,
% 各坐标范围均限定在span=[lb,ub], 网格数为gd，默认为25
%implicitsurf(f,xspan,yspan,zspan,gd):画隐函数曲面f(x,y,z)=0, 各坐标范围分别限定在xspan,yspan,zspan
%h=implicitsurf(...)：画隐函数曲面并输出句柄
%例一： implicitsurf(inline('x.*y+z.^2'),[-5 5])%注意*\^一定要设成点运算
%例二：f=@(x,y,z)x.^2+y.^2+0*z-1;%注意如果f中不含某个变量一定要加上诸如0*y的项。
%implicitsurf(f,[-1 1],10)
%例三：
%f=@(x,y,z)(x.^2 + (9/4)*y.^2 + z.^2 - 1).^3 - x.^2.*z.^3 - (9/80)*y.^2.*z.^3;
%g=@(x,y,z)(sqrt(x.^2+y.^2)-2).^2+z.^2-.09;
%h=implicitsurf(f,[-1.5 1.5],[-.8 .8],[-1.5 1.5],50);
%set(h,'AmbientStrength',.5);%可以设置各种效果
%hold on%可以添加图形
%h=implicitsurf(g,[-2.3,2.3],[-2.3,2.3],[-.3,.3]);
%colormap hsv;set(h,'AmbientStrength',.8,'FaceAlpha',.5);%可以设置各种效果
%axis off;axis equal;shading interp;camlight;lighting gouraud;
if nargin==2, ylimit=xlimit;zlimit=xlimit;gd=25;
elseif nargin==3, gd=ylimit;ylimit=xlimit;zlimit=xlimit;
elseif nargin==4, gd=25;
elseif nargin==5
else error('Error in input arguments')
end
x=linspace(xlimit(1),xlimit(2),gd); y=linspace(ylimit(1),ylimit(2),gd); z=linspace(zlimit(1),zlimit(2),gd);
[x,y,z]=meshgrid(x,y,z);val=f(x,y,z); [f,v]=isosurface(x,y,z,val,0);
if isempty(f), warning('There is no graph in the range.'); p=[];
else newplot;
p=patch('Faces',f,'Vertices',v,'CData',v(:,3),'facecolor','flat','EdgeColor','k');
isonormals(x,y,z,val,p);view(3);grid on
end
if nargout==0, else h=p; end
MATLAB 喊话，朗读句子 2010-09-23 20:45
function wav = matlabspeech(txt,voice,pace,fs)
% matlabspeech('I Love Maple and Matlab')
% matlabspeech('I Love Maple')
% matlabspeech('I Love Matlab')if ~ispc, error('要求使用Microsoft Win32 SAPI .'); end
if ~ischar(txt), error('第一个参数是字符串.'); end
SV = actxserver('SAPI.SpVoice'); TK = invoke(SV,'GetVoices');
if nargin > 1
for k = 0:TK.Count-1
if strcmpi(voice,TK.Item(k).GetDescription), SV.Voice = TK.Item(k); break;
elseif strcmpi(voice,'list'), disp(TK.Item(k).GetDescription); end
end
if nargin > 2, if isempty(pace), pace = 0; end
if abs(pace) > 10, pace = sign(pace)*10; end
SV.Rate = pace;
end
end
if nargin < 4 || ~ismember(fs,[8000,11025,12000,16000,22050,24000,32000,44100,48000]), fs = 16000; end
if nargout > 0
MS = actxserver('SAPI.SpMemoryStream');
MS.Format.Type = sprintf('SAFT%dkHz16BitMono',fix(fs/1000));
SV.AudioOutputStream = MS;
end
invoke(SV,'Speak',txt);
if nargout > 0
wav = reshape(double(invoke(MS,'GetData')),2,[])';
wav = (wav(:,2)*256+wav(:,1))/32768;
wav(wav >= 1) = wav(wav >= 1)-2;
delete(MS); clear MS;
end
delete(SV); clear SV TK; pause(0.2);
end
ha=axes('units','normalized');
% ha=axes('units','normalized','position',[0 0 1 1]);
uistack(ha,'down')
II=imread('xn3.bmp');%注意路径问题
image(II); colormap gray
set(ha,'handlevisibility','off','visible','off');
x=-pi:0.1:pi; y=x.*sin(x.*cos(x)).*tan(x); plot(x,y,'LineWidth',2)
set(gca,'color','none') %这里以前自己没有注意
title('设置axes背景为图片')
global rootBessel
maxv = 10; % 最大贝塞尔函数阶数-1，表示从0到9阶
maxs = 10; % 所需要计算贝塞尔函数的零点的数目
j = zeros(maxv, maxs); % 贝塞尔函数的根
incr = 4.0;
for v=0:maxv-1
h = v+1.9*v^(1/3)+1;
if (v==0) % 0阶贝塞尔函数的第一个零点
j(v+1,1) = fzero(@(x)besselj(v,x),2);
else % 1阶及以上阶贝塞尔函数的第一个零点
j(v+1,1) = fzero(@(x)besselj(v,x),h);
end
for s=2:maxs % 贝塞尔函数的第2个及后面的零点
j(v+1,s) = fzero(@(x)besselj(v,x),j(v+1,s-1)+incr);
end
end
rootBessel = j
%产生拟合曲线，并求某点导数
% hObject handle to btn_ployder (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
x=0.01:0.01:14;
steps=str2num(get(handles.edit_steps,'string')) ;
if length(steps)==1%必须是单个曲线
switch get(handles.popupmenu_pick,'value')
case 1, y=BESSELJ(steps(1),x);
case 2, y=BESSELY(steps(1),x);
end
a=polyfit(x,y,str2num(get(handles.edit_polyStep,'string')));%获取拟合多项式系数
a1=polyder(a);%多项式一阶导数
a2=polyder(a1);%多项式二阶导数
t=polyval(a,x); hold on;
plot(x,y,'b');%绘制贝塞尔曲线
plot(x,t,'r:');%绘制拟合曲线
hold off;
xValue=str2num(get(handles.edit_xValue,'string')) ;
yValue=polyval(a,xValue);%函数值
dy1Value=polyval(a1,xValue);%一阶导数值
dy2Value=polyval(a2,xValue);%二阶导数值
set(handles.edit_yValue,'string',num2str(yValue));%显示
set(handles.edit_dy1Value,'string',num2str(dy1Value));
set(handles.edit_dy2Value,'string',num2str(dy2Value));
end
syms x y;%定义符号
steps=str2num(get(handles.edit_steps,'string')) ;%获取阶数数组
if length(steps)==1%必须是单个曲线
switch get(handles.popupmenu_pick,'value')%获取贝塞尔函数的类型
case 1 , y=besselj(steps(1),x);
case 2, y=bessely(steps(1),x);
end
dy=diff(y);%一阶导函数
ddy=diff(dy);%二阶导函数
x=0.01:0.01:14; hold on;
plot(x,subs(dy),'m'); plot(x,subs(ddy),'g'); hold off;
end
