euclids.m
上传用户:xgxlgz
上传日期:2013-06-22
资源大小:6k
文件大小:19k
源码类别:

DVD

开发平台:

Matlab

  1. clear
  2. clc
  3. %此处是(72,64,9)RS码,本原多项式是x^7+x^6+1,生成多项式是g(x)=П (x-α^j),
  4.  %展开后的系数为[1,124,73,119,39,43,24,49,116],指数形式为:[ 0,86,106,49,28,58,124,113,36]
  5. L=10000;%编码个数,2进制表示
  6. m=7;
  7. N=2^m-1-55;
  8. K=N-8;
  9. msg=randint(L,1);%产生待编码的随机数
  10. for i=1:1:(length(msg)-1)/(m*K)  end,i;
  11. %若不够一个码组,后面补零
  12. msg=[msg;zeros((i+1)*m*K-length(msg),1)]; 
  13.   %指数转化成整数,a^0~a^126
  14.    exp_int=[1    2    4    8    16   32   64   65   67   71   79   95    127  63   126  61  122  53   106  21   42   84 ...     %22个数一行
  15.             105  19   38   76   89   115  39   78   93   123  55   110   29   58   116  41  82   101  11   22   44   88 ...
  16.             113  35   70   77   91   119  47   94   125  59   118  45    90   117  43   86  109  27   54   108  25   50 ...
  17.             100  9    18   36   72   81   99   7    14   28   56   112   33   66   69   75  87   111  31   62   124  57 ...
  18.             114  37   74   85   107  23   46   92   121  51   102  13    26   52   104  17  34   68   73   83   103  15 ...
  19.             30   60   120  49   98   5    10   20   40   80   97   3     6    12   24   48  96];
  20.    
  21.    %整数转化指数,1~127
  22.    int_exp=[0    1    121  2    115  122  73   3    67   116  40   123  99   74   109  4   103  68   23   117  19   41 ...
  23.             93   124  64   100  61   75   34   110  84   5    78   104  45   69   89   24  28   118  37   20    58  42 ... 
  24.             55   94   50   125  113  65   97   101  17   62   32   76   87   35   53   111 15   85   13   6     7   79 ...
  25.             8    105  80   46   9    70   106  90   81   25   47   29   10   119  71   38  107  21   91   59    82  43 ... 
  26.             26   56   48   95   30   51   11   126  120  114  72   66   39   98   108  102 22   18   92   63    60  33 ...
  27.             83   77   44   88   27   36   57   54   49   112  96   16   31   86   52   14  12];
  28.    
  29.     %将整数转化成m位的2进制
  30.     for i=1:N+55
  31.     a=exp_int(i);
  32.     for j=1:m
  33.         exp_int_2(i,m+1-j)=mod(a,2);
  34.         a=(a-exp_int_2(i,m+1-j))/2;
  35.     end,j;
  36.    end,i;
  37.     
  38. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  39. %编码
  40. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  41. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  42. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  43. %中间这一段是为了可以调用系统函数,为了测试方便
  44. for i=1:length(msg)/(m*K)
  45.     for j=1:K
  46.         %此处约定前面的数据为高位,即若收到d0d1d2d3,则认为d0为最高位
  47.         msg_1(i,j)=msg(K*m*i-K*m+m*j-m+1)*64+msg(K*m*i-K*m+m*j-m+2)*32+msg(K*m*i-K*m+m*j-m+3)*16+msg(K*m*i-K*m+m*j-m+4)*8+msg(K*m*i-K*m+m*j-m+5)*4+msg(K*m*i-K*m+m*j-m+6)*2+msg(K*m*i-K*m+m*j);
  48.     end,j;
  49. end,i;
  50. %将数据流转化到gf域,以便于rsenc函数的执行
  51. genpoly=rsgenpoly(N+55,K+55,193);
  52. msg_1=gf(msg_1,m,193);
  53. code=rsenc(msg_1,N,K,genpoly);
  54. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  55. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  56. %下面直接进行缩短编码
  57. for i=1:1:length(msg)/(m*K)
  58.     %编系统码字,前面K位直接赋值
  59.     reg0=[0,0,0,0,0,0,0]; reg1=[0,0,0,0,0,0,0]; reg2=[0,0,0,0,0,0,0]; reg3=[0,0,0,0,0,0,0];
  60.     reg4=[0,0,0,0,0,0,0]; reg5=[0,0,0,0,0,0,0]; reg6=[0,0,0,0,0,0,0]; reg7=[0,0,0,0,0,0,0];
  61.  for j=1:N
  62.      if(j>K) x=[0,0,0,0,0,0,0]; 
  63.      else x=xor(reg7,msg(i*m*K-m*K+m*j-m+1:i*m*K-m*K+m*j)');%计算反馈回路上的数
  64.      end%相当二选一器,计数器控制
  65.      
  66.      if(j>K)  my_code(N*m*i-N*m+m*j-m+1:N*m*i-N*m+m*j)=reg7;
  67.      else     my_code(N*m*i-N*m+m*j-m+1:N*m*i-N*m+m*j)=msg(i*m*K-m*K+m*j-m+1:i*m*K-m*K+m*j);
  68.      end      %相当一个二选一,由记数器控制
  69.  
  70.     %电路实现乘法结构    
  71.      b0(1)=mod(x(4)+x(5)+x(7),2);b0(2)=mod(x(1)+x(4)+x(6)+x(7),2);
  72.      b0(3)=mod(x(2)+x(5)+x(7),2);b0(4)=mod(x(1)+x(3)+x(6),2);
  73.      b0(5)=mod(x(1)+x(2)+x(4)+x(7),2);b0(6)=mod(x(2)+x(3)+x(5),2);
  74.      b0(7)=mod(x(3)+x(4)+x(6),2);% *a^36    [a4+a5+a7, a1+a4+a6+a7,    a2+a5+a7,    a1+a3+a6, a1+a2+a4+a7,    a2+a3+a5,    a3+a4+a6]
  75.     
  76.      b1(1)=mod(x(1)+x(6),2);b1(2)=mod(x(2)+x(6)+x(7),2);
  77.      b1(3)=mod(x(1)+x(3)+x(7),2);b1(4)=mod(x(2)+x(4),2);
  78.      b1(5)=mod(x(3)+x(5),2);b1(6)=mod(x(4)+x(6),2);
  79.      b1(7)=mod(x(5)+x(7),2);% *a^113   [a1+a6, a2+a6+a7, a1+a3+a7,    a2+a4,    a3+a5,    a4+a6,    a5+a7]
  80.      
  81.      b2(1)=x(5);b2(2)=mod(x(5)+x(6),2);
  82.      b2(3)=mod(x(6)+x(7),2);b2(4)=mod(x(1)+x(7),2);
  83.      b2(5)=x(2);b2(6)=x(3);
  84.      b2(7)=x(4);%*a^124   [a5, a5+a6, a6+a7, a1+a7,    a2,    a3,    a4]
  85.      
  86.      b3(1)=mod(x(2)+x(5)+x(6),2);b3(2)=mod(x(2)+x(3)+x(5)+x(7),2);
  87.      b3(3)=mod(x(1)+x(3)+x(4)+x(6),2);b3(4)=mod(x(1)+x(2)+x(4)+x(5)+x(7),2);
  88.      b3(5)=mod(x(2)+x(3)+x(5)+x(6),2);b3(6)=mod(x(3)+x(4)+x(6)+x(7),2);
  89.      b3(7)=mod(x(1)+x(4)+x(5)+x(7),2); %a^58    [a2+a5+a6, a2+a3+a5+a7,a1+a3+a4+a6, a1+a2+a4+a5+a7, a2+a3+a5+a6, a3+a4+a6+a7, a1+a4+a5+a7]
  90.      
  91.      b4(1)=mod(x(2)+x(4)+x(5)+x(6),2);b4(2)=mod(x(2)+x(3)+x(4)+x(7),2);
  92.      b4(3)=mod(x(1)+x(3)+x(4)+x(5),2);b4(4)=mod(x(1)+x(2)+x(4)+x(5)+x(6),2);
  93.      b4(5)=mod(x(1)+x(2)+x(3)+x(5)+x(6)+x(7),2);b4(6)=mod(x(2)+x(3)+x(4)+x(6)+x(7),2);
  94.      b4(7)=mod(x(1)+x(3)+x(4)+x(5)+x(7),2);  % a^28   [a2+a4+a5+a6,a2+a3+a4+a7,a1+a3+a4+a5,a1+a2+a4+a5+a6, a1+a2+a3+a5+a6+a7,a2+a3+a4+a6+a7,a1+a3+a4+a5+a7]
  95.       
  96.      b5(1)=mod(x(2)+x(4)+x(5)+x(7),2);b5(2)=mod(x(1)+x(2)+x(3)+x(4)+x(6)+x(7),2);
  97.      b5(3)=mod(x(2)+x(3)+x(4)+x(5)+x(7),2);b5(4)=mod(x(1)+x(3)+x(4)+x(5)+x(6),2);
  98.      b5(5)=mod(x(1)+x(2)+x(4)+x(5)+x(6)+x(7),2);b5(6)=mod(x(2)+x(3)+x(5)+x(6)+x(7),2);
  99.      b5(7)=mod(x(1)+x(3)+x(4)+x(6)+x(7),2);%a^49  [a2+a4+a5+a7, a1+a2+a3+a4+a6+a7,a2+a3+a4+a5+a7,a1+a3+a4+a5+a6, a1+a2+a4+a5+a6+a7,    a2+a3+a5+a6+a7,    a1+a3+a4+a6+a7]
  100.       
  101.      b6(1)=mod(x(1)+x(5)+x(6)+x(7),2);b6(2)=mod(x(1)+x(2)+x(5),2);
  102.      b6(3)=mod(x(1)+x(2)+x(3)+x(6),2);b6(4)=mod(x(1)+x(2)+x(3)+x(4)+x(7),2);
  103.      b6(5)=mod(x(2)+x(3)+x(4)+x(5),2);b6(6)=mod(x(3)+x(4)+x(5)+x(6),2);
  104.      b6(7)=mod(x(4)+x(5)+x(6)+x(7),2); %a^106     [a1+a5+a6+a7, a1+a2+a5, a1+a2+a3+a6, a1+a2+a3+a4+a7, a2+a3+a4+a5,  a3+a4+a5+a6,  a4+a5+a6+a7]
  105.      
  106.      b7(1)=mod(x(1)+x(2)+x(3)+x(5)+x(7),2);b7(2)=mod(x(1)+x(4)+x(5)+x(6)+x(7),2);
  107.      b7(3)=mod(x(2)+x(5)+x(6)+x(7),2);b7(4)=mod(x(1)+x(3)+x(6)+x(7),2);
  108.      b7(5)=mod(x(2)+x(4)+x(7),2);b7(6)=mod(x(1)+x(3)+x(5),2);
  109.      b7(7)=mod(x(1)+x(2)+x(4)+x(6),2);%a^86   [a1+a2+a3+a5+a7, a1+a4+a5+a6+a7, a2+a5+a6+a7, a1+a3+a6+a7,    a2+a4+a7,  a1+a3+a5,    a1+a2+a4+a6]
  110.  
  111.      reg7=xor(reg6,b7);reg6=xor(reg5,b6);reg5=xor(reg4,b5);reg4=xor(reg3,b4);
  112.      reg3=xor(reg2,b3);reg2=xor(reg1,b2);reg1=xor(reg0,b1);reg0=b0;
  113.  end,j;
  114. end,i;  
  115. for i=1:1:length(my_code)/m
  116.     my_code1(i)=my_code(m*i-m+1)*64+my_code(m*i-m+2)*32+my_code(m*i-m+3)*16+my_code(m*i-m+4)*8+my_code(m*i-m+5)*4+my_code(m*i-m+6)*2+my_code(m*i);
  117. end,i;
  118. %将码字转化成十进制,以便跟系统编码进行比较
  119. %
  120. %先将自己编写的码字按码组形式(72个码字构成一个码组)存储在数组中,该数组一行正好是一个码字
  121.  my_code3=reshape(my_code1,N,length(my_code1)/N);
  122.  my_code3=my_code3'
  123.  c=my_code3-code;
  124. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  125. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  126.     
  127. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  128. %加噪
  129. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  130. noise=rand(1,length(my_code))<0.007;  %产生随机噪声
  131. %noise=zeros(1,length(my_code));
  132. %noise(7)=1;
  133. for i=1:length(my_code)
  134.     my_code(i)=xor(my_code(i),noise(i));%加噪
  135. end,i;
  136. for i=1:(length(my_code)/m)
  137.     my_code2(i)=my_code(m*i-m+1)*64+my_code(m*i-m+2)*32+my_code(m*i-m+3)*16+my_code(m*i-m+4)*8+my_code(m*i-m+5)*4+my_code(m*i-m+6)*2+my_code(m*i-m+7);
  138. end,i;  
  139. %从开头到现在是获取加噪后的码流,以模拟信道
  140. %===============================================
  141.  for i=1:(length(noise)/m)
  142.     noise2(i)=noise(m*i-m+1)*64+noise(m*i-m+2)*32+noise(m*i-m+3)*16+noise(m*i-m+4)*8+noise(m*i-m+5)*4+noise(m*i-m+6)*2+noise(m*i-m+7);
  143.  end,i;  
  144.  noise2=reshape(noise2,N,length(my_code)/(m*N));
  145.  noise2=noise2';
  146.  sum_noise=noise2>0;
  147.  sum_noise=sum(sum_noise,2);   %计算加噪个数,也就是错码个数
  148.  my_code4=reshape(my_code2,N,length(my_code)/(m*N));
  149.  my_code4=my_code4';%加噪后的码字,一行即为一个码组  
  150.  %===============================================中间的式程中不用
  151.  
  152.  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  153.  %译码,BM迭代,时域法
  154.  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  155. %计算伴随式 
  156. for i=1:length(my_code)/(N*m) 
  157.     error_flag(i)=0;
  158. S=zeros(8,7);%纠4个错,需要计算8个伴随式,每个伴随式用7位表示
  159. for j=1:N
  160.     r(j)=my_code2(i*N-N+j);
  161.     for k=1:m
  162.         rr(1+m-k)=mod(r(j),2);
  163.         r(j)=(r(j)-rr(1+m-k))/2;
  164.     end,k;
  165.     for k=1:8
  166.         S0=S(k,1)*64+S(k,2)*32+S(k,3)*16+S(k,4)*8+S(k,5)*4+S(k,6)*2+S(k,7);
  167.         if(S0==0)  S(k,:)=rr;
  168.         else  S0_index=int_exp(S0);
  169.               S0_index=mod((S0_index+k),127);%当m0,h改变时,校验矩阵也是要变化的,见P243
  170.               SS=exp_int_2(S0_index+1,:);
  171.               S(k,:)=xor(SS,rr);
  172.         end;
  173.     end,k;
  174. end,j;
  175. %上面计算好了伴随式,下面求错误多项式
  176.     
  177. %==========================================================
  178. %==========================================================Euclids
  179.   Ri(1,:)=[0,0,0,0,0,0,0];Qi(1,:)=S(1,:);
  180.   Li(1,:)=[0,0,0,0,0,0,0];Ui(1,:)=[0,0,0,0,0,0,1];
  181.   for k=2:8%2*t+1=9
  182.       Ri(k,:)=[0,0,0,0,0,0,0];Qi(k,:)=S(k,:);
  183.       Li(k,:)=[0,0,0,0,0,0,0];Ui(k,:)=[0,0,0,0,0,0,0];
  184.   end,k;
  185.   Li(9,:)=[0,0,0,0,0,0,0];Ui(9,:)=[0,0,0,0,0,0,0];
  186.   Qi(9,:)=[0,0,0,0,0,0,0];Ri(9,:)=[0,0,0,0,0,0,1];%R0,Q0,L0,U0初始化
  187.   
  188.   deg_Ri=8;deg_Qi=0;  %Ri,Qi的deg初始化
  189. %=====================================
  190.   for k=9:-1:1
  191.       Qi0=Qi(k,1)*64+Qi(k,2)*32+Qi(k,3)*16+Qi(k,4)*8+Qi(k,5)*4+Qi(k,6)*2+Qi(k,7);
  192.       if(Qi0~=0) deg_Qi=k-1;break;
  193.       end
  194.   end,k;%搜索Qi的deg
  195.   
  196.  %====================================
  197.   
  198.  while((deg_Ri>=4)&&(deg_Qi>=4)) 
  199.      li=deg_Ri-deg_Qi;  %计算li
  200.      if(li>=0) sw=1;
  201.      else sw=0;
  202.      end;  %计算delta_i-1;
  203.      li=abs(li);
  204.      %========================
  205.      if(sw==1) 
  206.          for k=1:9
  207.              Ri2(k,:)=Ri(k,:);Qi2(k,:)=Qi(k,:);Li2(k,:)=Li(k,:);Ui2(k,:)=Ui(k,:);
  208.          end,k;
  209.      else
  210.          for k=1:9
  211.              Ri2(k,:)=Qi(k,:);Qi2(k,:)=Ri(k,:);Li2(k,:)=Ui(k,:);Ui2(k,:)=Li(k,:);
  212.          end,k;
  213.      end    %第一次进行交换
  214.      %=======================
  215.      deg_Ri=0;deg_Qi=0;a5=0;b5=0; %暂时没找到deg_Ri,deg_Qi之间的关系,故每次都重新搜索
  216.      for k=9:-1:1
  217.          Ri0=Ri2(k,1)*64+Ri2(k,2)*32+Ri2(k,3)*16+Ri2(k,4)*8+Ri2(k,5)*4+Ri2(k,6)*2+Ri2(k,7);
  218.          if(Ri0~=0) a5=Ri0;break;end;
  219.      end,k;
  220.      for k=9:-1:1
  221.           Qi0=Qi2(k,1)*64+Qi2(k,2)*32+Qi2(k,3)*16+Qi2(k,4)*8+Qi2(k,5)*4+Qi2(k,6)*2+Qi2(k,7);
  222.           if(Qi0~=0) b5=Qi0;break;end;
  223.      end,k;%搜索Ri,Qi的deg,以及首系数。如果为全,则赋值为0
  224.      %========================
  225.      
  226.      if(b5==0) 
  227.          for k=1:9 
  228.              Ri3(k,:)=[0,0,0,0,0,0,0];Li3(k,:)=[0,0,0,0,0,0,0];
  229.          end,k;
  230.      else b5_index=int_exp(b5);
  231.         for k=1:9 
  232.            r=Ri2(k,1)*64+Ri2(k,2)*32+Ri2(k,3)*16+Ri2(k,4)*8+Ri2(k,5)*4+Ri2(k,6)*2+Ri2(k,7);
  233.            ll=Li2(k,1)*64+Li2(k,2)*32+Li2(k,3)*16+Li2(k,4)*8+Li2(k,5)*4+Li2(k,6)*2+Li2(k,7);
  234.            if(r==0) Ri3(k,:)=[0,0,0,0,0,0,0];
  235.            else  r_index=int_exp(r);r_index=mod(r_index+b5_index,127);Ri3(k,:)=exp_int_2(r_index+1,:);
  236.            end;
  237.            if(ll==0) Li3(k,:)=[0,0,0,0,0,0,0];
  238.            else ll_index=int_exp(ll);ll_index=mod(ll_index+b5_index,127);Li3(k,:)=exp_int_2(ll_index+1,:);
  239.            end;
  240.        end,k;%计算公式的前半部分,R,L
  241.      end;
  242.     %=========================== 
  243.      if(a5==0) 
  244.          for k=1:9
  245.              Ri4(k,:)=[0,0,0,0,0,0,0];Li4(k,:)=[0,0,0,0,0,0,0];
  246.          end,k;
  247.      else a5_index=int_exp(a5);
  248.         for k=1:9-li
  249.            r=Qi2(k,1)*64+Qi2(k,2)*32+Qi2(k,3)*16+Qi2(k,4)*8+Qi2(k,5)*4+Qi2(k,6)*2+Qi2(k,7);
  250.            ll=Ui2(k,1)*64+Ui2(k,2)*32+Ui2(k,3)*16+Ui2(k,4)*8+Ui2(k,5)*4+Ui2(k,6)*2+Ui2(k,7);
  251.            if(r==0) Ri4(k+li,:)=[0,0,0,0,0,0,0];
  252.            else  r_index=int_exp(r); r_index=mod(r_index+a5_index,127);Ri4(k+li,:)=exp_int_2(r_index+1,:);
  253.            end 
  254.            if(ll==0) Li4(k+li,:)=[0,0,0,0,0,0,0];
  255.            else ll_index=int_exp(ll);ll_index=mod(ll_index+a5_index,127);Li4(k+li,:)=exp_int_2(ll_index+1,:);
  256.            end
  257.        end,k;%计算公式后半部分
  258.      end;
  259.    %=============================  
  260.      for k=1:li
  261.          Ri4(k,:)=[0,0,0,0,0,0,0];Li4(k,:)=[0,0,0,0,0,0,0];
  262.      end,k;
  263.      
  264.      for k=1:9
  265.          Ri(k,:)=xor(Ri3(k,:),Ri4(k,:));Li(k,:)=xor(Li3(k,:),Li4(k,:));
  266.          Qi(k,:)=Qi2(k,:);Ui(k,:)=Ui2(k,:);
  267.      end,k;%相加,并返回给Ri,Qi,Li,Ui
  268.      
  269.      for k=9:-1:1
  270.          Ri0=Ri(k,1)*64+Ri(k,2)*32+Ri(k,3)*16+Ri(k,4)*8+Ri(k,5)*4+Ri(k,6)*2+Ri(k,7);
  271.          if(Ri0~=0) deg_Ri=k-1;break;end;
  272.      end,k;
  273.      for k=9:-1:1
  274.           Qi0=Qi(k,1)*64+Qi(k,2)*32+Qi(k,3)*16+Qi(k,4)*8+Qi(k,5)*4+Qi(k,6)*2+Qi(k,7);
  275.           if(Qi0~=0) deg_Qi=k-1;break;end;
  276.      end,k;
  277.      if((deg_Ri<4)||(deg_Qi<4)) break;end;
  278.  end;
  279.  
  280.  for k=1:5
  281.      a2(k,:)=Li(k,:);w2(k,:)=Ri(k,:);
  282.  end,k;%将值赋给a2,w2(后续的chien搜索和Forney算法中用的是a2,w2)
  283. %==================================迭代计算结束,但结果不对
  284. for k=1:5
  285.     a20(k)=a2(k,1)*64+a2(k,2)*32+a2(k,3)*16+a2(k,4)*8+a2(k,5)*4+a2(k,6)*2+a2(k,7);
  286. end,k;
  287. if(a20(5)~=0) cofficent=4;
  288. else if(a20(4)~=0) cofficent=3;
  289. else if(a20(3)~=0) cofficent=2;
  290. else if(a20(2)~=0) cofficent=1;
  291. else cofficent=0;
  292. end;end;end;end;
  293. %
  294. %迭代算法到此结束,错误位置多项式存在a2中,错误值多项式存在w2中
  295. %==================================================
  296. %==================================================
  297. %钱搜索以及Forney算法
  298.    n=0;%n用来计错误个数
  299.    
  300. chien(1,:)=a2(1,:);chien(2,:)=a2(3,:);chien(3,:)=a2(5,:);chien(4,:)=a2(2,:);chien(5,:)=a2(4,:);
  301. chien_g_index=[0,2,4,0,2];
  302. %Multi=[0,0,0,0,0,0,1];
  303.     Multi=[0,1,0,1,1,0,1];%因为缩短了55个,可能应该从a^55开始 
  304. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  305. %%调整初值
  306. for l=1:5
  307.     chien0=chien(l,1)*64+chien(l,2)*32+chien(l,3)*16+chien(l,4)*8+chien(l,5)*4+chien(l,6)*2+chien(l,7);
  308.     if(chien0==0) chien(l,:)=[0,0,0,0,0,0,0];
  309.     else chien_index(l)=int_exp(chien0);
  310.          chien_index(l)=mod((chien_index(l)+chien_g_index(l)*55),127);
  311.          chien(l,:)=exp_int_2(chien_index(l)+1,:);
  312.      end;
  313. end,l;%先将初始值调整直55位置  
  314. for l=1:5
  315.       w20=w2(l,1)*64+w2(l,2)*32+w2(l,3)*16+w2(l,4)*8+w2(l,5)*4+w2(l,6)*2+w2(l,7);
  316.       if(w20==0)  w2=w2;
  317.       else  w2_index=int_exp(w20);
  318.            w2_index=mod(w2_index+(l-1)*55,127);
  319.            w2(l,:)=exp_int_2(w2_index+1,:);
  320.       end;
  321.    end,l;%计算w(a^i)的值,错误值位置也调整至55
  322. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  323. %对有用位置钱搜索
  324. for k=127+1-N:127
  325.     for l=1:5
  326.         chien0=chien(l,1)*64+chien(l,2)*32+chien(l,3)*16+chien(l,4)*8+chien(l,5)*4+chien(l,6)*2+chien(l,7);
  327.         if(chien0==0) chien(l,:)=[0,0,0,0,0,0,0];
  328.         else chien_index(l)=int_exp(chien0);
  329.             chien_index(l)=mod((chien_index(l)+chien_g_index(l)),127);
  330.             chien(l,:)=exp_int_2(chien_index(l)+1,:);
  331.         end;
  332.     end,l;%先相乘一次      %钱搜索位置前提55个
  333.     
  334.     w_sum=[0,0,0,0,0,0,0];
  335.     fenmu=[0,0,0,0,0,0,0];
  336.     for l=4:5
  337.         fenmu=xor(fenmu,chien(l,:));
  338.     end,l;%上面的再乘一个a,并算出a'的值
  339.   
  340.    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  341.    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  342.    %下面这段是计算钱搜索结果,是否为0,即是否检测到错误(0代表检测到错)
  343.    Multi0=Multi(1)*64+Multi(2)*32+Multi(3)*16+Multi(4)*8+Multi(5)*4+Multi(6)*2+Multi(7);%乘到奇数次幂的系数上
  344.    if(Multi0==0) Multi=[0,0,0,0,0,0,0];
  345.    else  Multi_index=int_exp(Multi0);
  346.          Multi_index=mod(Multi_index+1,127);%h=1
  347.          Multi=exp_int_2(Multi_index+1,:);
  348.    end;
  349.     fenmu0=fenmu(1)*64+fenmu(2)*32+fenmu(3)*16+fenmu(4)*8+fenmu(5)*4+fenmu(6)*2+fenmu(7);
  350.     Multi0=Multi(1)*64+Multi(2)*32+Multi(3)*16+Multi(4)*8+Multi(5)*4+Multi(6)*2+Multi(7);
  351.     if(fenmu0==0|Multi0==0) chien1=[0,0,0,0,0,0,0];%fenmu即为奇数次幂系数的导数
  352.     else fenmu_index=int_exp(fenmu0);
  353.         Multi_index=int_exp(Multi0);
  354.         chien1_index=mod(fenmu_index+Multi_index,127);
  355.         chien1=exp_int_2(chien1_index+1,:);
  356.      end;
  357.     
  358.     fai_sum=[0,0,0,0,0,0,0];
  359.     for l=1:3
  360.         fai_sum=xor(fai_sum,chien(l,:));
  361.     end,l;
  362.     fai_sum=xor(fai_sum,chien1);%计算出钱搜索的结果
  363.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  364.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  365.   
  366.   %================================================
  367.   %计算w(a^i)的值
  368.     for l=1:5
  369.       w20=w2(l,1)*64+w2(l,2)*32+w2(l,3)*16+w2(l,4)*8+w2(l,5)*4+w2(l,6)*2+w2(l,7);
  370.       if(w20==0)  w_sum=w_sum;
  371.       else  w2_index=int_exp(w20);
  372.            w2_index=mod((w2_index+l-1),127);
  373.            w2(l,:)=exp_int_2(w2_index+1,:);
  374.            w_sum=xor(w_sum,w2(l,:));
  375.        end;
  376.    end,l;%计算w(a^i)的值
  377.    %==================================================
  378.    
  379.    fai_sum0=fai_sum(1)*64+fai_sum(2)*32+fai_sum(3)*16+fai_sum(4)*8+fai_sum(5)*4+fai_sum(6)*2+fai_sum(7);
  380.    w_sum0=w_sum(1)*64+w_sum(2)*32+w_sum(3)*16+w_sum(4)*8+w_sum(5)*4+w_sum(6)*2+w_sum(7);  
  381.    if(fai_sum0~=0) M(m*N*i-m*N+m*(k-55)-m+1:m*N*i-m*N+m*(k-55))=my_code(m*N*i-m*N+m*(k-55)-m+1:m*N*i-m*N+m*(k-55));%没检测到错误,,码字直接输出
  382.    else 
  383.    n=n+1;%错误个数加1
  384.         if(fenmu0==0|w_sum0==0) M(m*N*i-m*N+m*(k-55)-m+1:m*N*i-m*N+m*(k-55))=my_code(m*N*i-m*N+m*(k-55)-m+1:m*N*i-m*N+m*(k-55)); %分子或者分母为0,直接输出
  385.         else
  386.         fenmu_index=int_exp(fenmu0);
  387.         w_sum0_index=int_exp(w_sum0);
  388.         error_index=mod((w_sum0_index-fenmu_index),127);
  389.         error=exp_int_2(error_index+1,:);%按公式计算错误值
  390.         M(m*N*i-m*N+m*(k-55)-m+1:m*N*i-m*N+m*(k-55))=xor(my_code(m*N*i-m*N+m*(k-55)-m+1:m*N*i-m*N+m*(k-55)),error);%错误值跟码字异或就是纠正后的输出
  391.     end;
  392. end;
  393. end,k;
  394. if(n<cofficent)
  395.   error_flag(i)=1;  %若检测到的错误数小于错误位置多项式的系数就认为错误个数超出了可纠范围,不能正确译码,给出标志
  396. end;
  397. end,i;%大循环,处理i个码组
  398.    
  399. for i=1:(length(M)/m)
  400.     M2(i)=M(m*i-m+1)*64+M(m*i-m+2)*32+M(m*i-m+3)*16+M(m*i-m+4)*8+M(m*i-m+5)*4+M(m*i-m+6)*2+M(m*i-m+7);
  401. end,i;  %转化成十进制,按码字显示
  402.  M2=reshape(M2,N,length(M2)/N);
  403.  M2=M2';
  404. c2=my_code3-M2;
  405. y=c2~=0;
  406. c2(:,N+1)=sum_noise;c2(:,N+2)=error_flag';c2(:,N+3)=sum(y,2)
  407. %c2前N列是译码后数据跟编码后数据(加错前)的比较,0表示正确的译码.N+1列是错码个数,
  408. %调试结果正确,编码译码都在里面
  409.  fid=fopen('c:decode_7264_out.txt','w');
  410.  for i=1:length(M2)
  411.     fprintf(fid,'%3d ',M2(i));
  412.     if(mod(mod(i,72),18)==0) fprintf(fid,'n');
  413.     end
  414.     if(mod(i,N)==0) fprintf(fid,'n');
  415.     end
  416.  end,i;        
  417.  fclose(fid);
  418. aaa=decode_7264(my_code);
  419. for i=1:(length(aaa)/m)
  420.     aaa2(i)=aaa(m*i-m+1)*64+aaa(m*i-m+2)*32+aaa(m*i-m+3)*16+aaa(m*i-m+4)*8+aaa(m*i-m+5)*4+aaa(m*i-m+6)*2+aaa(m*i-m+7);
  421. end,i;  %转化成十进制的码字
  422. aaa2=reshape(aaa2,N,length(aaa2)/N);
  423.  aaa2=aaa2';
  424.  ccc=M2-aaa2
  425.  
  426.  
  427.  %结果表明译码结果正确,采用Euclids算法和BM迭代法时后续的Forney算法是不一样的,前者是w(x^(-i))/a'(x^(-i)),后者是x^i*w(x^(-i))/a'(x^(-i))