function [result,error,errorDim]=zjfc(s,opt)
% 對正交試驗進行方差分析，s是輸入矩陣，opt是空列參數向量，給出s中是空白列的列序號
%s=[  1  1   1   1  1 1 1 83.4;
%     1  1   1   2  2 2 2 84;
%     1  2   2   1  1 2 2 87.3;
%     1  2   2   2  2 1 1 84.8;
%     2  1   2   1  2 1 2 87.3;
%     2  1   2   2  1 2 1 88;
%     2  2   1   1  2 2 1 92.3;
%     2  2   1   2  1 1 2 90.4;
%];
%opt=[3,7];
% s的最后一列是各個正交組合的試驗測量值，前幾列是正交表
[m,n]=size(s);
 p=max(s(:,1));%取水平數
 q=n-1;% 取列數
 sum0=zeros(p,q);
 for i=1:q
    for k=1:m 
         for j=1:p
          if(s(k,i)==j)
              sum0(j,i)=sum0(j,i)+s(k,n);%求和
          end
        end
     end
 end
totalsum=sum(s(:,n));
ss=sum0.*sum0;
levelsum=m/p; %水平重復數
ss=sum(ss./levelsum)-totalsum^2/m; %每一列的S
ssError=sum(ss(opt));
for i=1:q
    f(i)=p-1;    %自由度
end
fError=sum(f(opt));  %誤差自由度
ssbar=ss./f;
Errorbar=ssError/fError;
index=find(ssbar<Errorbar);
index1=find(index==opt);
index(index==index(index1))=[];%剔除重復
ssErrorNew=ssError+sum(ss(index));  %并入誤差
fErrorNew=fError+sum(f(index));  %新誤差自由度
F=(ss./f)/(ssErrorNew./fErrorNew);   %F值
errorDim=[opt,index];

errorDim=sort(errorDim);%誤差列的序號
result=[ss',f',ssbar',F'];
error=[ssError,fError;ssErrorNew,fErrorNew];