SNN—脉冲神经网络学习算法之Tempotron代码示例

脉冲神经网络之Tempotron代码示例上一篇从原理的角度大致介绍了脉冲神经网络的神经元模型以及Tempotron监督学习方法，这一章记录了Tempotron的代码实现。这..

SNN学习

1885人浏览 · 2021-08-12 15:07:03

SNN学习 · 2021-08-12 15:07:03 发布

脉冲神经网络之Tempotron代码示例

上一篇从原理的角度大致介绍了脉冲神经网络的神经元模型以及Tempotron监督学习方法，这一章记录了Tempotron的代码实现。这份代码是使用matlab编写，用脉冲神经网络实现对26个字母分类，我们从细节去解释该代码：

function TempotronClassify()
% Tempotron: a neuron that learns spike timing-based decisions
% Rober Gutig 2006 Nature Neuroscience
clear; clc;
NumImages=26;
for i=1:NumImages
    ImageName=strcat('Icon16X16\Letter-',char('A'+i-1),'-black-icon');% 从icon16X16文件夹中读取所有图片
    ImageMatrix=imread(ImageName,'bmp');% 读取图片为灰度图，保存在矩阵中，取反
    ImageMatrix=~ImageMatrix;  % make the white pixel be 0, and black be 1;
    TrainPtns(:,i)=image2ptn(ImageMatrix);
end

上面代码片段就是从Icon16X16\文件夹中读取图片文件，图片从A到Z，这个代码的图片数据很简单，只有0和1，由黑白两色构成字母，读取的数据以16*16矩阵形式保存，之后调用image2ptn()方法，将矩阵转换。这个方法实现的就是对外界刺激进行编码，也就是将图片数据转换成脉冲序列的形式。打开该方法查看代码：

%% convert a image to a spike train saved in a vector
RandParts=1;

spikeTrain=zeros(32,1);
if RandParts1
loadR=1;
AR=A(😃;
if loadR0
R = randperm(size(A,1)size(A,2));
save(‘RandIndex’,‘R’);
else
load(‘RandIndex’,‘R’);
end
numRandNeus=32;
for i=1:numRandNeus
IndexR=R((1+(i-1)8):(8+(i-1)8));
Bits=AR(IndexR);
Bits=Bits’;
spikeTime=bi2de(Bits,‘left-msb’); % 二进制转十进制
if spikeTime==0
spikeTime=2^8; % put 0 to the end of the interval
end
spikeTrain(i)=spikeTime;
end
end这里的rand函数目的只是打乱1616矩阵中的行，我是不太理解这个=-=，AR是个2561的数组，是将A（图片矩阵1616）展平（flatten）的结果，Bits=AR(IndexR);这句简单的说就是从AR中随机抽取8个数字，构成bits，然后将bits转换成spikeTime（十进制数字），代表脉冲发放的时间。综上原先的图片1616矩阵被转换成328矩阵，8代表8位二进制编码，范围在0_{255之间，也就是说时间窗口大小为256毫秒，然后将8为二进制数转为十进制，代表脉冲发放的时间，该时间在（0}255）之间。总结起来，就是对于一张1616的黑白图片，转换成了一个321的数组，数组中有32个时间窗口，每个元素代表脉冲发送时间。


 
 
   
   
    
    
   
   
   
   
    
    
     
     TrainPtns=TrainPtns*1e-3;  % scale to ms
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     nAfferents = size(TrainPtns,1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     nPtns = NumImages;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %nOutputs = 1;    %%%%%%%%%%   1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     nOutputs = 5;

TrainPtns是32*26的矩阵，也就是将所有转换后的spikeTrain拼接一起。nAfferents代表了输入脉冲数量，这里就是32，由这32个输入确定一个输出（字母）。输出就定义为5，也就是5个二进制数表示一个字母的编号。


 
 
   
   
    
    
   
   
   
   
    
    
     
     loadData=0;% 是否载入已保存的模型
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     V_thr = 1; V_rest = 0;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     T = 256e-3;         % pattern duration ms
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     dt = 1e-3;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     tau_m = 20e-3; % tau_m = 15e-3;？？？
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     tau_s = tau_m/4;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % K(t?ti)=V0(exp[?(t?ti)/τ]–exp[?(t?ti)/τs])
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     aa = exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s);
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     V0 = 1/max(exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s));
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     lmd = 2e-2;%1e-2/V0;   % optimal performance lmd=3e-3*T/(tau_m*nAfferents*V0)  1e-4/V0;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     maxEpoch = 200;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     mu = 0.99;  % momentum factor

