RS码编译matlab仿真2

君一一

已于 2022-08-15 10:11:45 修改

阅读量4.4k

点赞数 18

分类专栏：信道编码文章标签： matlab

于 2022-08-15 10:09:15 首次发布

本文链接：https://blog.csdn.net/qq_39731597/article/details/126341110

版权

信道编码专栏收录该内容

4 篇文章 3 订阅

订阅专栏

虽然跑出来了，但没有完全跑出来.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