Polar编码matlab程序

function TEB=polart2 N_capteurs=4; N_relais=N_capteurs; n=N_capteurs*2; SNR=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20]; Nb_pk=1000; G2=[1 0; 1 1]; if n==2 GG=G2; elseif n==4 GG=kron(G2,G2); elseif n==8 GG=kron(kron(G2,G2),G2); elseif n==16 GG=kron(kron(kron(G2,G2),G2),G2); elseif n==32 GG=kron(kron(kron(kron(G2,G2),G2),G2),G2); elseif n==64 GG=kron(kron(kron(kron(kron(G2,G2),G2),G2),G2),G2); elseif n==128 GG=kron(kron(kron(kron(kron(kron(G2,G2),G2),G2),G2),G2),G2); elseif n==256 GG=kron(kron(kron(kron(kron(kron(kron(G2,G2),G2),G2),G2),G2),G2),G2); elseif n==512 GG=kron(kron(kron(kron(kron(kron(kron(kron(G2,G2),G2),G2),G2),G2),G2),G2),G2); elseif n==1024 GG=kron(kron(kron(kron(kron(kron(kron(kron(kron(G2,G2),G2),G2),G2),G2),G2),G2),G2),G2); elseif n==2048 GG=kron(kron(kron(kron(kron(kron(kron(kron(kron(kron(G2,G2),G2),G2),G2),G2),G2),G2),G2),G2),G2); end bhat=zeros(n,n); bhat(1,1)=0.5; for a=2:2:n for b=1:n if (1<=b) && (b<=a/2) bhat(a,b)=(2*bhat(a/2,b)) - (bhat(a/2,b))^2; else if ((a/2)+1<=b) && (b<=a) bhat(a,b)=(bhat(a/2,b-(a/2)))^2; else bhat(a,b)=0; end end end end bhat=bhat(n,:); bhat; [val, index]=sort(bhat); lignes=index(1:n/2); lignes=sort(lignes); G=GG(lignes,:); coding_rate=1/2; for m=1:size(SNR,2) for i=1:Nb_pk sigma(m) = sqrt(10^(-SNR(:,m)/10)*(1/coding_rate)); % AWGN noise deviation u(i,:) = (sign(randn(1,size(G,1)))+1)/2; % random bits y(i,:) = mod(u(i,:)*G,2); % coding z(i,:) = 2*y(i,:)-1; channel=(1/sqrt(2))*(randn(1,n)+j*randn(1,n)); % channel=randn(1,n); % channel=ones(1,n); % BPSK modulation zd(i,:)=z(i,:).*channel; % zd(i,:)=zd(i,:)+ (sigma(m)*randn(1,n)); zd(i,:)=zd(i,:)+sigma(m)*(1/sqrt(2))*(randn(1,size(G,2))+j*randn(1,size(G,2))); % AWGN transmission ze(i,:)=conj(channel).*zd(i,:); L1=4*real(ze(i,1))/sigma(m)^2; L2=4*real(ze(i,2))/sigma(m)^2; L3=4*real(ze(i,3))/sigma(m)^2; L4=4*real(ze(i,4))/sigma(m)^2; L5=4*real(ze(i,5))/sigma(m)^2; L6=4*real(ze(i,6))/sigma(m)^2; L7=4*real(ze(i,7))/sigma(m)^2; L8=4*real(ze(i,8))/sigma(m)^2; if (L1>709) L1=709; end if (L2>709) L2=709; end if (L3>709) L3=709; end if (L4>709) L4=709; end if (L5>709) L5=709; end if (L6>709) L6=709; end if (L7>709) L7=709; end if (L8>709) L8=709; end LLR4 = log((exp(L3)+exp(L7))/(1+ exp((L3)*(L7)))) + log((exp(L4)+exp(L8))/(1+ exp((L4)*(L8)))) if (LLR4>0) u4=1; else u4=0; end LLR61=log((exp(L3)+exp(L1))/(1+exp(L3+L1))); LLR62=log((exp(L7)+exp(L5))/(1+exp(L7+L5))); LLR63=log((exp(L4)+exp(L2))/(1+exp(L4+L2))); LLR64=log((exp(L8)+exp(L6))/(1+exp(L8+L6))); LLR6 = LLR61 + LLR62 + LLR63 + LLR64 if (LLR6>0) u6=1; else u6=0; end LT3= L3+ ((1-2*u6)*L1) T3=exp(LT3) LT7= L7+ ((1-2*u6)*L5) T7=exp(LT7) LT4= L4+ ((1-2*u6)*L2) T4=exp(LT4) LT8= L8+ ((1-2*u6)*L6) T8=exp(LT8) LS7= L7+ ((1-2*u4)*LT3) S7=exp(LS7) LS8= L8+ ((1-2*u4)*LT4) S8=exp(LS8) LLR7=log((exp(LS7)+exp(LS8))/(1+ exp((LS7)*(LS8)))); if (LLR7>0) u7=1; else u7=0; end LLR8=LS8+((1-2*u7)*LS7) if (LLR8>0) u8=1; else u8=0; end u_hat=[u4 u6 u7 u8]; o = xor(u_hat,u(i,:)); numerror(m,i) = sum(o); BER(m,i) = numerror(m,i)/size(G,1); end; TEB(m)=(1/Nb_pk)*sum(BER(m,:),2) end; semilogy(SNR,TEB(:)); xlabel('SNR (en dB)'); ylabel('TEB'); title('TEB des codes polaires'); grid on; return;