上面定义了阈值电压V_thr，这是按经验来取，这里取1，复位电压则取0.定义了时间窗口大小为256毫秒，单位时间是1毫秒，aa则是为了核函数K的计算简化定义。即预先定义公式：

K(t-ti)=V0(exp[-(t-ti)/τ]–exp[-(t-ti)/τs])

V0也被预先定义。最大训练次数为200次（maxEpoch），lmd类似于学习率，而mu则是动量，该代码用了动量优化模型。


 
 
   
   
    
    
   
   
   
   
    
    
     
     if loadData ==0 %初始化网络
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         weights = 1e-2*randn(nAfferents,nOutputs);  % 1e-3*randn(nAfferents,1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         save('weights0','weights');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     else
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         load('weights0','weights');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end

随机初始化权重，或者读取保存的权重文件。


 
 
   
   
    
    
   
   
   
   
    
    
     
     Class = de2bi(1:26,'left-msb'); Class=Class';
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     correctRate=zeros(1,maxEpoch);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     dw_Past=zeros(nAfferents,nPtns,nOutputs);  % momentum for accelerating learning.上一个权重的更新，用于动量计算
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     for epoch=1:maxEpoch    
    
    
   
   

   
   
    
    
   
   
   
   
    
        
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         Class_Tr = false(nOutputs,nPtns);  % actual outputs of training
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         for pp=1:nPtns                 
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             for neuron=1:nOutputs
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 Vmax=0; tmax=0;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 fired=false;        
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 Vm1=zeros(1,256); indx1= 1; % trace pattern 1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 for t=dt:dt:T
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     Vm = 0; 
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     if fired==false
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Tsyn=find(TrainPtns(:,pp)
     
     <=t+0.1*dt);    % no cut window
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     else
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     Tsyn=
     
     find(TrainPtns(:,pp)<=
     
     t_fire+0.1*dt); % 
     
     shut 
     
     down 
     
     inputs
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     if ~
     
     isempty(
     
     Tsyn)                    
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A1=
     
     TrainPtns(:,pp);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A2=
     
     A1(Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     K =
     
     V0*(exp(-(t-A2)/tau_m)-exp(-(t-A2)/tau_s)); % 
     
     the 
     
     kernel 
     
     value 
     
     for 
     
     each 
     
     fired 
     
     afferent
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A1=
     
     weights(:,neuron);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     firedWeights=
     
     A1(Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     Vm = 
     
     Vm + 
     
     firedWeights'*
     
     K ;
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     Vm = 
     
     Vm + 
     
     V_rest;
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     if 
     
     Vm>=V_thr && fired==false % fire
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         fired=true;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         t_fire=t;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Class_Tr(neuron,pp)=true;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     if Vm>Vmax
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Vmax=Vm; tmax=t;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     end
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     %if pp==1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Vm1(indx1)=Vm;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         indx1=indx1+1;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     %end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 end
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 %if pp==1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     figure(1); plot(dt:dt:T,Vm1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     title(strcat('Image ',char('A'+pp-1),'; neuron: ',num2str(neuron))); drawnow;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 %end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 if Vmax
     
     <=0
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     tmax=
     
     max(TrainPtns(:,pp));
    
    
   
   

   
   
    
    
   
   
   
   
    
                
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
                
    
    
   
   

   
   
    
    
   
   
   
   
    
                
     
     if 
     
     Class_Tr(
     
     neuron,
     
     pp)~=
     
     Class(neuron,pp) % 
     
     error
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     Tsyn=
     
     find(TrainPtns(:,pp)<=
     
     tmax+0.1*dt); 
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     if ~
     
     isempty(
     
     Tsyn)                    
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A1=
     
     TrainPtns(:,pp);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A2=
     
     A1(Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     K =
     
     V0*(exp(-(tmax-A2)/tau_m)-exp(-(tmax-A2)/tau_s)); % 
     
     the 
     
     kernel 
     
     value 
     
     for 
     
     each 
     
     fired 
     
     afferent
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A1=
     
     weights(:,neuron);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     dwPst=
     
     dw_Past(:,pp,neuron);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     if 
     
     fired==
     
     false    % 
     
     LTP
    
    
   
   

   
   
    
    
   
   
   
   
    
                            
     
     Dw=
     
     lmd*K;
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     else           % 
     
     LTD
    
    
   
   

   
   
    
    
   
   
   
   
    
                            
     
     Dw=
     
     -1.1*lmd*K;
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A1(
     
     Tsyn) = 
     
     A1(Tsyn) + 
     
     Dw + 
     
     mu*
     
     dwPst(
     
     Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     weights(
     
     :,
     
     neuron)=
     
     A1;
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     dwPst(
     
     Tsyn)=
     
     Dw;
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     dw_Past(
     
     :,
     
     pp,
     
     neuron) = 
     
     dwPst;
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     end                
    
    
   
   

   
   
    
    
   
   
   
   
    
                
     
     end            
    
    
   
   

   
   
    
    
   
   
   
   
    
                
    
    
   
   

   
   
    
    
   
   
   
   
    
            
     
     end  % 
     
     end 
     
     of 
     
     one 
     
     neuron 
     
     computation
    
    
   
   

   
   
    
    
   
   
   
   
    
            
    
    
   
   

   
   
    
    
   
   
   
   
    
       
     
     end % 
     
     end 
     
     of 
     
     one 
     
     image
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
        %
     
     CC=
     
     isequal(Class,Class_Tr);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
        %
     
     correctRate(
     
     epoch)=
     
     sum(Class==
     
     Class_Tr)/length(Class);
    
    
   
   

   
   
    
    
   
   
   
   
    
       
     
     CC = 
     
     bi2de(Class_Tr','
     
     left-msb');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     save('
     
     TrainedWt','
     
     weights');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     figure(
     
     2); 
     
     plot(
     
     1:maxEpoch,
     
     correctRate,'
     
     -b.');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end

class代表最后的类别，由5位二进制编码表示。最外层循环是最大训练次数，第二层nPtns代表了26个样本，最后一层循环nOutputs为5，代表对每个输出神经元进行一次计算。之后将Vmax，Tmax（代表最大膜电位出现时间）初始化为0，fired表示是否发出脉冲，这里初始化为false，Vm1记录了此输出神经元膜电位（在一个时间窗口内）的变化，之后进入循环，循环一个时间窗口长度，对一个时间窗口内的每个单位时间，若该神经元已经发送过脉冲，则shut down所有脉冲输入；若该神经元还未发出过脉冲（fired = false)，就执行Tsyn=find(TrainPtns(:,pp)<=t+0.1*dt);也就是在当前样本（pp）的脉冲序列中，找到比时间t+0.1*dt小的脉冲发送时间，也就是在t这个时间点之前有脉冲输入发生，此时A2就保存了该脉冲发送的时间（也就是ti），之后就可调用K（t - ti）来计算此时的核函数；下一步就是找到发送这个脉冲的输入神经元和该输出神经元之间连接的权重数值，用K*Weights就代表了该脉冲输入对膜电位的影响，最后将该贡献值累加到Vm上，得到此刻膜电位大小。

对整个时间窗口计算完后，将Vm加上复位电压，如果该值超过了阈值电压，则发放脉冲；若该值超过了已有的Vmax，则更新Vmax。之后进入判断语句：

if Class_Tr(neuron,pp)~=Class(neuron,pp) 神经元的输出和实际输出不相符时，分为两个情况：

（1）实际输出为0，但发放了脉冲。

此时应该抑制脉冲发放，于是权重应该减小：

Dw=-1.1*lmd*K;

（2）实际输出为1，但没有发放脉冲：

此时应该增大神经元的刺激，权重应该增加：

Dw=lmd*K;
最后保存此次循环的权重。于是单层的脉冲神经网路监督学习就这样完成了。下面附上完整代码：


 
 
   
   
    
    
   
   
   
   
    
    
     
     function TempotronClassify()
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % Tempotron: a neuron that learns spike timing-based decisions
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % Rober Gutig 2006 Nature Neuroscience
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     clear; clc;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     NumImages=26;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     for i=1:NumImages
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         ImageName=strcat('Icon16X16\Letter-',char('A'+i-1),'-black-icon');% 从icon16X16文件夹中读取所有图片
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         ImageMatrix=imread(ImageName,'bmp');% 读取图片为灰度图，保存在矩阵中，取反
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         ImageMatrix=~ImageMatrix;  % make the white pixel be 0, and black be 1;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         TrainPtns(:,i)=image2ptn(ImageMatrix);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     TrainPtns=TrainPtns*1e-3;  % scale to ms
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     nAfferents = size(TrainPtns,1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     nPtns = NumImages;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %nOutputs = 1;    %%%%%%%%%%   1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     nOutputs = 5;
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     loadData=0;% 是否载入已保存的模型
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     V_thr = 1; V_rest = 0;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     T = 256e-3;         % pattern duration ms
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     dt = 1e-3;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     tau_m = 20e-3; % tau_m = 15e-3;？？？
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     tau_s = tau_m/4;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % K(t?ti)=V0(exp[?(t?ti)/τ]–exp[?(t?ti)/τs])
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     aa = exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s);
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     V0 = 1/max(exp(-(0:dt:3*tau_m)/tau_m)-exp(-(0:dt:3*tau_m)/tau_s));
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     lmd = 2e-2;%1e-2/V0;   % optimal performance lmd=3e-3*T/(tau_m*nAfferents*V0)  1e-4/V0;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     maxEpoch = 200;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     mu = 0.99;  % momentum factor
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % generate patterns (each pattern consists one spik-e per afferent)
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     if loadData ==0 %初始化网络
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         weights = 1e-2*randn(nAfferents,nOutputs);  % 1e-3*randn(nAfferents,1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         save('weights0','weights');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     else
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         load('weights0','weights');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %Class = logical(eye(nOutputs));     % desired class label for each pattern
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %Class = false(1,26); Class(26)=true;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     Class = de2bi(1:26,'left-msb'); Class=Class';
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     correctRate=zeros(1,maxEpoch);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     dw_Past=zeros(nAfferents,nPtns,nOutputs);  % momentum for accelerating learning.上一个权重的更新，用于动量计算
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     for epoch=1:maxEpoch    
    
    
   
   

   
   
    
    
   
   
   
   
    
        
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         Class_Tr = false(nOutputs,nPtns);  % actual outputs of training
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         for pp=1:nPtns 
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
      %       Class_Tr = false(nOutputs,1);  % actual outputs of training
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             for neuron=1:nOutputs
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 Vmax=0; tmax=0;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 fired=false;        
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 Vm1=zeros(1,256); indx1= 1; % trace pattern 1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 for t=dt:dt:T
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     Vm = 0; 
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     if fired==false
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Tsyn=find(TrainPtns(:,pp)
     
     <=t+0.1*dt);    % no cut window
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     else
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     Tsyn=
     
     find(TrainPtns(:,pp)<=
     
     t_fire+0.1*dt); % 
     
     shut 
     
     down 
     
     inputs
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     if ~
     
     isempty(
     
     Tsyn)                    
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A1=
     
     TrainPtns(:,pp);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A2=
     
     A1(Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     K =
     
     V0*(exp(-(t-A2)/tau_m)-exp(-(t-A2)/tau_s)); % 
     
     the 
     
     kernel 
     
     value 
     
     for 
     
     each 
     
     fired 
     
     afferent
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     A1=
     
     weights(:,neuron);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     firedWeights=
     
     A1(Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
                        
     
     Vm = 
     
     Vm + 
     
     firedWeights'*
     
     K ;
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     Vm = 
     
     Vm + 
     
     V_rest;
    
    
   
   

   
   
    
    
   
   
   
   
    
                    
     
     if 
     
     Vm>=V_thr && fired==false % fire
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         fired=true;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         t_fire=t;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Class_Tr(neuron,pp)=true;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     if Vm>Vmax
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Vmax=Vm; tmax=t;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     end
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     %if pp==1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         Vm1(indx1)=Vm;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         indx1=indx1+1;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     %end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 end
    
    
   
   

   
   
    
    
   
   
   
   
    
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 %if pp==1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     figure(1); plot(dt:dt:T,Vm1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     title(strcat('Image ',char('A'+pp-1),'; neuron: ',num2str(neuron))); drawnow;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 %end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 if Vmax
     
     <=0
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     tmax=max(TrainPtns(:,pp));
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 if Class_Tr(neuron,pp)~=Class(neuron,pp) % error
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     Tsyn=find(TrainPtns(:,pp)<=tmax+0.1*dt); 
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     if ~isempty(Tsyn)                    
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         A1=TrainPtns(:,pp);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         A2=A1(Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         K =V0*(exp(-(tmax-A2)/tau_m)-exp(-(tmax-A2)/tau_s)); % the kernel value for each fired afferent
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         A1=weights(:,neuron);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         dwPst=dw_Past(:,pp,neuron);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         if fired==false    % LTP
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                             Dw=lmd*K;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         else           % LTD
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                             Dw=-1.1*lmd*K;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         A1(Tsyn) = A1(Tsyn) + Dw + mu*dwPst(Tsyn);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         weights(:,neuron)=A1;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         dwPst(Tsyn)=Dw;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                         dw_Past(:,pp,neuron) = dwPst;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                     end                
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 end            
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             end  % end of one neuron computation
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
        end % end of one image
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
        %CC=isequal(Class,Class_Tr);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
        %correctRate(epoch)=sum(Class==Class_Tr)/length(Class);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
        CC = bi2de(Class_Tr','left-msb');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     save('TrainedWt','weights');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     figure(2); plot(1:maxEpoch,correctRate,'-b.');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %%将图片编码为脉冲序列并保存在向量中
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     function spikeTrain=image2ptn(A)  
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %% convert a image to a spike train saved in a vector
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     RandParts=1;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % A1=A';
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % B=[A1(:);A(:)];
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % numPixels=length(B);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % numInputNeurons=numPixels/8; % 64 neurons
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % spikeTrain=zeros(numInputNeurons,1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % for i=1:numInputNeurons
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %     Bits=B((1+(i-1)*8):(8+(i-1)*8));
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %     Bits=Bits';
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %     spikeTime=bi2de(Bits,'left-msb');    
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %     if spikeTime==0
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %         spikeTime=2^8;  % put 0 to the end of the interval
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     %     spikeTrain(i)=spikeTime;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     % end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     spikeTrain=zeros(32,1);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     if RandParts==1
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         loadR=1;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         AR=A(:);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         if loadR==0
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             R = randperm(size(A,1)*size(A,2));
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             save('RandIndex','R');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         else
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             load('RandIndex','R');
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         numRandNeus=32;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         for i=1:numRandNeus
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             IndexR=R((1+(i-1)*8):(8+(i-1)*8));
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             Bits=AR(IndexR);
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             Bits=Bits';
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             spikeTime=bi2de(Bits,'left-msb');    % 二进制转十进制
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             if spikeTime==0
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
                 spikeTime=2^8;  % put 0 to the end of the interval
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
             spikeTrain(i)=spikeTime;
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
         end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end
    
    
   
   

   
   
    
    
   
   
   
   
    
    
     
     end