1 clc;clear
 2 %%这个文件，是最终版本，
 3 %% 此版块内，n=4，N由1到10,间隔为1，
 4 %%(i+0.5)*deierka)这一行，中的i有问题，已经用平移量解决.
 5 %%
 6 kesai=0.1;                   %过量噪声
 7 seigema=sqrt(1+kesai);
 8 b_lim=7*seigema;          %  a_lim用于后面计算S的过程里，就是区间长度的1半，根据文章来的。
 9 length=1000;                 %产生随机数组的次数，即100
10 for n=4:4:16                %外层循环是n的循环，n是零差检测精度
11     i_min=-2^(n-1);i_max=2^(n-1)-1;
12     PingYi=1-i_min;         %为了把角标变为正数，所用的平移量。比如要把-8到7平移成1到16，那么平移量就是1-（-8）=9
13 for N=0.1:0.1:10                  % 内层循环是N的循环，N采样范围
14     deierka=N/2^(n-1);  m=0;         %是间隔，m用于计算后面的a_ba
15   %%这一块在计算P和H
16     for i=1:i_max-i_min+1
17         a(i)=-N+(i-1)*deierka;   %deierka就是步进间隔
18     end
19     for i=1:i_max-i_min+1
20         if i==1
21            P(i)=(1/2)*erfc((N-0.5*deierka)/(sqrt(2)*seigema));
22         elseif i>1&&i<i_max-i_min+1
23            P(i)=(1/2)*erf(((i-PingYi+0.5)*deierka)/(sqrt(2)*seigema))-(1/2)*erf(((i-PingYi-0.5)*deierka)/(sqrt(2)*seigema));
24         elseif i==i_max-i_min+1
25            P(i)=(1/2)*erfc((N-1.5*deierka)/(sqrt(2)*seigema));
26         end
27     end
28           index=int8(10*N);
29         H(index)=0;  R(index)=0  %清零操作
30         for i=1:i_max-i_min+1
31             diedai=P(i)*log2(P(i))
32             if P(i)==0
33                 diedai=0;
34             end
35             H(index)=H(index)-diedai;
36         end
37   %%    这一块在计算S(a：E)
38       b_min_buf=0;  b_max_buf=0;  b_ba_buf=0;
39       buffer1=0;  buffer2=0  %  两个buffer用于计算式子C.3中的两个求和
40      for k=1:length
41          b=-b_lim+2*b_lim*rand(1,1000);
42          b_min_buf=b_min_buf+min(b);  b_max_buf=b_max_buf+max(b);
43          b_ba_buf=b_ba_buf+mean(b);
44      end
45          b_min=b_min_buf/length;  b_max=b_max_buf/length;
46          b_ba=mean(a);
47          for i=2:i_max-i_min   %不从1开始。另外，倒数第二个是结尾
48            if a(i)<=b_ba
49                 buffer1=buffer1+P(i)*(a(i)-b_ba-0.5*deierka)^2;
50            else
51                 buffer2=buffer2+P(i)*(a(i)-b_ba+0.5*deierka)^2;
52            end
53          end
54          Vx_ba(index)= P(1)*(b_min-b_ba)^2+P(i_max-i_min+1)*(b_max-b_ba)^2+buffer1+buffer2;
55          S(index)=((Vx_ba(index)+1)/2)*log2((Vx_ba(index)+1)/2)-((Vx_ba(index)-1)/2)*log2((Vx_ba(index)-1)/2);
56     %  这一块在计算R，
57           R(index)=H(index)-S(index);
58
59 end
60          t=0.1:0.1:10;
61          if n==4
62             plot(t,R,'r*');
63          end
64          if n==8
65             plot(t,R,'*');
66          end
67          if n==12
68             plot(t,R,'d');
69          end
70          if n==16
71             plot(t,R,'p');
72          end
73          legend('n=4','n=8','n=12','n=16');
74          xlim([1,10]);ylim([0,24]);
75          xlabel('Sample range');
76          ylabel('R_dis（ai:E');
77          hold on
78 end
79 grid