精华内容
下载资源
问答
  • ![图片说明](https://img-ask.csdn.net/upload/201705/31/1496218447_927648.png) 就是这个mu,代码太庞大了,我也不知道怎么一个一个找这个mu param这个参数如图所示 ![图片说明]...
  • 但是,如果使用分号结束行,MATLAB 会执行计算,但不会显示任何输出。当生成大型矩阵时,此功能尤其有用。例如, A = magic(100); 输入长语句 如果语句无法容纳在一行中,请使用省略号(三个句点)...,后跟...

    取消输出

    如果您在仅键入语句后按 Return 或 Enter,MATLAB 会在屏幕上自动显示结果。但是,如果使用分号结束行,MATLAB 会执行计算,但不会显示任何输出。当生成大型矩阵时,此功能尤其有用。例如,

    A = magic(100);

    输入长语句

    如果语句无法容纳在一行中,请使用省略号(三个句点)...,后跟 Return 或 Enter 以指示该语句在下一行继续。例如,

    s = 1 -1/2 + 1/3 -1/4 + 1/5 - 1/6 + 1/7 ...
          - 1/8 + 1/9 - 1/10 + 1/11 - 1/12;

    =+ 和 - 符号周围的空白是可选的,但可提高可读性。

    命令行编辑

    使用键盘上的各个箭头键和控制键可以重新调用、编辑和重用先前键入的语句。例如,假定您错误地输入了

    rho = (1 + sqt(5))/2

    sqrt 的拼写不正确。MATLAB 会给出以下错误信息

    Undefined function 'sqt' for input arguments of type 'double'.

    下标

    A 的行 i 和列 j 中的元素通过 A(i,j) 表示。例如,A(4,2) 表示第四行和第二列中的数字。在幻方矩阵中,A(4,2) 为 15。因此,要计算 A 第四列中的元素的总和,请键入

    A(1,4) + A(2,4) + A(3,4) + A(4,4)

    此下标生成

    ans =
         34

    但这不是计算某列总和的最佳方法。

    此外,还可以使用单一下标 A(k) 引用矩阵的元素。单一下标是引用行和列向量的常见方法。但是,也可以对满二维矩阵应用单一下标。在这种情况下,数组被视为一个由原始矩阵的列构成的长列向量。因此,在幻方矩阵中,A(8) 是另一种引用存储在 A(4,2) 中的值 15 的方法。

    如果尝试使用矩阵外部元素的值,则会生成错误:

    t = A(4,5)
    
    索引超出矩阵维度。

    相反,如果将值存储在矩阵外部元素中,则会增大大小以便容纳新元素:

    X = A;
    X(4,5) = 17
    
    X =
        16     3     2    13     0
         5    10    11     8     0
         9     6     7    12     0
         4    15    14     1    17

    冒号运算符

    冒号 : 是最重要的 MATLAB® 运算符之一。它以多种不同形式出现。表达式

    1:10

    是包含从 1 到 10 之间的整数的行向量:

    1     2     3     4     5     6     7     8     9    10

    要获取非单位间距,请指定增量。例如,

    100:-7:50

    100    93    86    79    72    65    58    51

    0:pi/4:pi

    0    0.7854    1.5708    2.3562    3.1416

    包含冒号的下标表达式引用部分矩阵:

    A(1:k,j)

    表示 A 第 j 列中的前 k 个元素。因此,

    sum(A(1:4,4))

    计算第四列的总和。但是,执行此计算有一种更好的方法。冒号本身引用矩阵行或列中的所有元素,而关键字 end 引用最后一个行或列。因此,

    sum(A(:,end))

    计算 A 最后一列中的元素的总和:

    ans =
         34

    为什么 4×4 幻方矩阵的幻数和等于 34?如果将介于 1 到 16 之间的整数分为四个总和相等的组,该总和必须为

    sum(1:16)/4

    当然,也即

    ans =
         34

    串联

    串联是连接小矩阵以便形成更大矩阵的过程。实际上,第一个矩阵是通过将其各个元素串联起来而构成的。成对的方括号 [] 即为串联运算符。例如,从 4×4 幻方矩阵 A 开始,组成

    B = [A  A+32; A+48  A+16]

    结果会生成一个 8×8 矩阵,这是通过连接四个子矩阵获得的:

    B =
    
        16     3     2    13    48    35    34    45
         5    10    11     8    37    42    43    40
         9     6     7    12    41    38    39    44
         4    15    14     1    36    47    46    33
        64    51    50    61    32    19    18    29
        53    58    59    56    21    26    27    24
        57    54    55    60    25    22    23    28
        52    63    62    49    20    31    30    17

    此矩阵是一个接近于幻方矩阵的矩阵。此矩阵的元素是经过重新排列的整数 1:64。此矩阵的列总和符合 8×8 幻方矩阵的要求:

    sum(B)
    
    ans =
       260   260   260   260   260   260   260   260

    但是其行总和 sum(B')' 并不完全相同。要使其成为有效的 8×8 幻方矩阵,需要进行进一步操作。

    删除行和列

    只需使用一对方括号即可从矩阵中删除行和列。首先

    X = A;

    然后,要删除 X 的第二列,请使用

    X(:,2) = []

    这会将 X 更改为

    X =
        16     2    13
         5    11     8 
         9     7    12
         4    14     1

    如果您删除矩阵中的单个元素,结果将不再是矩阵。因此,以下类似表达式

    X(1,2) = []

    将会导致错误。但是,使用单一下标可以删除一个元素或元素序列,并将其余元素重构为一个行向量。因此

    X(2:2:10) = []

    生成

    X =
        16     9     2     7    13    12     1

    标量扩展

    可以采用多种不同方法将矩阵和标量合并在一起。例如,通过从每个元素中减去标量而将其从矩阵中减去。幻方矩阵的元素平均值为 8.5,因此

    B = A - 8.5

    形成一个列总和为零的矩阵:

    B =
          7.5     -5.5     -6.5      4.5
         -3.5      1.5      2.5     -0.5
          0.5     -2.5     -1.5      3.5
         -4.5      6.5      5.5     -7.5
    
    sum(B)
    
    ans =
         0     0     0     0

    通过标量扩展,MATLAB 会为范围中的所有索引分配一个指定标量。例如,

    B(1:2,2:3) = 0

    将 B 的某个部分清零:

    B =
          7.5        0        0      4.5
         -3.5        0        0     -0.5
          0.5     -2.5     -1.5      3.5
         -4.5      6.5      5.5     -7.5

    逻辑下标

    根据逻辑和关系运算创建的逻辑向量可用于引用子数组。假定 X 是一个普通矩阵,L 是一个由某个逻辑运算生成的同等大小的矩阵。那么,X(L) 指定 X 的元素,其中 L 的元素为非零。

    通过将逻辑运算指定为下标表达式,可以在一个步骤中完成这种下标。假定您具有以下数据集:

    x = [2.1 1.7 1.6 1.5 NaN 1.9 1.8 1.5 5.1 1.8 1.4 2.2 1.6 1.8];

    NaN 是用于缺少的观测值的标记,例如,无法响应问卷中的某个项。要使用逻辑索引删除缺少的数据,请使用 isfinite(x),对于所有有限数值,该函数为 true;对于 NaN 和 Inf,该函数为 false:

    x = x(isfinite(x))
    x =
      2.1 1.7 1.6 1.5 1.9 1.8 1.5 5.1 1.8 1.4 2.2 1.6 1.8

    现在,存在一个似乎与其他项很不一样的观测值,即 5.1。这是一个离群值。下面的语句可删除离群值,在本示例中,即比均值大三倍标准差的元素:

    x = x(abs(x-mean(x)) <= 3*std(x))
    x =
      2.1 1.7 1.6 1.5 1.9 1.8 1.5 1.8 1.4 2.2 1.6 1.8

    标量扩展对于另一示例,请使用逻辑索引和标量扩展将非质数设置为 0,以便高亮显示丢勒幻方矩阵中的质数的位置。(请参阅 magic 函数。)

    A(~isprime(A)) = 0
    
    A =
         0     3     2    13
         5     0    11     0
         0     0     7     0
         0     0     0     0

    find 函数

    find 函数可用于确定与指定逻辑条件相符的数组元素的索引。find 以最简单的形式返回索引的列向量。转置该向量以便获取索引的行向量。例如,再次从丢勒的幻方矩阵开始。(请参阅 magic 函数。)

    k = find(isprime(A))'

    使用一维索引选取幻方矩阵中的质数的位置:

    k =
         2     5     9    10    11    13

    使用以下命令按 k 确定的顺序将这些质数显示为行向量

    A(k)
    
    ans =
         5     3     2    11     7    13

    将 k 用作赋值语句的左侧索引时,会保留矩阵结构:

    A(k) = NaN
    
    A =
        16   NaN   NaN   NaN
       NaN    10   NaN     8
         9     6   NaN    12
         4    15    14     1
    展开全文
  • 输入, a : 一个 Matlab 结构体,例如 a(1).name='red', a(2).name='blue'; field : 搜索字段的名称,例如 'name' value : 搜索,例如 'blue' 输出, index : 与搜索匹配的结构索引 例子, a(1).name='\u84dd\u...
  • % 查找集合中某个元素的位置 C = {'xlh','gyl','xyh'}; c = strcmp(C,'gyl');...% Struct结构体某个字段的位置 S = repmat(struct('name',[],'age',[],'sex',[]),3,1);S(1).name = 'xlh'; S(1).age = '20
    
    % 查找集合中某个元素的位置
    C = {'xlh','gyl','xyh'};
    c = strcmp(C,'gyl'); % logical array
    ind0 = find(c==1);
    
    % Struct结构体某个字段值的位置 
    S = repmat(struct('name',[],'age',[],'sex',[]),3,1);
    
    S(1).name = 'xlh'; S(1).age = '20';S(1).sex = '男';
    S(2).name = 'gyl'; S(2).age = '19';S(2).sex = '女';
    S(3).name = 'xyh'; S(3).age = '5' ;S(3).sex = '男';
    
    s = strcmp({S.name},'gyl');% {S.name} :convert to cell
    ind = find(s==1);
    展开全文
  • 输出连续的深度神经网络 1. 前言  本文设计并实现了输出连续的深度神经网络。可用于自动构图特征线位置判断等需要连续的场合。特征线的位置可以是垂直的,也可以是水平的。即特征线有垂直和水平两种特征线。...

    输出连续值的深度神经网络

    1.     前言

            本文设计并实现了输出连续值的深度神经网络。可用于自动构图特征线位置判断等需要连续值的场合。特征线的位置可以是垂直的,也可以是水平的。即特征线有垂直和水平两种特征线。

          所设计的深度神经网络基于Deeplearning Toolbox中的神经神经网络修改而来。主要修改内容为将输出节点数由原来的10个调整为3个。每个节点的输出值也进行重新定义。原节点输出值表示属于某个分类的概率值。而新的节点输出值表示特征线所位于的图像中位置在整幅图像中的比例值。

           深度神经网络共5层,包含两层卷积层。每层卷积层后会加一层sigmoid函数层作为卷积层的激活函数。

          本文所涉及的代码已经在UFLDL教程的相关它模块中被调试通过并保存。并且与CNN_matlab_code一文中的代码相对应。

    2.     经验总结

    深度神经网络就像是一个具备很强学习能力的学生,对深度神经网络进行训练的样本就像是教学生用的教材。选用什么教材,能教会学生什么样的技能。在教学生的过程中,还要注意方式方法。对于同样的教材和同样的学生,如果采用不同的教学方式方法,可能获得完全不同的教学结果。

    对于自动构图而言,如果需要使得同一个深度神经网络能识别水平和垂直两种特征线位置,首先要收集足够的具备两种特征线的图片样本。其次就是设计3个输出节点的输出值的含义。有两种方案,能取得截然不同的训练效果。第一种方案:令节点1表示水平特征线位置在整幅图像中的比例,令节点2表示垂直特征线位置在整幅图像中的比例,那么训练效果会很差;节点1和节点2的输出值无法准确反映水平或垂直特征线位置在整幅图像中的比例。第二种方案:令节点1区分水平特征线还是垂直特征线;如果节点1输出值为1,那么表示是水平特征线;如果节点1输出值为0,那么表示是垂直特征线。而节点2所输出的值表示特征线的位置在整幅图像中的比例;当节点1输出的值为1时,节点2输出值表示水平特征线的位置在整幅图像中的比例;当节点1输出的值为0时,节点2输出值表示垂直特征线的位置在整幅图像中的比例;这种训练效果良好;节点1输出的值能准确反映图像中包含水平特征线还是垂直特征线;节点2输出的值能准确反映图像中包含水平特征线或垂直特征线的位置在图像中所处的比例。

    对于深度神经网络,如果出现过拟合现象,那么可减少深度神经网络的层数,从而减少网络中可调整的参数,使得神经网络能正确处理训练样本之外的样本。如果出现欠拟合现象,那么可增加深度神经网络的层数,从而增加网络中可调整的参数,使得神经网络能准确拟合所有的训练样本。

    当深度网络的层数较少时,sigmoid函数能很好地使得网络收敛。当深度网络的层数较多时,sigmoid函数无法使得网络很好地收敛。此时需要采用Relu等改进的激活函数以及层之间的归一化处理等方式,使得深度神经网络尽快收敛。

    sample3_1.jpg 至 sample3_8.jpg如下图所示:

             

     上述子图像从左到右依次为sample3_1.jpg 至sample3_8.jpg。图像的构图的水平特征线依次位于整幅图像的1/8处、2/8处、3/8处、... ,7/8处,8/8处。

     

    3.     代码更新部分

    3.1   测试代码

    %1. 2017.04.11:perfect test results: the continuous output results incided withtheir

    %exptected values.

    %2. 2017.04.12:test my samples

    %3. horitizatal or vertical characteritic identification;

    %4. reduce the cnn parameters to solve the problem of overfitting

    %5. reduce the precision of parameter to single precision

    clearall

    clc

    closeall

     

    fork=1:8

        temp=imread(['img_sample_h',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k)=img1(:,:,1);

        train_y(1:3,k)=single(k)/10;

    end

     

    fork=1:8

        temp=imread(['sample_folder\sample3_',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k+8)=img1(:,:,1);

    %     train_y(1:3,k)=double(k)/10;

        train_y(1:3,k+8)=single(k)/10;

    end

     

    fork=1:8

        temp=imread(['sample_folder\sample_',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k+16)=img1(:,:,1);

    %     train_y(1:3,k)=double(k)/10;

        train_y(1:3,k+16)=single(k)/10;

    end

     

    fork=1:8

        temp=imread(['sample_folder\sample4_',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k+24)=img1(:,:,1);

    %     train_y(1:3,k)=double(k)/10;

        train_y(1:3,k+24)=single(k)/10;

    end

     

    train_x1= single(reshape(train_x,96,96,32))/255;

    train_y1= single(train_y);

     

    figure

    fori=1:32

    subplot(4,8,i)

    imshow(train_x1(:,:,i))

    title(sprintf('h=%.2f',train_y1(1,i)))

    end

     

    %v_cnn=cnn

    load cnn_epoch_500.mat

    %num_train_samples=30000;

    opts.alpha= 1;

    opts.batchsize=4;

    %opts.numepochs = 1;

    opts.numepochs=1000;

     

    %cnn = cnntrain(cnn, train_x1, train_y1, opts);

    cnn= cnntrain(cnn, train_x1, train_y1, opts);

    save('cnn_epoch_500.mat','cnn')

     

     

    %test_x(:,:,1)=imread('img_sample_h3.jpg');

    %test_x(:,:,2)=imread('img_sample_h34.jpg');

    %test_x(:,:,3)=imread('img_sample_h37.jpg');

    %test_x(:,:,4)=imread('img_sample_h4.jpg');

     

    test_x(:,:,1)=imread('sample_folder\sample3_4.jpg');

    test_x(:,:,2)=imread('sample_folder\sample3_8.jpg');

    test_x(:,:,3)=imread('sample_folder\sample_6.jpg');

    test_x(:,:,4)=imread('sample_folder\sample_3.jpg');

     

    test_y(1:3,1)=0.3;

    test_y(1:3,2)=0.34;

    test_y(1:3,3)=0.37;

    test_y(1:3,4)=0.4;

     

    %test_x1 = double(reshape(test_x,28,28,4))/255;

    %test_y1 = double(test_y);

    test_x1= single(reshape(test_x,96,96,4))/255;

    test_y1= single(test_y);

    [er,bad] = cnntest(cnn, test_x1, test_y1);

    3.2 训练部分

    % 1. 2017.04.11: perfecttest results: the continuous output results incided with their expected values.

    % 2. 2017.04.12:test mysamples

    % 3. horitizatal orvertical characteritic identification;

    % 4. reduce the cnnparameters to solve the problem of overfitting

    % 5. reduce the precisionof parameter to single precision

    clear all

    clc

    close all

     

    for k=1:8

        temp=imread(['img_sample_h',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k)=img1(:,:,1);

    %     train_y(1:3,k)=double(k)/10;

        train_y(1:3,k)=single(k)/10;

    end

     

    for k=1:8

        temp=imread(['sample_folder\sample3_',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k+8)=img1(:,:,1);

    %     train_y(1:3,k)=double(k)/10;

        train_y(1:3,k+8)=single(k)/10;

    end

     

    for k=1:8

        temp=imread(['sample_folder\sample_',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k+16)=img1(:,:,1);

    %     train_y(1:3,k)=double(k)/10;

        train_y(1:3,k+16)=single(k)/10;

    end

     

    for k=1:8

        temp=imread(['sample_folder\sample4_',num2str(k),'.jpg']);

        img1=repmat(temp,1,1,2);

        train_x(:,:,k+24)=img1(:,:,1);

    %     train_y(1:3,k)=double(k)/10;

        train_y(1:3,k+24)=single(k)/10;

    end

     

    train_x1 =single(reshape(train_x,96,96,32))/255;

    train_y1 = single(train_y);

     

    figure

    for i=1:32

    subplot(4,8,i)

    imshow(train_x1(:,:,i))

    title(sprintf('h=%.2f',train_y1(1,i)))

    end

     

    rand('state',0)

    % cnn.layers = {

    %     struct('type', 'i') %input layer

    %     struct('type', 'c', 'outputmaps', 6, 'kernelsize',5) %convolution layer

    %     struct('type', 's', 'scale', 2) %subsampling layer

    %     struct('type', 'c', 'outputmaps', 12,'kernelsize', 5) %convolution layer

    %     struct('type', 's', 'scale', 2)%subsampling layer   

    % };

     

    cnn.layers = {

        struct('type', 'i') %input layer

        struct('type', 'c', 'outputmaps', 6, 'kernelsize', 5) %convolutionlayer

        struct('type', 's', 'scale', 2) %sub samplinglayer

        struct('type', 'c', 'outputmaps', 12, 'kernelsize', 5) %convolutionlayer

        struct('type', 's', 'scale', 2) %subsamplinglayer   

        struct('type', 'c', 'outputmaps', 7, 'kernelsize',4) %convolutionlayer

        struct('type', 's', 'scale', 3) %subsamplinglayer   

    };

     

    cnn = cnnsetup(cnn, train_x,train_y);

     

    opts.alpha = 1;

    opts.batchsize =4;

    opts.numepochs =1000;

     

    cnn = cnntrain(cnn, train_x1,train_y1, opts);

    save('cnn_epoch_500.mat','cnn')

     

    test_x(:,:,1)=imread('sample_folder\sample3_4.jpg');

    test_x(:,:,2)=imread('sample_folder\sample3_8.jpg');

    test_x(:,:,3)=imread('sample_folder\sample_6.jpg');

    test_x(:,:,4)=imread('sample_folder\sample_3.jpg');

     

    test_y(1:3,1)=0.3;   % the number of output nodes is equal to 3.

    test_y(1:3,2)=0.34;

    test_y(1:3,3)=0.37;

    test_y(1:3,4)=0.4;

     

    test_x1 =single(reshape(test_x,96,96,4))/255;

    test_y1 = single(test_y);

     

    [er, bad] = cnntest(cnn,test_x1, test_y1);

     

    3.2   cnnsetup

    这部分中的onum应该修改为输出节点的数目。

    function net = cnnsetup(net, x, y)

        assert(~isOctave() ||compare_versions(OCTAVE_VERSION, '3.8.0', '>='), ['Octave 3.8.0 orgreater is required for CNNs as there is a bug in convolution in previousversions. See http://savannah.gnu.org/bugs/?39314. Your version is ' myOctaveVersion]);

        inputmaps = 1;

        mapsize = size(squeeze(x(:, :, 1)));

      

        for l = 1 : numel(net.layers)  %  layer

            ifstrcmp(net.layers{l}.type,'s')

                mapsize = mapsize /net.layers{l}.scale;

               assert(all(floor(mapsize)==mapsize), ['Layer ' num2str(l)' size must be integer. Actual: 'num2str(mapsize)]);

                for j = 1 : inputmaps

                    net.layers{l}.b{j} = 0;

                end

            end

            if strcmp(net.layers{l}.type,'c')

                mapsize = mapsize -net.layers{l}.kernelsize + 1;

                fan_out = net.layers{l}.outputmaps* net.layers{l}.kernelsize ^ 2;

                for j = 1 : net.layers{l}.outputmaps %  output map

                    fan_in = inputmaps *net.layers{l}.kernelsize ^ 2;

                    for i = 1 : inputmaps %  input map

                        net.layers{l}.k{i}{j} =(rand(net.layers{l}.kernelsize) - 0.5) * 2 * sqrt(6 / (fan_in + fan_out));

                    end

                    net.layers{l}.b{j} = 0;

                end

                inputmaps =net.layers{l}.outputmaps;

            end

        end

        %'onum' is the number of labels, that's why it is calculated using size(y, 1).If you have 20 labels so the output of the network will be 20 neurons.

        %'fvnum' is the number of output neurons at the last layer, the layer justbefore the output layer.

        %'ffb' is the biases of the output neurons.

        %'ffW' is the weights between the last layer and the output neurons. Note thatthe last layer is fully connected to the output layer, that's why the size ofthe weights is (onum * fvnum)

        fvnum = prod(mapsize) * inputmaps;

    %     onum =size(y, 1);

         % the onum is same with thenout in the test_example_CNN2.m

        onum = 3;

     

        net.ffb = zeros(onum, 1);

        net.ffW = (rand(onum, fvnum) - 0.5) * 2 *sqrt(6 / (onum + fvnum));

    end

    4.     采用Relu激活函数代替Sigmoid函数

    4.1 cnnff

    function net = cnnff(net, x)

        n = numel(net.layers);

        net.layers{1}.a{1} = x;

        inputmaps = 1;

     

        for l = 2 : n  % for each layer

            ifstrcmp(net.layers{l}.type,'c')

                %  !!below canprobably be handled by insane matrix operations

                for j = 1 :net.layers{l}.outputmaps  % for each output map

                    %  create tempoutput map

                    z = zeros(size(net.layers{l -1}.a{1}) - [net.layers{l}.kernelsize - 1 net.layers{l}.kernelsize - 1 0]);

                    for i = 1 :inputmaps  %  for each inputmap

                        %  convolve withcorresponding kernel and add to temp output map

                        z = z + convn(net.layers{l- 1}.a{i}, net.layers{l}.k{i}{j},'valid');

                    end

                    %  add bias, passthrough nonlinearity

    %                 net.layers{l}.a{j} = sigm(z +net.layers{l}.b{j});

                      net.layers{l}.a{j} =sigm_new(z + net.layers{l}.b{j});

                end

                %  set number ofinput maps to this layers number of outputmaps

                inputmaps =net.layers{l}.outputmaps;

            elseifstrcmp(net.layers{l}.type,'s')

                %  downsample

                for j = 1 :inputmaps

                    z = convn(net.layers{l -1}.a{j}, ones(net.layers{l}.scale) / (net.layers{l}.scale ^ 2),'valid');   % !! replace with variable

                    net.layers{l}.a{j} = z(1 :net.layers{l}.scale : end, 1 : net.layers{l}.scale : end, :);

                end

            end

        end

     

        %  concatenate all end layer featuremaps into vector

        net.fv = [];

        for j = 1 : numel(net.layers{n}.a)

            sa = size(net.layers{n}.a{j});

            net.fv = [net.fv;reshape(net.layers{n}.a{j}, sa(1) * sa(2), sa(3))];

        end

        %  feedforward into outputperceptrons

        net.o = sigm(net.ffW * net.fv +repmat(net.ffb, 1, size(net.fv, 2)));

    end

     

    4.2 cnnbp

    function net = cnnbp(net, y)

        n = numel(net.layers);

        net.e = net.o - y;

       

        %  loss function

        net.L = 1/2* sum(net.e(:) .^ 2) /size(net.e, 2);

     

        %% backprop deltas

        net.od = net.e .* (net.o .* (1 -net.o));   %  output delta

        net.fvd = (net.ffW' * net.od);           %  feature vectordelta

        if strcmp(net.layers{n}.type,'c')      % only conv layers has sigm function

              net_fv=zeros(size(net.fv));

              net_fv(net.fv>0)=1;

              net_fv(net.fv==0)=0.5;

              net.fvd = net.fvd .* net_fv;

        end

     

        %  reshape feature vector deltasinto output map style

        sa = size(net.layers{n}.a{1});

        fvnum = sa(1) * sa(2);

        for j = 1 : numel(net.layers{n}.a)

            net.layers{n}.d{j} =reshape(net.fvd(((j - 1) * fvnum + 1) : j * fvnum, :), sa(1), sa(2), sa(3));

        end

     

        for l = (n - 1) : -1 : 1

            ifstrcmp(net.layers{l}.type,'c')

                for j = 1 :numel(net.layers{l}.a)

                    net.layers{l}.d{j} =net.layers{l}.a{j} .* (1 - net.layers{l}.a{j}) .* (expand(net.layers{l +1}.d{j}, [net.layers{l + 1}.scale net.layers{l + 1}.scale 1]) / net.layers{l +1}.scale ^ 2);

                end

            elseifstrcmp(net.layers{l}.type,'s')

                for i = 1 :numel(net.layers{l}.a)

                    z =zeros(size(net.layers{l}.a{1}));

                    for j = 1 :numel(net.layers{l + 1}.a)

                         z = z + convn(net.layers{l+ 1}.d{j}, rot180(net.layers{l + 1}.k{i}{j}),'full');

                    end

                    net.layers{l}.d{i} = z;

                end

            end

        end

     

        %% calc gradients

        for l = 2 : n

            ifstrcmp(net.layers{l}.type,'c')

                for j = 1 : numel(net.layers{l}.a)

                    for i = 1 :numel(net.layers{l - 1}.a)

                        net.layers{l}.dk{i}{j} =convn(flipall(net.layers{l - 1}.a{i}), net.layers{l}.d{j},'valid') / size(net.layers{l}.d{j}, 3);

                    end

                    net.layers{l}.db{j} =sum(net.layers{l}.d{j}(:)) / size(net.layers{l}.d{j}, 3);

                end

            end

        end

        net.dffW = net.od * (net.fv)' /size(net.od, 2);

        net.dffb = mean(net.od, 2);

     

        function X = rot180(X)

            X = flipdim(flipdim(X, 1), 2);

        end

    end

     

    4.3 sigm_new

    function X = sigm_new(P)

        X = max(P,0);

    end

     

    5.     Appendix

    Deep Learning论文笔记之(五)CNN卷积神经网络代码理解

    zouxy09@qq.com

    http://blog.csdn.net/zouxy09

     

             自己平时看了一些论文,但老感觉看完过后就会慢慢的淡忘,某一天重新拾起来的时候又好像没有看过一样。所以想习惯地把一些感觉有用的论文中的知识点总结整理一下,一方面在整理过程中,自己的理解也会更深,另一方面也方便未来自己的勘察。更好的还可以放到博客上面与大家交流。因为基础有限,所以对论文的一些理解可能不太正确,还望大家不吝指正交流,谢谢。

     

          本文的代码来自githupDeep Learning的toolbox,(在这里,先感谢该toolbox的作者)里面包含了很多Deep Learning方法的代码。是用Matlab编写的(另外,有人翻译成了C++Python的版本了)。本文中我们主要解读下CNN的代码。详细的注释见代码。

           在读代码之前,最好先阅读下我的上一个博文:

            Deep Learning论文笔记之(四)CNN卷积神经网络推导和实现

                http://blog.csdn.net/zouxy09/article/details/9993371

          里面包含的是我对一个作者的CNN笔记的翻译性的理解,对CNN的推导和实现做了详细的介绍,看明白这个笔记对代码的理解非常重要,所以强烈建议先看懂上面这篇文章。

     

             下面是自己对代码的注释:

    cnnexamples.m

    [plain] view plain copy

    1.              clear all; close all; clc;  

    2.              addpath('../data');  

    3.              addpath('../util');  

    4.              load mnist_uint8;  

    5.                

    6.              train_x = double(reshape(train_x',28,28,60000))/255;  

    7.              test_x = double(reshape(test_x',28,28,10000))/255;  

    8.              train_y = double(train_y');  

    9.              test_y = double(test_y');  

    10.             

    11.           %% ex1   

    12.           %will run 1 epoch in about 200 second and get around 11% error.   

    13.           %With 100 epochs you'll get around 1.2% error  

    14.             

    15.           cnn.layers = {  

    16.               struct('type', 'i') %input layer  

    17.               struct('type', 'c', 'outputmaps', 6, 'kernelsize', 5) %convolution layer  

    18.               struct('type', 's', 'scale', 2) %sub sampling layer  

    19.               struct('type', 'c', 'outputmaps', 12, 'kernelsize', 5) %convolution layer  

    20.               struct('type', 's', 'scale', 2) %subsampling layer  

    21.           };  

    22.             

    23.           % 这里把cnn的设置给cnnsetup,它会据此构建一个完整的CNN网络,并返回  

    24.           cnn = cnnsetup(cnn, train_x, train_y);  

    25.             

    26.           % 学习率  

    27.           opts.alpha = 1;  

    % 每次挑出一个batchsize的batch来训练,也就是每用batchsize个样本就调整一次权值,而不是把所有样本都输入了,计算所有样本的误差了才调整一次权值  

    28.           opts.batchsize = 50;   

    29.           % 训练次数,用同样的样本集。我训练的时候:  

    30.           % 1的时候 11.41% error  

    31.           % 5的时候 4.2% error  

    32.           % 10的时候 2.73% error  

    33.           opts.numepochs = 10;  

    34.             

    35.           % 然后开始把训练样本给它,开始训练这个CNN网络  

    36.           cnn = cnntrain(cnn, train_x, train_y, opts);  

    37.             

    38.           % 然后就用测试样本来测试  

    39.           [er, bad] = cnntest(cnn, test_x, test_y);  

    40.             

    41.           %plot mean squared error  

    42.           plot(cnn.rL);  

    43.           %show test error  

    44.           disp([num2str(er*100) '% error']);  

     

    cnnsetup.m

    [plain] view plain copy

    1.              function net = cnnsetup(net, x, y)  

    2.              inputmaps = 1;  

    3.              % B=squeeze(A) 返回和矩阵A相同元素但所有单一维都移除的矩阵B,单一维是满足size(A,dim)=1的维。  

    4.                  % train_x中图像的存放方式是三维的reshape(train_x',28,28,60000),前面两维表示图像的行与列,  

    5.                  % 第三维就表示有多少个图像。这样squeeze(x(:, :, 1))就相当于取第一个图像样本后,再把第三维  

    6.                  % 移除,就变成了28x28的矩阵,也就是得到一幅图像,再size一下就得到了训练样本图像的行数与列数了  

    7.                  mapsize = size(squeeze(x(:, :, 1)));  

    8.                

    % 下面通过传入net这个结构体来逐层构建CNN网络  

    % n = numel(A)返回数组A中元素个数  

    % net.layers中有五个struct类型的元素,实际上就表示CNN共有五层,这里范围的是5  

    9.                  for l = 1 : numel(net.layers)   %  layer  

    10.                   if strcmp(net.layers{l}.type, 's') % 如果这层是 子采样层  

    11.           % subsampling层的mapsize,最开始mapsize是每张图的大小28*28  

    12.           % 这里除以scale=2,就是pooling之后图的大小,pooling域之间没有重叠,所以pooling后的图像为14*14, 注意这里的右边的mapsize保存的都是上一层每张特征map的大小,它会随着循环进行不断更新  

    13.                       mapsize = floor(mapsize / net.layers{l}.scale);  

    14.                       for j = 1 : inputmaps % inputmap就是上一层有多少张特征图  

    15.                           net.layers{l}.b{j} = 0; % 将偏置初始化为0  

    16.                       end  

    17.                   end  

    18.                   if strcmp(net.layers{l}.type, 'c') % 如果这层是 卷积层  

    19.           % 旧的mapsize保存的是上一层的特征map的大小,那么如果卷积核的移动步长是1,那用kernelsize*kernelsize大小的卷积核卷积上一层的特征map后,得到的新的map的大小就是下面这样  

    20.                       mapsize = mapsize - net.layers{l}.kernelsize + 1;  

    % 该层需要学习的参数个数。每张特征map是一个(后层特征图数量)*(用来卷积的patch图的大小)。 因为是通过用一个核窗口在上一个特征map层中移动(核窗口每次移动1个像素),遍历上一个特征map 层的每个神经元。核窗口由kernelsize*kernelsize个元素组成,每个元素是一个独立的权值,所以就有kernelsize*kernelsize个需要学习的权值,再加一个偏置值。另外,由于是权值共享,也就是说同一个特征map层是用同一个具有相同权值元素的kernelsize*kernelsize的核窗口去感受输入上一个特征map层的每个神经元得到的,所以同一个特征map,它的权值是一样的,共享的,权值只取决于核窗口。然后,不同的特征map提取输入上一个特征map层不同的特征,所以采用的核窗口不一样,也就是权值不一样,所以outputmaps个特征map就有(kernelsize*kernelsize+1)* outputmaps那么多的权值了。但这里fan_out只保存卷积核的权值W,偏置b在下面独立保存。

    21.                       fan_out = net.layers{l}.outputmaps * net.layers{l}.kernelsize ^ 2;  

    22.                       for j = 1 : net.layers{l}.outputmaps  %  output map  

    23.            % fan_out保存的是对于上一层的一张特征map,我在这一层需要对这一张特征map提取outputmaps种特征,提取每种特征用到的卷积核不同,所以fan_out保存的是这一层输出新的特征需要学习的参数个数。而,fan_in保存的是,我在这一层,要连接到上一层中所有的特征map,然后用fan_out保存的提取特征的权值来提取他们的特征。也即是对于每一个当前层特征图,有多少个参数链到前层  

    24.                           fan_in = inputmaps * net.layers{l}.kernelsize ^ 2;  

    25.                           for i = 1 : inputmaps  %  input map  

    26.           % 随机初始化权值,也就是共有outputmaps个卷积核,对上层的每个特征map,都需要用这么多个卷积核  

    27.           % 去卷积提取特征。  

    28.           % rand(n)是产生n×n的 0-1之间均匀取值的数值的矩阵,再减去0.5就相当于产生-0.5到0.5之间的随机数,再 *2 就放大到 [-1, 1]。然后再乘以后面那一数,why?  

    % 反正就是将卷积核每个元素初始化为[-sqrt(6 / (fan_in + fan_out)), sqrt(6 / (fan_in + fan_out))]之间的随机数。因为这里是权值共享的,也就是对于一张特征map,所有感受野位置的卷积核都是一样的,所以只需要保存的是 inputmaps * outputmaps 个卷积核。  

    29.           net.layers{l}.k{i}{j} = (rand(net.layers{l}.kernelsize) - 0.5) * 2 * sqrt(6 / (fan_in + fan_out));

    30.                           end  

    31.                           net.layers{l}.b{j} = 0; % 将偏置初始化为0  

    32.                       end  

    33.            % 只有在卷积层的时候才会改变特征map的个数,pooling的时候不会改变个数。这层输出的特征map个数就是输入到下一层的特征map个数  

    34.                       inputmaps = net.layers{l}.outputmaps;   

    35.                   end  

    36.               end  

    37.                 

    38.               % fvnum 是输出层的前面一层的神经元个数。  

    39.               % 这一层的上一层是经过pooling后的层,包含有inputmaps个特征map。每个特征map的大小是mapsize。所以,该层的神经元个数是 inputmaps * (每个特征map的大小)  

    40.               % prod: Product of elements.  

    41.               % For vectors, prod(X) is the product of the elements of X  

    42.               % 在这里 mapsize = [特征map的行数 特征map的列数],所以prod后就是 特征map的行*列  

    43.               fvnum = prod(mapsize) * inputmaps;  

    44.               % onum 是标签的个数,也就是输出层神经元的个数。你要分多少个类,自然就有多少个输出神经元  

    45.               onum = size(y, 1);  

    46.             

    47.               % 这里是最后一层神经网络的设定  

    48.               % ffb 是输出层每个神经元对应的基biases  

    49.               net.ffb = zeros(onum, 1);  

    50.               % ffW 输出层前一层 与 输出层 连接的权值,这两层之间是全连接的  

    51.               net.ffW = (rand(onum, fvnum) - 0.5) * 2 * sqrt(6 / (onum + fvnum));  

    52.           end  

     

    cnntrain.m

    [plain] view plain copy

    1.              function net = cnntrain(net, x, y, opts)  

    2.                  m = size(x, 3); % m 保存的是 训练样本个数  

    3.                  numbatches = m / opts.batchsize;  

    4.                  % rem: Remainder after division. rem(x,y) is x - n.*y 相当于求余  

    5.                  % rem(numbatches, 1) 就相当于取其小数部分,如果为0,就是整数  

    6.                  if rem(numbatches, 1) ~= 0  

    7.                      error('numbatches not integer');  

    8.                  end  

    9.                    

    10.               net.rL = [];  

    11.               for i = 1 : opts.numepochs  

    12.                   % disp(X) 打印数组元素。如果X是个字符串,那就打印这个字符串  

    13.                   disp(['epoch ' num2str(i) '/' num2str(opts.numepochs)]);  

    14.                   % tic 和 toc 是用来计时的,计算这两条语句之间所耗的时间  

    15.                   tic;  

    16.                   % P = randperm(N) 返回[1, N]之间所有整数的一个随机的序列,例如  

    17.                   % randperm(6) 可能会返回 [2 4 5 6 1 3]  

    18.                   % 这样就相当于把原来的样本排列打乱,再挑出一些样本来训练  

    19.                   kk = randperm(m);  

    20.                   for l = 1 : numbatches  

    21.                       % 取出打乱顺序后的batchsize个样本和对应的标签  

    22.                       batch_x = x(:, :, kk((l - 1) * opts.batchsize + 1 : l * opts.batchsize));  

    23.                       batch_y = y(:,    kk((l - 1) * opts.batchsize + 1 : l * opts.batchsize));  

    24.             

    25.                       % 在当前的网络权值和网络输入下计算网络的输出  

    26.                       net = cnnff(net, batch_x); % Feedforward  

    27.                       % 得到上面的网络输出后,通过对应的样本标签用bp算法来得到误差对网络权值  

    28.                       %(也就是那些卷积核的元素)的导数  

    29.                       net = cnnbp(net, batch_y); % Backpropagation  

    30.                       % 得到误差对权值的导数后,就通过权值更新方法去更新权值  

    31.                       net = cnnapplygrads(net, opts);  

    32.                       if isempty(net.rL)  

    33.                           net.rL(1) = net.L; % 代价函数值,也就是误差值  

    34.                       end  

    35.                       net.rL(end + 1) = 0.99 * net.rL(end) + 0.01 * net.L; % 保存历史的误差值,以便画图分析  

    36.                   end  

    37.                   toc;  

    38.               end  

    39.                 

    40.           end  

     

    cnnff.m

    [plain] view plain copy

    function net = cnnff(net, x)  

    1.               n = numel(net.layers); % 层数  

    2.               net.layers{1}.a{1} = x; % 网络的第一层就是输入,但这里的输入包含了多个训练图像  

    3.                  inputmaps = 1; % 输入层只有一个特征map,也就是原始的输入图像  

    4.                

    5.                  for l = 2 : n   %  for each layer  

    6.                      if strcmp(net.layers{l}.type, 'c') % 卷积层  

    7.                    %  !!below can probably be handled by insane matrix operations  

    8.                     % 对每一个输入map,或者说我们需要用outputmaps个不同的卷积核去卷积图像  

    9.                          for j = 1 : net.layers{l}.outputmaps   %  for each output map  

    10.                           %  create temp output map  

    11.                           % 对上一层的每一张特征map,卷积后的特征map的大小就是   

    12.                           % (输入map宽 - 卷积核的宽 + 1)* (输入map高 - 卷积核高 + 1)  

    % 对于这里的层,因为每层都包含多张特征map,对应的索引保存在每层map的第三维  

    % 所以,这里的z保存的就是该层中所有的特征map了  

    z = zeros(size(net.layers{l - 1}.a{1}) - [net.layers{l}.kernelsize - 1 net.layers{l}.kernelsize - 1 0]);                  for i = 1 : inputmaps   %  for each input map  

    13.            %  convolve with corresponding kernel and add to temp output map  

    14.            % 将上一层的每一个特征map(也就是这层的输入map)与该层的卷积核进行卷积,然后将对上一层特征map的所有结果加起来。也就是说,当前层的一张特征map,是用一种卷积核去卷积上一层中所有的特征map,然后所有特征map对应位置的卷积值的和。另外,有些论文或者实际应用中,并不是与全部的特征map链接的,有可能只与其中的某几个连接。

    15.                    z = z + convn(net.layers{l - 1}.a{i}, net.layers{l}.k{i}{j}, 'valid');  

    16.               end  

    17.            %  add bias, pass through nonlinearity  

    18.            % 加上对应位置的基b,然后再用sigmoid函数算出特征map中每个位置的激活值,作为该层输出特征map  

    19.                           net.layers{l}.a{j} = sigm(z + net.layers{l}.b{j});  

    20.                       end  

    21.                       %  set number of input maps to this layers number of outputmaps  

    22.                       inputmaps = net.layers{l}.outputmaps;  

    23.                   elseif strcmp(net.layers{l}.type, 's') % 下采样层  

    24.                       %  downsample  

    25.                       for j = 1 : inputmaps  

    26.                           %  !! replace with variable  

    27.                           % 例如我们要在scale=2的域上面执行mean pooling,那么可以卷积大小为2*2,每个元素都是1/4的卷积核  

    28.             z = convn(net.layers{l - 1}.a{j}, ones(net.layers{l}.scale) / (net.layers{l}.scale ^ 2), 'valid'); 

    29.            % 因为convn函数的默认卷积步长为1,而pooling操作的域是没有重叠的,所以对于上面的卷积结果, 最终pooling的结果需要从上面得到的卷积结果中以scale=2为步长,跳着把mean pooling的值读出来  

    30.                           net.layers{l}.a{j} = z(1 : net.layers{l}.scale : end, 1 : net.layers{l}.scale : end, :);  

    31.                       end  

    32.                   end  

    33.               end  

    34.             

    35.               %  concatenate all end layer feature maps into vector  

    36.               % 把最后一层得到的特征map拉成一条向量,作为最终提取到的特征向量  

    37.               net.fv = [];  

    38.               for j = 1 : numel(net.layers{n}.a) % 最后一层的特征map的个数  

    39.                   sa = size(net.layers{n}.a{j}); % 第j个特征map的大小  

    40.                   % 将所有的特征map拉成一条列向量。还有一维就是对应的样本索引。每个样本一列,每列为对应的特征向量  

    41.                   net.fv = [net.fv; reshape(net.layers{n}.a{j}, sa(1) * sa(2), sa(3))];  

    42.               end  

    43.               %  feedforward into output perceptrons  

    44.               % 计算网络的最终输出值。sigmoid(W*X + b),注意是同时计算了batchsize个样本的输出值  

    45.               net.o = sigm(net.ffW * net.fv + repmat(net.ffb, 1, size(net.fv, 2))); 

    46.           end  

     

    cnnbp.m

    [plain] view plain copy

    1.              function net = cnnbp(net, y)  

    2.                  n = numel(net.layers); % 网络层数  

    3.                

    4.                  %  error  

    5.                  net.e = net.o - y;   

    6.                  %  loss function  

    7.                  % 代价函数是 均方误差  

    8.                  net.L = 1/2* sum(net.e(:) .^ 2) / size(net.e, 2);  

    9.                

    10.               %%  backprop deltas  

    11.               % 这里可以参考 UFLDL 的 反向传导算法 的说明  

    12.               % 输出层的 灵敏度 或者 残差  

    13.               net.od = net.e .* (net.o .* (1 - net.o));   %  output delta  

    14.               % 残差 反向传播回 前一层  

    15.               net.fvd = (net.ffW' * net.od);              %  feature vector delta  

    16.               if strcmp(net.layers{n}.type, 'c')         %  only conv layers has sigm function  

    17.                   net.fvd = net.fvd .* (net.fv .* (1 - net.fv));  

    18.               end  

    19.             

    20.               %  reshape feature vector deltas into output map style  

    21.               sa = size(net.layers{n}.a{1}); % 最后一层特征map的大小。这里的最后一层都是指输出层的前一层  

    22.               fvnum = sa(1) * sa(2); % 因为是将最后一层特征map拉成一条向量,所以对于一个样本来说,特征维数是这样  

    23.               for j = 1 : numel(net.layers{n}.a) % 最后一层的特征map的个数  

    24.                   % 在fvd里面保存的是所有样本的特征向量(在cnnff.m函数中用特征map拉成的),所以这里需要重新  

    25.                   % 变换回来特征map的形式。d 保存的是 delta,也就是 灵敏度 或者 残差  

    26.            net.layers{n}.d{j} = reshape(net.fvd(((j - 1) * fvnum + 1) : j * fvnum, :), sa(1), sa(2), sa(3));      end  

    27.             

    28.               % 对于 输出层前面的层(与输出层计算残差的方式不同)  

    29.               for l = (n - 1) : -1 : 1  

    30.                   if strcmp(net.layers{l}.type, 'c')  

    31.                       for j = 1 : numel(net.layers{l}.a) % 该层特征map的个数  

    % net.layers{l}.d{j} 保存的是 第l层 的 第j个 map 的 灵敏度map。 也就是每个神经元节点的delta的值。expand的操作相当于对l+1层的灵敏度map进行上采样。然后前面的操作相当于对该层的输入a进行sigmoid求导  

    32.            % 这条公式请参考 Notes on Convolutional Neural Networks  

    33.            % for k = 1:size(net.layers{l + 1}.d{j}, 3)  

    34.            % net.layers{l}.d{j}(:,:,k) = net.layers{l}.a{j}(:,:,k) .* (1 - net.layers{l}.a{j}(:,:,k)) .*  kron(net.layers{l + 1}.d{j}(:,:,k), ones(net.layers{l + 1}.scale)) / net.layers{l + 1}.scale ^ 2;  

    35.            % end  

    36.             net.layers{l}.d{j} = net.layers{l}.a{j} .* (1 - net.layers{l}.a{j}) .* (expand(net.layers{l + 1}.d{j}, [net.layers{l + 1}.scale net.layers{l + 1}.scale 1]) / net.layers{l + 1}.scale ^ 2);  

    37.                       end  

    38.                   elseif strcmp(net.layers{l}.type, 's')  

    39.                       for i = 1 : numel(net.layers{l}.a) % 第l层特征map的个数  

    40.                           z = zeros(size(net.layers{l}.a{1}));  

    41.                           for j = 1 : numel(net.layers{l + 1}.a) % 第l+1层特征map的个数  

    42.                                z = z + convn(net.layers{l + 1}.d{j}, rot180(net.layers{l + 1}.k{i}{j}), 'full');  

    43.                           end  

    44.                           net.layers{l}.d{i} = z;  

    45.                       end  

    46.                   end  

    47.               end  

    48.             

    49.               %%  calc gradients  

    50.               % 这里与 Notes on Convolutional Neural Networks 中不同,这里的 子采样 层没有参数,也没有激活函数,所以在子采样层是没有需要求解的参数的  

    51.               for l = 2 : n  

    52.                   if strcmp(net.layers{l}.type, 'c')  

    53.                       for j = 1 : numel(net.layers{l}.a)  

    54.                           for i = 1 : numel(net.layers{l - 1}.a)  

    55.                               % dk 保存的是 误差对卷积核 的导数  

    56.                               net.layers{l}.dk{i}{j} = convn(flipall(net.layers{l - 1}.a{i}), net.layers{l}.d{j}, 'valid') / size(net.layers{l}.d{j}, 3);  

    57.                           end  

    58.                           % db 保存的是 误差对于bias基 的导数  

    59.                           net.layers{l}.db{j} = sum(net.layers{l}.d{j}(:)) / size(net.layers{l}.d{j}, 3);  

    60.                       end  

    61.                   end  

    62.               end  

    63.               % 最后一层perceptron的gradient的计算  

    64.               net.dffW = net.od * (net.fv)' / size(net.od, 2);  

    65.               net.dffb = mean(net.od, 2);  

    66.             

    67.               function X = rot180(X)  

    68.                   X = flipdim(flipdim(X, 1), 2);  

    69.               end  

    70.           end  

     

    cnnapplygrads.m

    [plain] view plain copy

    1.              function net = cnnapplygrads(net, opts)  

    2.                  for l = 2 : numel(net.layers)  

    3.                      if strcmp(net.layers{l}.type, 'c')  

    4.                          for j = 1 : numel(net.layers{l}.a)  

    5.                              for ii = 1 : numel(net.layers{l - 1}.a)  

    6.                                  % 这里没什么好说的,就是普通的权值更新的公式:W_new = W_old - alpha * de/dW(误差对权值导数)  

    7.                            net.layers{l}.k{ii}{j} = net.layers{l}.k{ii}{j} - opts.alpha * net.layers{l}.dk{ii}{j};  

    8.                              end  

    9.                          end  

    10.                       net.layers{l}.b{j} = net.layers{l}.b{j} - opts.alpha * net.layers{l}.db{j};  

    11.                   end  

    12.               end  

    13.             

    14.               net.ffW = net.ffW - opts.alpha * net.dffW;  

    15.               net.ffb = net.ffb - opts.alpha * net.dffb;  

    16.           end  

    cnntest.m

    [plain] view plain copy

    1.              function [er, bad] = cnntest(net, x, y)  

    2.                  %  feedforward  

    3.                  net = cnnff(net, x); % 前向传播得到输出  

    4.                  % [Y,I] = max(X) returns the indices of the maximum values in vector I  

    5.                  [~, h] = max(net.o); % 找到最大的输出对应的标签  

    6.                  [~, a] = max(y);     % 找到最大的期望输出对应的索引  

    7.                  bad = find(h ~= a);  % 找到他们不相同的个数,也就是错误的次数  

    8.                

    9.                  er = numel(bad) / size(y, 2); % 计算错误率  

    10.           end  

     

    88888888888888888888888888888888888888888888888888888888888888888888888888888888

    Deep Learning论文笔记之(四)CNN卷积神经网络推导和实现

    zouxy09@qq.com

    http://blog.csdn.net/zouxy09

           本文的论文来自:

    Notes on Convolutional NeuralNetworks,Jake Bouvrie

            这个主要是CNN的推导和实现的一些笔记,再看懂这个笔记之前,最好具有CNN的一些基础。这里也先列出一个资料供参考:

    [1] Deep Learning(深度学习)学习笔记整理系列之(七)

    [2] LeNet-5, convolutional neural networks

    [3]卷积神经网络

    [4] Neural Network for Recognition ofHandwritten Digits

    [5] Deep learning:三十八(Stacked CNN简单介绍)

    [6] Gradient-based learning appliedto document recognition.

    [7]Imagenet classification withdeep convolutional neural networks.

    [8] UFLDL中的卷积特征提取池化

           另外,这里有个matlabDeep Learning的toolbox,里面包含了CNN的代码,在下一个博文中,我将会详细注释这个代码。这个笔记对这个代码的理解非常重要。

            下面是自己对其中的一些知识点的理解:

     

    Notes on Convolutional Neural Networks

    一、介绍

            这个文档讨论的是CNNs的推导和实现。CNN架构的连接比权值要多很多,这实际上就隐含着实现了某种形式的规则化。这种特别的网络假定了我们希望通过数据驱动的方式学习到一些滤波器,作为提取输入的特征的一种方法。

            本文中,我们先对训练全连接网络的经典BP算法做一个描述,然后推导2D CNN网络的卷积层和子采样层的BP权值更新方法。在推导过程中,我们更强调实现的效率,所以会给出一些Matlab代码。最后,我们转向讨论如何自动地学习组合前一层的特征maps,特别地,我们还学习特征maps的稀疏组合。

     

    二、全连接的反向传播算法

            典型的CNN中,开始几层都是卷积和下采样的交替,然后在最后一些层(靠近输出层的),都是全连接的一维网络。这时候我们已经将所有两维2D的特征maps转化为全连接的一维网络的输入。这样,当你准备好将最终的2D特征maps输入到1D网络中时,一个非常方便的方法就是把所有输出的特征maps连接成一个长的输入向量。然后我们回到BP算法的讨论。(更详细的基础推导可以参考UFLDL反向传导算法)。

    2.1Feedforward Pass前向传播

            在下面的推导中,我们采用平方误差代价函数。我们讨论的是多类问题,共c类,共N个训练样本。

            这里表示第n个样本对应的标签的第k维。表示第n个样本对应的网络输出的第k个输出。对于多类问题,输出一般组织为one-of-c的形式,也就是只有该输入对应的类的输出节点输出为正,其他类的位或者节点为0或者负数,这个取决于你输出层的激活函数。sigmoid就是0tanh就是-1.

            因为在全部训练集上的误差只是每个训练样本的误差的总和,所以这里我们先考虑对于一个样本的BP。对于第n个样本的误差,表示为:

           传统的全连接神经网络中,我们需要根据BP规则计算代价函数E关于网络每一个权值的偏导数。我们用l来表示当前层,那么当前层的输出可以表示为:

     

           输出激活函数f(.)可以有很多种,一般是sigmoid函数或者双曲线正切函数。sigmoid将输出压缩到[0, 1],所以最后的输出平均值一般趋于0。所以如果将我们的训练数据归一化为零均值和方差为1,可以在梯度下降的过程中增加收敛性。对于归一化的数据集来说,双曲线正切函数也是不错的选择。

     

    2.2Backpropagation Pass反向传播

            反向传播回来的误差可以看做是每个神经元的基的灵敏度sensitivities(灵敏度的意思就是我们的基b变化多少,误差会变化多少,也就是误差对基的变化率,也就是导数了),定义如下:(第二个等号是根据求导的链式法则得到的)

            因为∂u/∂b=1,所以∂E/∂b=∂E/∂u=δ,也就是说bias基的灵敏度∂E/∂b=δ和误差E对一个节点全部输入u的导数∂E/∂u是相等的。这个导数就是让高层误差反向传播到底层的神来之笔。反向传播就是用下面这条关系式:(下面这条式子表达的就是第l层的灵敏度,就是)

       公式(1

            这里的表示每个元素相乘。输出层的神经元的灵敏度是不一样的:

            最后,对每个神经元运用delta(即δ)规则进行权值更新。具体来说就是,对一个给定的神经元,得到它的输入,然后用这个神经元的delta(即δ)来进行缩放。用向量的形式表述就是,对于第l层,误差对于该层每一个权值(组合为矩阵)的导数是该层的输入(等于上一层的输出)与该层的灵敏度(该层每个神经元的δ组合成一个向量的形式)的叉乘。然后得到的偏导数乘以一个负学习率就是该层的神经元的权值的更新了:

      公式(2

            对于bias基的更新表达式差不多。实际上,对于每一个权值(W)ij都有一个特定的学习率ηIj

     

    三、Convolutional Neural Networks卷积神经网络

    3.1Convolution Layers卷积层

            我们现在关注网络中卷积层的BP更新。在一个卷积层,上一层的特征maps被一个可学习的卷积核进行卷积,然后通过一个激活函数,就可以得到输出特征map。每一个输出map可能是组合卷积多个输入maps的值:

           这里Mj表示选择的输入maps的集合,那么到底选择哪些输入maps呢?有选择一对的或者三个的。但下面我们会讨论如何去自动选择需要组合的特征maps。每一个输出map会给一个额外的偏置b,但是对于一个特定的输出map,卷积每个输入maps的卷积核是不一样的。也就是说,如果输出特征map j和输出特征map k都是从输入map i中卷积求和得到,那么对应的卷积核是不一样的。

    3.1.1Computing the Gradients梯度计算

            我们假定每个卷积层l都会接一个下采样层l+1。对于BP来说,根据上文我们知道,要想求得层l的每个神经元对应的权值的权值更新,就需要先求层l的每一个神经节点的灵敏度δ(也就是权值更新的公式(2))。为了求这个灵敏度我们就需要先对下一层的节点(连接到当前层l的感兴趣节点的第l+1层的节点)的灵敏度求和(得到δl+1),然后乘以这些连接对应的权值(连接第l层感兴趣节点和第l+1层节点的权值)W。再乘以当前层l的该神经元节点的输入u的激活函数f的导数值(也就是那个灵敏度反向传播的公式(1)的δl的求解),这样就可以得到当前层l每个神经节点对应的灵敏度δl了。

          然而,因为下采样的存在,采样层的一个像素(神经元节点)对应的灵敏度δ对应于卷积层(上一层)的输出map的一块像素(采样窗口大小)。因此,层l中的一个map的每个节点只与l+1层中相应map的一个节点连接。

         为了有效计算层l的灵敏度,我们需要上采样upsample这个下采样downsample层对应的灵敏度map(特征map中每个像素对应一个灵敏度,所以也组成一个map),这样才使得这个灵敏度map大小与卷积层的map大小一致,然后再将层lmap的激活值的偏导数与从第l+1层的上采样得到的灵敏度map逐元素相乘(也就是公式(1)。

           在下采样层map的权值都取一个相同值β,而且是一个常数。所以我们只需要将上一个步骤得到的结果乘以一个β就可以完成第l层灵敏度δ的计算。

           我们可以对卷积层中每一个特征map j重复相同的计算过程。但很明显需要匹配相应的子采样层的map(参考公式(1))

           up(.)表示一个上采样操作。如果下采样的采样因子是n的话,它简单的将每个像素水平和垂直方向上拷贝n次。这样就可以恢复原来的大小了。实际上,这个函数可以用Kronecker乘积来实现:

           好,到这里,对于一个给定的map,我们就可以计算得到其灵敏度map了。然后我们就可以通过简单的对层l中的灵敏度map中所有节点进行求和快速的计算bias基的梯度了:

       公式(3

           最后,对卷积核的权值的梯度就可以用BP算法来计算了(公式(2)。另外,很多连接的权值是共享的,因此,对于一个给定的权值,我们需要对所有与该权值有联系(权值共享的连接)的连接对该点求梯度,然后对这些梯度进行求和,就像上面对bias基的梯度计算一样:

           这里,中的在卷积的时候与逐元素相乘的patch,输出卷积map(u,v)位置的值是由上一层的(u, v)位置的patch与卷积核k_ij逐元素相乘的结果。

          咋一看,好像我们需要煞费苦心地记住输出map(和对应的灵敏度map)每个像素对应于输入map的哪个patch。但实际上,在Matlab中,可以通过一个代码就实现。对于上面的公式,可以用Matlab的卷积函数来实现:

           我们先对delta灵敏度map进行旋转,这样就可以进行互相关计算,而不是卷积(在卷积的数学定义中,特征矩阵(卷积核)在传递给conv2时需要先翻转(flipped)一下。也就是颠倒下特征矩阵的行和列)。然后把输出反旋转回来,这样我们在前向传播进行卷积的时候,卷积核才是我们想要的方向。

     

    3.2Sub-sampling Layers子采样层

            对于子采样层来说,有N个输入maps,就有N个输出maps,只是每个输出map都变小了。

           down(.)表示一个下采样函数。典型的操作一般是对输入图像的不同nxn的块的所有像素进行求和。这样输出图像在两个维度上都缩小了n倍。每个输出map都对应一个属于自己的乘性偏置β和一个加性偏置b

     

    3.2.1Computing the Gradients梯度计算

            这里最困难的是计算灵敏度map。一旦我们得到这个了,那我们唯一需要更新的偏置参数βb就可以轻而易举了(公式(3)。如果下一个卷积层与这个子采样层是全连接的,那么就可以通过BP来计算子采样层的灵敏度maps

            我们需要计算卷积核的梯度,所以我们必须找到输入map中哪个patch对应输出map的哪个像素。这里,就是必须找到当前层的灵敏度map中哪个patch对应与下一层的灵敏度map的给定像素,这样才可以利用公式(1那样的δ递推,也就是灵敏度反向传播回来。另外,需要乘以输入patch与输出像素之间连接的权值,这个权值实际上就是卷积核的权值(已旋转的)。

          在这之前,我们需要先将核旋转一下,让卷积函数可以实施互相关计算。另外,我们需要对卷积边界进行处理,但在Matlab里面,就比较容易处理。Matlab中全卷积会对缺少的输入像素补0

          到这里,我们就可以对bβ计算梯度了。首先,加性基b的计算和上面卷积层的一样,对灵敏度map中所有元素加起来就可以了:

          而对于乘性偏置β,因为涉及到了在前向传播过程中下采样map的计算,所以我们最好在前向的过程中保存好这些maps,这样在反向的计算中就不用重新计算了。我们定义:

    这样,对β的梯度就可以用下面的方式计算:

     

    3.3Learning Combinations of Feature Maps学习特征map的组合

            大部分时候,通过卷积多个输入maps,然后再对这些卷积值求和得到一个输出map,这样的效果往往是比较好的。在一些文献中,一般是人工选择哪些输入maps去组合得到一个输出map。但我们这里尝试去让CNN在训练的过程中学习这些组合,也就是让网络自己学习挑选哪些输入maps来计算得到输出map才是最好的。我们用αij表示在得到第j个输出map的其中第i个输入map的权值或者贡献。这样,第j个输出map可以表示为:

            需要满足约束:

            这些对变量αij的约束可以通过将变量αij表示为一个组无约束的隐含权值cijsoftmax函数来加强。(因为softmax的因变量是自变量的指数函数,他们的变化率会不同)。

            因为对于一个固定的j来说,每组权值cij都是和其他组的权值独立的,所以为了方面描述,我们把下标j去掉,只考虑一个map的更新,其他map的更新是一样的过程,只是map的索引j不同而已。

            Softmax函数的导数表示为:

           这里的δKronecker delta。对于误差对于第l层变量αi的导数为:

            最后就可以通过链式规则去求得代价函数关于权值ci的偏导数了:

     

    3.3.1Enforcing Sparse Combinations加强稀疏性组合

            为了限制αi是稀疏的,也就是限制一个输出map只与某些而不是全部的输入maps相连。我们在整体代价函数里增加稀疏约束项Ω(α)。对于单个样本,重写代价函数为:

    然后寻找这个规则化约束项对权值ci求导的贡献。规则化项Ω(α)αi求导是:

            然后,通过链式法则,对ci的求导是:

            所以,权值ci最后的梯度是:

     3.4Making it Fast with MATLAB

           CNN的训练主要是在卷积层和子采样层的交互上,其主要的计算瓶颈是:

    1)前向传播过程:下采样每个卷积层的maps

    2)反向传播过程:上采样高层子采样层的灵敏度map,以匹配底层的卷积层输出maps的大小;

    3sigmoid的运用和求导。

            对于第一和第二个问题,我们考虑的是如何用Matlab内置的图像处理函数去实现上采样和下采样的操作。对于上采样,imresize函数可以搞定,但需要很大的开销。一个比较快速的版本是使用Kronecker乘积函数kron。通过一个全一矩阵ones来和我们需要上采样的矩阵进行Kronecker乘积,就可以实现上采样的效果。对于前向传播过程中的下采样,imresize并没有提供在缩小图像的过程中还计算nxn块内像素的和的功能,所以没法用。一个比较好和快速的方法是用一个全一的卷积核来卷积图像,然后简单的通过标准的索引方法来采样最后卷积结果。例如,如果下采样的域是2x2的,那么我们可以用2x2的元素全是1的卷积核来卷积图像。然后再卷积后的图像中,我们每个2个点采集一次数据,y=x(1:2:end,1:2:end),这样就可以得到了两倍下采样,同时执行求和的效果。

            对于第三个问题,实际上有些人以为Matlab中对sigmoid函数进行inline的定义会更快,其实不然,MatlabC/C++等等语言不一样,Matlabinline反而比普通的函数定义更非时间。所以,我们可以直接在代码中使用计算sigmoid函数及其导数的真实代码。

     

    展开全文
  • MATLAB颜色图中,小于某个值的所有点设为白色 图1 原始图 如图1所示MATLAB输出的二维颜色图,这个图使用pcolor(x,y,Er)指令产生,x,y是对应的坐标,x=0:0.5:50;y=0:05:50;Er(101*101矩阵)表示(x,y)点计算出来的...

    1.问题提出

    如图1所示MATLAB输出的二维颜色图,这个图使用pcolor(x,y,Er)指令产生,x,y是对应的坐标,x=0:0.5:50;y=0:05:50;Er(101*101矩阵)表示(x,y)点计算出来的误差值。现在如果想把Er中小于0.05的所有点设置显示为白色,该如何做?
    图1

    图1 原始图

    2.MATLAB工程图形处理方式

    在解决这个问题之前,先来说一下MATLAB对彩色图形的处理方式。MATLAB中的图形常见的有以下几类:

    1. Truecolor:数码相机的格式,广泛用于计算机图形。
    2. Indexed 和 scaled indexed;经常用来显示科学或者工程数据,使用的不同的颜色深浅可以代表不同的数据大小。
    3. Grayscale:经常用在图像处理和图像分析算法中;
    4. Binary:经常用做为一个封装来表示图形的分割结果或者是感兴趣的区域。

    那么图1这种图形就属于第2种,使用颜色比例代表不同的数据大小。在上述图形中,Er中不同的值,是如何对应到不同的颜色的呢?使用pcolor画图时,需要介绍两个关键的数据,一个是colormap,一个是Cdata。

    3.什么是colormap

    先来看colormap。当用figure产生一个图窗时,里面就包含了一个预先设置的颜色数据矩阵,这个颜色数据矩阵的尺寸一般是64*3,每一行为介于0.0到1.0之间的值组成的RGB颜色数据。可以用MAPX=colormap,获取当前figure的颜色数据矩阵,如图2所示。实际上MATLAB预先提供了一些预定义的颜色数据矩阵,这些颜色数据矩阵的名称和对应的色阶如图3所示。

    图2

    图2,颜色数据矩阵

    图3

    图3 MATLAB预定义的颜色矩阵名称与色阶

    通常,MATLAB默认使用的颜色矩阵就是jet,大家看图1中的颜色条和jet色阶是不是一样的?都是从深蓝开始到深红结束。这些预定义的颜色矩阵都是64行(即64种颜色)。可以使用jet(n)指令,重新设置jet的行数,如jet(256),这样就把jet颜色扩充到256种颜色。需要注意的是,虽然扩展到了256种颜色,但是还是从深蓝开始到深红结束。只不过原来这两种颜色中间划分了64种颜色,而现在划分了256种颜色,划分的颜色更细腻而已。

    我们可以把新的颜色矩阵应用到figure中,例如:

    MAP=jet(256);
    colormap(MAP);
    

    这样,figure中使用的颜色矩阵就是被扩充到256色的jet颜色,而不是原来64色的jet颜色。

    4.什么是CData

    下面再来看CData数据。我们知道,图1中的颜色代表了Er的数值大小,即Er数值越大,颜色越红,Er数值越小,颜色越蓝。那么这种数据大小是如何对应到不同的颜色中的呢?实际上MATLAB是使用Er的值作为索引值,在colormap颜色矩阵中进行查表,从而知道该点对应的颜色。在MATLAB中,这种索引机制有两种,一种是直接索引direct,另一种是scaled,即按比例索引。到底用哪种索引方式,由当前图形句柄的’CDatamapping’属性决定。

    h=pcolor(x,y,Er);
    get(h,’ CDatamapping’);
    

    输出为scaled,这表明当前使用的索引方式是按比例索引。

    那么CData数据是什么呢?实际上当我们用h=pcolor(x,y,Er)画图时,真正使用的索引数据并不是Er,而是CData数据。CData是一个矩阵,大小和Er相同,在调用pcolor时,由Er按比例计算得到(即Er中数值大的对应的CData中的值就大,Er中数值小的对应的CData中的值就小)。一般而言,初始得到的CData的值和Er的值是相同的。可以使用Cdatat=get(h,‘CData’),获取当前的CData数据,如图4所示。

    图4

    图4 CData数据

    5.Scaled方式

    索引方式使用scaled方式时,CData的最小值对应colormap中颜色矩阵的第1行的颜色,CData的最大值对应colormap中的最后1行颜色,CData的中间值,就按比例关系,计算出对应的colormap中的行数,即该点的颜色。例如,假设colormap中共有64行颜色,CData的最大值为cmax(对应第64行),最小值为cmin(对应第1行),若CData有一个数据,值为cind,该值在colormap对应的行数为ind。那么按照比例原则,可列出如下公式:(cmax-cind)/(cind-cmin)=(64-ind)/(ind-1)。如果令k=(cmax-cind)/(cind-cmin),那么可计算出ind=(64+k)/(k+1),因为是行数,计算出来的值取整数。

    这样,根据CData中的数值大小,计算出在colormap中对应的行数,CData的值越小,对应的行号越小,CData的值越大,对应的行号越大。而CData是由Er按比例得到的,因此就实现了Er中不同的值,对应的颜色不同,即图1中的颜色图。

    另外,scaled比例的最大值和最小值一般就是CData中的最大值和最小值,但这个值也可以由axes的Clim属性控制。可以用cm=caxis获取当前scaled比例的最大值和最小值。cm(1)是最小值,cm(2)是最大值。可以使用caxis([newmin,newmax])来设置比例的最小值和最大值。此时CData中所有比newmin小的值,都对应colormap的第一行颜色,CData中所有比newmax大的值,都对应colormap的最后一行颜色。

    如果索引方式是direct,那么就用CData的真实值进行索引颜色,如Cdata的值是1.2,那对应的颜色就是colormap中的第1行,如CData的值是1000.4(超过了colormap的最大行数),那么就对应colormap中的最后一行。

    事实上,在获取CData的数据之后,我们还可以对获取的CData数据进行调整,然后把调整之后的CData数据重新应用到图形中。使用如下指令使用新的CData数据:set(h,‘CData’,Cdatat)。Cdatat是新调整后的CData数据。

    做一个梳理:
    当使用h=pcolor(x,y,Er)画图时,完成了下列工作:
    (1)根据Er,按比例计算CData矩阵;
    (2)CData的最小值对应colormap中的第1行颜色,CData的最大值对应colormap中最后一行颜色,CData的中间值按照比例关系进行计算对应的颜色。
    (3)索引出来每个点的颜色之后,使用该颜色对该点着色,得到图1。

    那么回过头来,看当初的问题。对于图1所示的图,如何把Er中小于0.05的数值显示为白色呢?一种思路是这样:
    (A)获取当前图形的colormap颜色数据;
    (B)在颜色数据的最前面插入一行白色[1,1,1];
    (C)获取当前CData值;
    (D)把所有的Er<=0.05对应的CData值设为CData的最小值,这样这些点就对应到了colormap中的第一行颜色(白色)
    (E)把修改后的CData数据应用到图形中。

    6.pcolor函数不能胜任

    但是这种方法会带来一些问题。如果Er中的数据大小分布比较均匀,用这种方法或许是可行的。但是如果Er中有突变数据,用这种方法就有可能带来一些问题。比如Er中有个别点的值非常大,Er中大于0.05的点(比如0.051),通过比例公式计算出来的行号有可能是1,而对应到了白色。显然按照要求,0.051对应的点不应显示白色,这样就造成了一些误差。有时甚至会带来严重的错误(比如0.051这种值很多,得到的结论就严重错误)。

    这个问题的根源是由于Er数值对应的颜色是按照比例关系计算得到的,而colormap中的颜色数量不够大(一般小于256,再多的行数意义也不大),但Er的尺寸可能非常大(导致Er中的每一点的值和colormap中的颜色不是一一对应的),而且里面的数据有可能有突变。对于Er中稍微大于0.05的值,按照比例关系索引出来的颜色有可能是白色,从而造成误差,甚至严重错误。

    7.patch登场

    如何解决这个问题呢?要改用另一个函数画图,使用patch函数画图。
    patch函数的一种用法是patch(X,Y,C);用来构建一个或者多个可填充的多边形,其使用X和Y作为每个点的坐标值,patch将会按顺序连接每个点。如果要得到一个多边形,将X和Y设置为向量;如果要得到多个多边形,将X和Y设置为矩阵,每一列对应一个多边形。C决定多边形的颜色,可以是系统认定的字符,也可以是一个数值,也可以是RGB向量。

    对于图1中的问题,实际上是100100个小方块,我们给这100100个小方块着色,其中Er值小于0.05对应的小方块着色为白色,大于0.05的方块按照比例关系在colormap中索引相应的颜色来着色。

    为了完成这个工作,需要做两件事:
    (1)构造这100*100个小方块的顶点数据矩阵X和Y;
    (2)设置每个小方块相应的颜色数据矩阵C。

    先看第一个。根据patch帮助手册,X和Y的列数对应整幅图中多边形的个数,X和Y的行数,对应每个多边形的顶点数。图1中有10000个小方块,每个小方块有4个顶点,把每个小方块4个顶点的x坐标存放在X每一列中,所以X矩阵尺寸为4*10000;同样,每个小方块4个顶点的y坐标存放在Y每一列中,如果已知每个小方块左下角顶点的坐标为(x0,y0),那么按照逆时针方向,其余三个顶点坐标依次是(x0+dx,y0)、(x0+dx,y0+dy)、(x0,y0+dy),如图5所示。已知x坐标和y坐标的范围及步长,就能把X和Y构造出来。
    图5

    图5 小方块四个顶点坐标

    再看如何确定C。这里摘抄一段patch的帮助。

    When X and Y are matrices, if C is a 1xn, where n is the number of columns in X and Y, then each face j=1:n is flat colored by the colormap index C(j). Note the special case of a 1x3 C is always assumed to be an RGB triplet ColorSpec and specifies the same flat color for each face. If C is a matrix the same size as X and Y, then it specifies the colors at the vertices as colormap indices and bilinear interpolation is used to color the faces. If C is 1xnx3, where n is the number of columns of X and Y, then each face j is flat colored by the RGB triplet C(1,j,:). If C is mxnx3, where X and Y are mxn, then each vertex (X(i,j),Y(i,j)) is colored by the RGB triplet C(i,j,:) and the face is colored using interpolation.

    关注一下黄色底纹部分:如果C是1xnx3矩阵(这里n是X和Y的列数,也就是小方块数),那么第j个小方块的颜色就被C(1,j,:)对应的RGB颜色着色(以flat方式,整个小方块着单色,如果是以interpolation方式,那么这个小方块以插值方式着色,小方块内是渐变色。显然这里我们要用flat方式)。

    我们可以定义一个C(1,n,3)矩阵,其中n是小方块数(10000),如果该方块对应的Er值小于0.05,那么就对C(1,j,:)直接赋值[1,1,1](白色),如果Er大于0.05,则按照前文所述的scaled方式,计算该点对应的colormap中的行数,把colormap中该行的颜色赋值给C(1,j,:)。这种方法就不会存在pcolor中的误差,是精确的显示结果。最终显示的结果如图6所示。

    在这里插入图片描述

    图6 图1中小于0.05的点着色白色

    附带MATLAB程序如下:

    clc;
    clear;
    load eremback.mat;%载入Er数据
    err_level=0.05;%误差门限
    xt=0:0.5:50;%生成x坐标
    yt=0:0.5:50;%生成y坐标
    n=length(xt);
    X=zeros(4,(n-1)*(n-1));
    Y=zeros(4,(n-1)*(n-1));
    C=zeros(1,(n-1)*(n-1),3);
    emax=max(max(Er));%Er的最大值
    emin=min(min(Er));%Er的最小值
    MAPX=jet((n-1)*(n-1));%对jet颜色矩阵扩展。这一步要不要关系不大
    lenMAPX=length(MAPX);
    for i=1:n-1
        for j=1:n-1
            X(1,100*(i-1)+j)=xt(i);
            X(2,100*(i-1)+j)=xt(i)+0.5;%构造X位置矩阵
            X(3,100*(i-1)+j)=xt(i)+0.5;
            X(4,100*(i-1)+j)=xt(i);
            
            Y(1,100*(i-1)+j)=yt(j);
            Y(2,100*(i-1)+j)=yt(j);%构造Y位置矩阵
            Y(3,100*(i-1)+j)=yt(j)+0.5;
            Y(4,100*(i-1)+j)=yt(j)+0.5;
            
            if Er(i,j)<=err_level;
                C(1,100*(i-1)+j,:)=[1,1,1];%如果该位置对应的Er值小于err_level,着白色
            else%Er大于0.05的值,从colormap中索引确定颜色
                kk=(emax-Er(i,j))/(Er(i,j)-emin);
                ind=round((lenMAPX+kk)/(kk+1));
                C(1,100*(i-1)+j,:)=MAPX(ind,:);%把colormap中的颜色赋值给C。
            end
        end
    end
    h=patch(X,Y,C);
    set(h,'edgecolor','none');
    xlabel('x/m');
    ylabel('y/m');
    title('Er');
    colorbar;
    caxis([emin,emax]);
    
    展开全文
  • 在排序向量“x”中搜索“v”并找到索引和关于向量 x ... 输入: x:数值向量, x 应该已经按升序排序(例如 2,7,20,...120) v:要在 x 中搜索的数值输出: i:v 相对于 x 的索引。 cv:x 中与 v 相等或最接近的
  • matlab显示曲线图中某个点的坐标

    万次阅读 2020-11-27 16:38:44
    在画好一幅曲线图后,有时候我们需要找出某个坐标的点,可以通过打开数据游标功能,此时鼠标点击的位置就会显示出坐标值,如下图: 此时保存下来的图也会带有显示的坐标。 如果需要显示多个坐标的,按住...
  • matlab取目标轮廓坐标

    2018-06-03 17:39:33
    化轮廓图像,经过8邻域扫描,得到顺时针方向各个轮廓点的坐标。部分参考网络,给出了解释与理解。
  • %将Mydata矩阵输出到指定的txt文件中 fid = fopen('data.txt','wt'); [m n] = size(Mydata); for i=1:m, for j=1:n, if j==i, %对角线上的元素置零 fprintf(fid,'%g, ',0); elseif j==n, %输出到行尾,换行
  • matlab输出图片

    2021-09-15 21:08:56
    for循环中,保存图片,并命名不同... fileName = ['特征1_幅值标准差/',titleName,'.png']; print(fileName,figure1,'-dpng','-r600'); 特征1_幅值标准差/是要保存的文件夹名字,titleName是图片要保存的名字 ...
  • Matlab 取结构体中某一特定字段的所有
  • matlab中变量的保存到excel中

    千次阅读 2018-10-07 15:11:20
    xlswrite('E:\我院\2018年秋各学科\计算机体系结构\chapter2_4_x.xlsx', x) xlswrite('E:\我院\2018年秋各学科\计算机体系结构\chapter2_4_data.xlsx', cost_of_die) 参考:...
  • 最近处理模式输出,总是要在B矩阵中寻找最接近A矩阵的最小值的数据和位置 minA=min(min(A)); 这里二维矩阵A的最小值就是minA [row,column]=find(A==min(min(A))); [row1,column1]对应minA在A中的位置。 ...
  • 查找矩阵中唯一元素为7的位置 注:元素的线性索引就是以矩阵的列为基准对元素进行标号,第一列排完,排第二列,具体说明如下图所示: % 开始 % 定义矩阵 A = [1 8 3;4 5 8;7 8 9]; % 查找元素7的线性索引 ind =...
  • 用于检查给定文件是否已存在的函数。 如果找到匹配的文件,它会创建一个可用于避免覆盖的新文件名。 新文件名是通过在文件名后附加一个像“_001”这样的序列来创建... 如果需要,枚举的起始和使用的位数是可调整的。
  • MATLAB教程(1) MATLAB 基础知识

    万次阅读 多人点赞 2017-10-26 20:57:32
    去年看过一点点MATLAB,很久不用,遗忘惊人。为了加深自己的印象,扎实基础,现将官网上的基础教程做简单的翻译。 首先,以下从九个部分简单介绍基础入门知识。第一部分:MATLAB显示桌面的基本布局...
  • 素数返回小于或等于某个值的素数的整个列表。因此,您可能会调用质数,仅表示您生成的质数太少而无法获得所需的特定质数。更糟糕的是,假设您想找到P(1e8)?仅生成列表的最后一个元素,生成100,000,000个素数的...
  • 您的算法非常慢,因为如果有(…)必须在第一次...在时间和记忆中是线性的一种解决方案可能是首先计算从该点到最后的最大的向量,其可以在一次向后传递中计算(你可以称之为累积最大,即cannot be vectorized).之后...
  • 当数据文件比较大时,靠直接地观察不能轻松地得到想要的结果,但matlab提供了很多简单的函数来对文件进行操作  文件操作首先要使用fopen或者load打开或者加载文件 ,下面是我最初使用的代码(带注释): data=...
  • 1. 先将那个图片显示出来 img_a = imread('test.tiff'); %写文件名字,文件放在matlab可检索的目录下 imshow(img_a); ...3. 在目标图像上点一下,显示该点的RGB,三者相等表...
  • MATLAB 绘图 横纵坐标制定内容显示(xlabel、ylabel)&amp;amp;amp;amp;MATLAB自定义输出图片...Matlab输出图片文件结果 功能代码 坐标斜体 xlabel('\itSin-Cos') 制定Y轴内容 set(haxe...
  • 一、对象句柄获取、 1、句柄、 2、创建对象时获取句柄、 3、函数获取句柄、 4、获取 / 设置 对象属性、 二、获取对象属性、 1、获取 线 对象属性、 2、获取 坐标轴 对象属性、
  • 每个序列有 4 列,其中: - 第一列:重复的- 第二列:序列的第一个的位置- 第三列:序列最后一个的位置- 第 4 列:序列的长度 [VALUES,INPOS,FIPOS,LEN] = findseq(...) OUT 作为单独的输出。 详情见...
  • MATLAB中求矩阵中最大所在的位置

    万次阅读 2017-09-11 17:24:11
    clear all a0=100*rand(1,9); a=reshape(a0,3,3); a_max=max(a(:)); [x,y]=find(a==max(a(:))); 求矩阵a中最大的元素,以及元素所在的位置
  • Matlab数值剔除

    千次阅读 2020-10-17 11:13:08
    因为输入输出都是多的,成矩阵排列,删除某一个意味着要删除相关的所人,而不仅仅是该本身。花了一天时间,终于在matlab时实现了。 数据处理 如下图所示,某一行数据中有一个超出了要求,哪么整行都要删除。...
  • I = imread('D:\456\Ga.BMP'); figure, imshow(I); I_gray = rgb2gray(I); figure, imshow(I_gray); % I_reverse = imcomplement(I_gray); % figure, imshow(I_reverse); I_reverse2 = 255 - I_gray;...
  • A = [1 2 1 1 2 3 4 1 1 3 3 3]; count = 0; for s = 1:length(A)-2 B = A(s:s+2);%获得连续三个数 if (length(find(B >= 2))==3) %判断三个元素是否大于等于2 count = count+1; %计数 ...count
  • 按指定的格式将变量的值输出到屏幕或指定文件  fid 为文件句柄,若缺省,则输出到屏幕  1 for standard output (the screen) or 2 for standard error. If FID is omitted, output goes to the screen. ...
  • 我在 Mathworks 上找到的所有教程都使用“helperCWTTimeFreqPlot”函数绘制了标度图,但对于较新版本的 Matlab,此函数被替换为没有输出参数的 CWT,没有输出参数的 CWT 根据采样频率自动调整频率和时间,这有时...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 9,651
精华内容 3,860
关键字:

matlab输出某个值

matlab 订阅