我要加入 登录
声振论坛 返回首页

ChaChing的个人空间 http://home.vibunion.com/?63979 [收藏] [复制] [分享] [RSS]

日志

matlab source

已有 2297 次阅读2011-12-27 11:15 |

搜集至hyyly520的百度空间的资料 http://hi.baidu.com/hyyly520/blog/category/matlab%20%B1%E0%B3%CC  (2011.12.27)
 
雅可比(jacobian)迭代法求线性方程组Ax=b的解  2011-10-29 21:17
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};
end 

D=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

matlab 文件夹中所有的文件路径 包括子文件夹的文件路径  2011-05-28 10:24
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
二维矩阵 求极值点 2010-11-06 13:04
% 一维极值点
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

求离散数据极值(点)  2010-11-05 12:56

%方法一
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); %方法2

CC = 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 实现二维和三维隐函数绘图 2010-09-20 19:02
自定义二维隐函数绘图
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;
end
if 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 by

x=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; end

2.隐函数曲面表面图函数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

MATLAB 设置axes背景为图片  2010-09-20 20:29

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背景为图片')

MATLAB 贝塞尔函数零点  2010-07-21 23:14

   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

Matlab 通过多项式拟合求贝塞尔函数的导函数  2010-07-21 23:03

%产生拟合曲线,并求某点导数
% 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

MATLAB 直接求贝塞尔函数的导函数  2010-07-21 23:01

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

 

评论 (0 个评论)

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 我要加入

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-28 20:08 , Processed in 0.042036 second(s), 15 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部