文档采用MATLAB发布,仿真没有跑完。
CCSDS标准的LDPC编译码仿真
本脚本完成了CCSDS标准(o1)版本中适用于深空通信任务的LDPC编译码过程的仿真, 同时给出了在信息位长度为1024,码率为1/2时的误码率仿真结果
撰写人:*** 最后修改日期:2015-03-07 软件版本:MATLAB(R) 2014a
和之前版本(LDPC_Simulation.m)对比,本程序添加功能包括
- 完成了校验矩阵构造(ccsdscheckmatrix.m)
- 通过校验矩阵构造(准循环)生成矩阵 (ccsdsgeneratematrix.m)
- 按照ccsds所述,对编码后数据进行截取,完成符合ccsds文档的ldpc编译码仿真
程序修改包括
- 之前文件夹内的译码函数均不能对全0项进行正确的处理
- 问题将得到合适的修正,同时优化了输出显示
Contents
初始设置
清空工作区,数据,关闭所有窗口
clc;clear all;close all;
仿真参数设置
EbN0_dB = 1; % EbN0_dB 误码率范围,向量; FRAMES_NUM = 1; % FRAMES_NUM 最大仿真帧数目; MAX_ITER_NUM = 30; % MAX_ITER_NUM 算法最大迭代次数; MAX_ERROR_FRAME = 200; % MAX_ERROR_FRAME 仿真最大错误帧数目; bitError = zeros(1,length(EbN0_dB)); % bitError 误码数目统计,向量; BER = bitError; % BER 误码率,向量; frameError = bitError; % frameError 误帧数目; iterNumTotal = zeros(1,length(EbN0_dB)); % iterNumTotal 总迭代次数统计,向量; INFO_LENGTH = 1024; % INFO_LENGTH 信息位长度 RATE = 1/2; % RATE 码率 SIZE_M = 512; % M矩阵大小(参考CCSDS文档)
打开(新建)文件,保存运行数据
FILE_NAME = ['LDPC_CCSDS_' datestr(now,'yyyymmdd') '.txt']; fid = fopen(FILE_NAME,'at+'); fprintf(fid,'日期 %s ',datestr(now,'yyyymmdd')); fprintf(fid,'%s ','CCSDS的LDPC标准中section3中的码,06年版'); fprintf(fid,'信息位长度 = %d, ',INFO_LENGTH); fprintf(fid,'码率 = %d, ',RATE);
H为校验矩阵;G为生成矩阵。通过调用函数实现
H = ccsdscheckmatrix(SIZE_M,RATE); G = ccsdsgeneratematrix(H,SIZE_M,RATE);
译码准备
[r_mark,c_mark] = find(H~=0); HColNum = sum(H); HRowNum = cell(1,size(H,1)); for rowH = 1:size(H,1) HRowNum{rowH} = find(r_mark==rowH); end
误码率仿真
在不同EbN0下,对误码率进行仿真
for nEbN0 = 1:length(EbN0_dB)
多帧取统计平均
for nF=1:FRAMES_NUM
编码过程
message = randi([0 1],1,INFO_LENGTH); encodeData = mod(message*G,2);
调制
transmitSignal = 2*encodeData - 1; %映射 0—-1; 1—+1 transmitSignalPower = sqrt(var(transmitSignal)); transmitSignal = transmitSignal/transmitSignalPower;%归一化
AWGN 信道,SNR和EbN0换算关系
SNR_dB = EbN0_dB((nEbN0)) + 10*log10(2)+10*log10(RATE);
SNR = 10^(SNR_dB/10);
noise = randn(1,length(transmitSignal));
noise = noise/sqrt(SNR); %SNR换算,加噪声
receiveSignal = transmitSignal + noise;
根据CCSDS文档实际上有M个符号没有被传送,此处取其值为0,做仿真上的等效
receiveSignal(end-1*SIZE_M+1:end-0*SIZE_M) = 0;
译码,可调用不同的函数
[iterNum,recoverData] = ... ldpcdecoderllr(H,HRowNum,HColNum,receiveSignal,SNR,MAX_ITER_NUM); %ldpcdecoderbp1(H,HRowNum,HColNum,receiveSignal,SNR,MAX_ITER_NUM); %ldpcdecoderminsum(H,HRowNum,HColNum,receiveSignal,SNR,MAX_ITER_NUM); % 文件输出 if(nEbN0==1 && nF==1) fprintf(fid,'译码程序为ldpcdecoderllr '); fprintf(fid,'最大迭代次数 = %d, ',MAX_ITER_NUM); fprintf(fid,' ------------------------------- '); fprintf(fid,'EbN0 总帧数 误帧数 误比特数 平均迭代次数 误比特率 误帧率 '); % MATLAB的readtable笨笨的不认识中文,所以打一行给他看 fprintf(fid,'EbN0 T_Frame E_Frames E_Bits A_IterNums BER FER '); end
误比特数,误帧数,迭代次数统计
bitError(nEbN0) = bitError(nEbN0) + ... sum(abs(message - recoverData(1:length(message)))); frameError(nEbN0) = frameError(nEbN0) + ... (sum(abs(encodeData - recoverData(1:length(encodeData))))~=0); iterNumTotal(nEbN0) = iterNumTotal(nEbN0) + iterNum;
停止条件和文件输出
if ( frameError(nEbN0)>=MAX_ERROR_FRAME || nF==FRAMES_NUM) BER(nEbN0) = bitError(nEbN0)/nF/length(message); fprintf(fid,'%3.2g %5d %4d ',EbN0_dB(nEbN0),nF,frameError(nEbN0)); fprintf(fid,'%8d %8.6g ',bitError(nEbN0),iterNumTotal(nEbN0)/nF); fprintf(fid,'%e %e ',BER(nEbN0),frameError(nEbN0)/nF); break; end if (mod(nf,100)==0) BER(nEbN0) = bitError(nEbN0)/nF/320; fprintf(fid,' ------------------------------- '); fprintf('Eb/No = %e ',EbN0_dB(nEbN0)); fprintf('总帧数 = %d, ',nF); fprintf('误帧数 = %d, ',frameError(nEbN0)); fprintf('误比特数 = %d, ',bitError(nEbN0)); fprintf('最大迭代次数 = %d, ',MAX_ITER_NUM); fprintf('平均迭代次数 = %e, ',iterNumTotal(nEbN0)/nF); fprintf('误码率 = %g, ',BER(nEbN0)); fprintf('误帧率 = %g ',frameError(nEbN0)/nF); end
end
end
仿真结果
通过编码后实际误码率对比图,不要看下面的图,图没有仿真出来的,太慢了
fclose(fid); readtable(FILE_NAME,'HeaderLines',6,'Delimiter',' ') semilogy(EbN0_dB,BER,'r'); xlabel('Eb/N0(dB)'); ylabel('BER'); grid on;
ans = EbN0 T_Frame E_Frames E_Bits A_IterNums BER FER ____ _______ ________ ______ __________ ___ ___ 1 1 0 0 27 0 0