虽然跑出来了,但没有完全跑出来.jpg
原来是想看看自带的rs编解码函数长什么样,但解码部分不会看,好吧
利用自带的GF域来编了一遍,自带的gf类相关计算函数都在@gf文件夹里,最近去看其他东西了暂时还没来得及研究。
函数列表
函数名 | 功能 |
---|---|
Berlekamp2(syn,M,N,K) | Berlekamp迭代算法求错误位置多项式 |
bin2GFdec(input,k,m) | 二进制转GF域 |
errorcorrect (msg,errorst,pos) | 错误位置 |
errornum2( errorpoly,syn,N,GF_index) | 错误数值 |
GFdec2bin( rs_code,m ) | GF域转二进制 |
GFindex(M) | GF域幂次索引 |
rsdec2(dec_msg,N,K) | rs解码 |
rsenc2(msg, N, K) | rs编码 |
rsgenpoly2(N, K) | 生成多项式(没用到,用的自带的) |
syndrome2(msg,N,K) | 校正子计算 |
主程序
会跑很久,信息长度缩短能快些,但那样差距不太能体现出来。
clc;
clear;
len=100000;
m=8; %扩域GF(2^m)
n=2^m-1; %2^m进制编码长度;
k=223; %信息位长度;
t=(n-k)/2; %纠错数
msg=randi([0 1], len, 1);
mod_bpsk = comm.BPSKModulator;
demod = comm.BPSKDemodulator;
%rs编码
demsg=bin2GFdec(msg,k,m);
rs_code=rsenc2(demsg,n,k);
data=GFdec2bin( rs_code,m );
%BPSK调制-高斯信道
modSignal = step(mod_bpsk,data); %RS编码
modSignal_2=step(mod_bpsk,msg); %无编码
N=1;
for i=1:10
chan = comm.AWGNChannel('NoiseMethod','Signal to noise ratio (SNR)','SNR',i);
noisySignal = step(chan,modSignal);
receivedData = step(demod,noisySignal);
noisySignal_2 = step(chan,modSignal_2);
receivedData_2 = step(demod,noisySignal_2);
%无编码
errorRate = comm.ErrorRate;
errorStats = step(errorRate,msg,receivedData_2);
BER_BPSK(N)=errorStats(1);
fprintf('no RS Error rate = %f\nNumber of errors = %d\n', errorStats(1), errorStats(2))
%rs译码1
dec_msg=bin2GFdec(receivedData,n,m);
decrs_code = rsdec2(dec_msg,n,k);
decrsbin_code=GFdec2bin( decrs_code,m );
recmsg=decrsbin_code(1:len);
errorRate = comm.ErrorRate;
errorStats = step(errorRate,msg,recmsg);
BER_BPSK_RS(N)=errorStats(1);
fprintf('Error rate = %f\nNumber of errors = %d\n', errorStats(1), errorStats(2))
fprintf('\n')
Y(N)=i;
N=N+1;
end
figure
semilogy(Y,BER_BPSK,'-o',Y,BER_BPSK_RS,'-o');
legend('无编码','(255,223)RS码');
grid on
xlabel('SNR dB');
ylabel('BER');
函数
RS编码
function code = rsenc2(msg, N, K)
% Find fundamental parameters m and t
M = msg.m;
t = (N-K)/2;
T2 = 2*t; % number of parity symbols
[m_msg, n_msg] = size(msg);
genpoly = rsgenpoly(N,K,msg.prim_poly); % 生成多项式
parity = gf(zeros(m_msg,T2),M,msg.prim_poly);
%生成多项式最高项为1
genpoly = genpoly(2:T2+1);
% Encoding
% Each row gives the parity symbols for each message word
for j=1:size(msg,2),
parity = [parity(:,2:T2) zeros(m_msg,1)] + (msg(:,j)+parity(:,1))*genpoly; %gf多项式除法
end
code = [msg parity];
RS译码
function decrs_code = rsdec2(dec_msg,N,K)
M=dec_msg.m;
GF_index = GFindex(M);
for i=1:size(dec_msg,1)
msg=dec_msg(i,:);
syn = syndrome2(msg,N,K);
errorpoly = Berlekamp2(syn,M,N,K);
[errorst pos] = errornum2( errorpoly,syn,N,GF_index);
errccor = errorcorrect ( msg,errorst,pos);
rmsg=errccor(1:K);
cor(i,:)=rmsg.x;
end
decrs_code=gf(cor,M);
end
生成多项式
function [genpoly,t] = rsgenpoly2(N, K)
t = floor((N-K)/2); %向下取整
t2 = N-K;
m = log2(N+1);
b = 1;
alpha = gf(2,m);
genpoly = 1;
for k=1:t2
evalc('genpoly = conv(genpoly,[1 alpha.^(k)]);');
end
校正子计算
跑了下发现这里用的时间最长,引用polyval的次数太多了…
function syn = syndrome2(msg,N,K)
M=msg.m;
t2=N-K;
alpha=gf(2,M);
for i=1:t2
a=alpha.^i;
sym=polyval(msg,a);
p(i)=sym.x;
end
syn=gf(p,M);
end
伯利坎普迭代
function errorpoly = Berlekamp2(syn,M,N,K)
t=(N-K)/2;
u=1; %迭代步数初始化
errorpoly=gf(1,M); %错误位置多项式初始化
du=syn(1); %迭代步差
lup=0;
lu=0; %errorpoly最高幂次
p=-1;
dp=gf(1,M);; %ρ初始化
errp=gf(1,M);
for u=1:2*t
errp2=errorpoly;
xerrp=errp;
if du==gf(0,M) %du=0
errorpoly=errorpoly;
else
if u-1-p>0
xerrp(size(errp,2)+u-1-p)=0;
end
xerrp=du*dp.^(-1)*xerrp;
if length(xerrp)>length(errorpoly) %多项式加法长度保持一致
a=length(xerrp)-length(errorpoly);
errorpoly(a+1:length(errorpoly)+a)=errorpoly;
errorpoly(1:a)=gf(0,M);
elseif length(xerrp)<length(errorpoly)
a=length(errorpoly)-length(xerrp);
xerrp(a+1:length(xerrp)+a)=xerrp;
xerrp(1:a)=gf(0,M);
end
errorpoly=errorpoly+xerrp;
end
if (p-lup<u-1-lu)&&(du~=gf(0,M)) %选ρ
dp=du;
p=u-1;
errp=errp2;
lup=lu;
end
lu=size(errorpoly,2)-1;
if u<2*t
du=syn(u+1);
for b=1:lu
if size(errorpoly,2)==1
du=du;
else
du=du+errorpoly(end-b)*syn(u+1-b);
end
end
end
end
end
错误数值
function [errorst pos] = errornum2( errorpoly,syn,N,GF_index)
%求错误数值
%errorpoly-错误位置多项式
%syn-校正子
%GF-GF域
%errorst-错误多项式
%pos-错误位置
M=errorpoly.m;
root=roots(errorpoly);
syn1=syn(end:-1:1);
if length(root)==0 %如果没有根,则错误多项式为0
errorst=0;
pos=0;
else
Z=conv(errorpoly,syn1);
Z0=Z(end-length(errorpoly)+2:end);
for i=1:length(root)
ZB=polyval(Z0,root(i));
sigma=gf(1,M);
for k=1:length(root)
if root(k)~=root(i)
sigma=sigma*(1+(root(k)).^(-1)*root(i));
end
end
sigma=sigma*(root(i)).^(-1);
error=ZB/sigma;
errornum(i)=error.x;
pos1=root(i);
if pos1==gf(1,M)
pos(i)=N;
else
pos(i)=GF_index(pos1.x+1);
end
end
errorst=gf(errornum,M);
end
end
错误位置矫正
function errccor = errorcorrect (msg,errorst,pos)
errccor=msg;
if pos==0
else
for i=1:length(errorst)
errccor(pos(i))=msg(pos(i))+errorst(i);
end
end
end
二进制转GF域
function msg= bin2GFdec(input,k,m)
%二进制分帧并转GF域
L=length(input);
z=fix(L/(k*m));
r=mod(L,k*m);
r1=k*m-mod(L,k*m);
if r~=0 %若原始信息数除不尽信息位,则用0填满
input(L+1:L+r1)=0;
z=z+1;
end
inputGF=zeros(length(input)/m,1);
for i=1:length(input)/m
for x=1:m
inputGF(i)=2^(m-x)*input((i-1)*m+x)+inputGF(i);
end
end
msgtemp=reshape(inputGF,k,z);
msg=msgtemp';
msg=gf(msg,m);
end
GF域转二进制
function data = GFdec2bin( rs_code,m )
a=size(rs_code,1);
b=size(rs_code,2);
data=zeros(a*b*m,1);
for x=1:a
for y=1:b
rs_dec=rs_code(x,y);
rs_dec=dec2bin(rs_dec.x,m);
for i=1:m
data((x-1)*b*m+(y-1)*m+i)=str2num(rs_dec(i));
end
end
end
end
GF域幂次索引
function GF_index = GFindex(M)
GF_index=zeros(2^M,1);
GF_index(1)=-1;
GF_index(2)=0;
alpha=gf(2,M);
for i=1:2^M-2
a=alpha.^i;
GF_index(a.x+1)=i;
end
end
仿真结果
信噪比大于4的场合误码率比起无编码就显著下降了,从6开始其实是0。
感想:该复习数电了。