rgb=imread('g:\pp.bmp','bmp');
pic=imresize(rgb,[64 64]);
st=zeros(1,4096);
for i=0:63
    for j=0:63
         st(1,i*64+j+1)=pic(i+1,j+1); %str(1,(64*(i-1)+j))=pic(i,j);
   end
end                                                   
s=zeros(1024,4);
for i=0:1023    
    for j=0:3        
        s(i+1,j+1)=st(i*4+j+1);    
    end
end
g=[1 0 1 1];
codecrc=encode(s,7,4,'cyclic',g);%codecrc 7168=4096*7/4
codebch=bchenco(codecrc,15,7);%7168*15/7=1024*15=15360

code=zeros(1,15360);
 for i=0:1023
    for j=0:14
         code(1,i*15+j+1)=codebch(i+1,j+1); %str(1,(64*(i-1)+j))=pic(i,j);
   end
end          

aa=wavread('22');
[sqnr256,aquan256,code256]=u_pcm(aa,256);%注意量化位數的選取！！！這里選取8，2^8=256
A1=aquan256;
AA=A1';
A(1:130000)=AA(20001:150000);
for i=1:122880 % 15360*8=122880,錄音的時候采樣頻率是22.05KHZ
     Ae(i)=A(i);
end

for k=1:15360
    Aef(k,(1:8))=Ae(((k-1)*8+1):((k-1)*8+8));
    
    y(k,(1:8))=dct(Aef(k,(1:8)));
    
    yy(k,2)=y(k,2)*(1+0.8*code(k));
    
    yy(k,1)=y(k,1);
    yy(k,(3:8))=y(k,(3:8));
    
    Aee(k,(1:8))=idct( yy(k,(1:8)));
    Aey(1,(((k-1)*8+1):((k-1)*8+8)))=Aee(k,(1:8));
end

for  i=1:122880;%64*64*8=32768錄音的時候采樣頻率是22.05KHZ
     Aw(i)=Aey(i);
end
for i=122881:130000;%64*64*8=32768錄音的時候采樣頻率是22.05KHZ
     Aw(i)=A(i);
end

p=0.0001;
N=130000;%注意N的選取
u=randn(1,N);
a=sqrt(p)
u=u*a; 
As=Aw+u;%加入白噪聲以后，待檢測信號就是音頻信號與噪聲相加

for q=1:122880% 48*8=384
     Ase(q)=As(q);
end

for j=1:15360
     Aeff(j,(1:8))=Ase(((j-1)*8+1):((j-1)*8+8));
     
     yyy(j,(1:8))=dct(Aeff(j,(1:8)));
     
     strr(j)=[yyy(j,2)-y(j,2)]/[0.8*y(j,2)];
     strrr=abs(strr);
end
minstr=min(strrr);
maxstr=max(strrr)
b=strrr;
for  p=1:15360 %此處沒有錯 
     if  b(p)>0.60  % threshoud 是自己設定的門限值，可以取得0.5到0.7之間，但是0.6效果最好。
           b(p)=1; %似乎是1表示白，0表示
     else  
           b(p)=0;
     end 
end 
bb=zeros(1024,15);
for i=0:1023    
    for j=0:14        
        bb(i+1,j+1)=b(i*15+j+1);    
    end
end
codebchr=bchdeco(bb,7,1);
codecrcr=decode(codebchr,7,4,'cyclic',g); %輸出糾正的碼字

sr=zeros(1,4096);
for i=0:1023
    for j=0:3
         sr(1,i*4+j+1)=codecrcr(i+1,j+1); %str(1,(64*(i-1)+j))=pic(i,j);
   end
end                                                   
picc=zeros(64,64);
for i=0:63    
    for j=0:63        
        picc(i+1,j+1)=sr(i*64+j+1);    
    end
end
figure(1)
subplot(211),imshow(pic);title('pic');
subplot(212),imshow(picc);title('picc');
figure(2)
subplot(2,1,1);plot(Ae);title('Ae');
subplot(2,1,2);plot(As);title('As');