a=4;b=4;c=0.1;d=9;%常数初始
k=0.9;%对相似度的限制
k1=0.011;%对幅度值的限制
n=0;%类数
l=0;%记录输入的元素个数
h=0;%输入样本与聚类中心的幅度值
out=zeros(89,89);% 输出矩阵
w1=zeros(89,89);w2=zeros(89,89);%向上向下劝值初始
C=zeros(89,89);% 聚类中心
C1=zeros(89,89);C2=zeros(89,89);%中心迭代中间变量
for x=1:89% 输入原始数据
     l=l+1;
s=zeros(size(x));z=zeros(size(x));v=zeros(size(x));u=zeros(size(x));p=zeros(size(x));q=zeros(size(x));g=zeros(size(x));%各层初始

y=zeros(size(x));%F2层输出
t=zeros(size(x));%F2层输入

%使F1层达到稳定
if n==0 %F2层未被激活,第一个元素处理
    n=n+1;%类数加一个
    for m=1:length(x); %
        s1(m)=x(m);
        z1(m)=s1(m)/sqrt(max(cumsum(s1.^2)));
        v1(m)=f(z1(m),n);
        u1(m)=v1(m)/sqrt(max(cumsum(v1.^2)));
        p1(m)=u1(m);
        q1(m)=p1(m)/sqrt(max(cumsum(p1.^2)));
        while s1(m)~=s(m) & z1(m)~=z(m) & v1(m)~=v(m) & u1(m)~=u(m) & p1(m)~=p(m) & q1(m)~=q(m)
            s1(m)=x(m)+a*u(m);
            z1(m)=s1(m)/sqrt(max(cumsum(s1.^2)));
            v1(m)=f(z1(m),n)+b*f(q(m),n);
            u1(m)=v1(m)/sqrt(max(cumsum(v1.^2)));
            p1(m)=u1(m)+max(cumsum(g.*w2(:,m)));
            q1(m)=p1(m)/sqrt(max(cumsum(p1.^2)));
            s(m)=s1(m);z(m)=z1(m);v(m)=v1(m);u(m)=u1(m);p(m)=p1(m);q(m)=q1(m);
        end;
    end;
    out(1,1)=x;
    g(1)=d;
    
    for m=1:length(x)%修改权值
       C(1,m)=x(m);
       w2(1,m)=x(m)/sqrt(max(cumsum(x.^2)));
       w1(m,1)=w2(1,m);
    end;    
end;
%F2层已被激活    
if l>1   %输入元素个数大于1

    for i=1:n%类
         for m=1:length(x) %维
             s1(m)=x(m)+a*u(m);
             z1(m)=s1(m)/sqrt(max(cumsum(s1.^2)));
             v1(m)=f(z1(m),n)+b*f(q(m),n);
             u1(m)=v1(m)/sqrt(max(cumsum(s1.^2)));
             p1(m)=u1(m)+max(cumsum(g(:).*w2(:,m)));
             q1(m)=p1(m)/sqrt(max(cumsum(p1.^2)));
             while s1(m)~=s(m) & z1(m)~=z(m) & v1(m)~=v(m) & u1(m)~=u(m) & p1(m)~=p(m) & q1(m)~=q(m)
            
                 s1(m)=x(m)+a*u(m);
                 z1(m)=s1(m)/sqrt(max(cumsum(s1.^2)));
                 v1(m)=f(z1(m),n)+b*f(q(m),n);
                 u1(m)=v1(m)/sqrt(max(cumsum(v1.^2)));
                 p1(m)=u1(m)+max(cumsum(g.*w2(:,m)));
                 q1(m)=p1(m)/sqrt(max(cumsum(p1.^2)));
                 s(m)=s1(m);z(m)=z1(m);v(m)=v1(m);u(m)=u1(m);p(m)=p1(m);q(m)=q1(m);
             end;
         end;  
       r=(u+c*p)/(sqrt(max(cumsum(u.^2)))+sqrt(max(cumsum((c*p).^2))));
       if sqrt(max(cumsum(r.^2)))<k & i<n
           return;
       end;
        if sqrt(max(cumsum(r.^2)))<k & i==n
            n=n+1;
            out(l,n)=x;
            for k=1:n
                g(k)=0;
            end;
            g(n)=d;
             for k=1:length(x)
                 C(n,k)=x(k);
                 w2(n,k)=x(k)/sqrt(max(cumsum(x.^2)));
                 w1(k,n)=w2(n,k);
             end;
         end;
        if sqrt(max(cumsum(r.^2)))>=k & abs(sqrt(max(cumsum(x.^2)))-sqrt(max(cumsum(C(i,:).^2))))>k1 & i<n
            return;
        end;
        if sqrt(max(cumsum(r.^2)))>=k & abs(sqrt(max(cumsum(x.^2)))-sqrt(max(cumsum(C(i,:).^2))))>k1 & i==n
            n=n+1;
            out(l,n)=x;
            for k=1:n
                g(k)=0;
            end;
            g(n)=d;
             for k=1:length(x)
                 C(n,k)=x(k);
                 w2(n,k)=x(k)/sqrt(max(cumsum(x.^2)));
                 w1(k,n)=w2(n,k);
             end;
        end;
        if sqrt(max(cumsum(r.^2)))>=k & abs(sqrt(max(cumsum(x.^2)))-sqrt(max(cumsum(C(i,:).^2))))<=k1
            J=0;
            dJ=0;
            for j=1:l
                if out(j,i)~=0
                     J=J+sqrt(max(cumsum(out(j,i)-C(i,:)).^2))^2;
                     dJ=dJ-2*sqrt(max(cumsum(out(j,i)-C(i,:))));%对J求导
                end;
            end;
                C1(i,:)=C(i,:)+dJ;
                while C1(i,:)~=C(i,:)
                    for j=1:l
                        if out(j,i)~=0
                            J=J+sqrt(max(cumsum(out(j,i)-C(i,:)).^2))^2;
                            dJ=dJ-2*sqrt(max(cumsum(out(j,i)-C(i,:))));%对J求导
                        end;
                    end;
                    C1(i,:)=C1(i,:)+dJ;
                    C2(i,:)=C(i,:);
                    C(i,:)=C1(i,:);
                    C1(i,:)=C2(i,:);
                end;
                for k=1:length(x)
                    w2(i,k)=C(i,k)/sqrt(max(cumsum(C(i,:).^2)));
                    w1(k,i)=w2(i,k);
                end;
                for k=1:n
                    g(k)=0;
                end;
                g(i)=d;
                out(l,i)=x;
         end;
        end;
    end;
end;
out





function f=f(x,n)
d=9;h=1/(d*sqrt(n));
for i=1:length(x)
    if x(i)>=h
        f(i)=x(i);
    else
        f(i)=0;
    end;
end;