image thumbnail

Channel Noise Estimation Using Particle based Belief Propagation for LDPC decoding in AWGN and BSC

by

 

Channel Noise Estimation Using Particle based Belief Propagation for LDPC decoding in AWGN and BSC

demo001_PBP_AWGN_BSC.m
% Authors:
% 
% Shuang Wang, Lijuan Cui and Samuel Cheng
% 
% Copyright:
% 
% Copyright ? 2007-2011 Shuang Wang
% 
% Questions:
% 
% For any questions, please contact me by shuangwang AT ou dot edu
%
% reference
%  L. Cui, S. Wang, S. Cheng, M. Yeary, "Adaptive Binary Slepian-Wolf Decoding using Particle Based Belief Propagation", Communications, IEEE Transactions on, to appear.
%  S. Wang, L. Cui, S. Cheng, Y. Zhai, M. Yeary, Q. Wu, "Noise Adaptive LDPC Decoding Using Particle Filtering," Communications, IEEE Transactions on, Vol 59. pp. 913 - 916, April 2011.
%

warning off all; clearvars; % clear java; 
clear import; warning on; %#ok<WNON>   % clear java;
clear import; warning off all;
addpath(genpath('.'));
p = mfilename('fullpath'); p = p(1:numel(p)-numel(mfilename)); cd(p);
javarmpath('../JavaBPForMatlab/JavaBPForMatlab.jar');
javarmpath('./JavaPBPForMatlab.jar');
warning on;  %#ok<WNON>

% if( matlabpool('size') == 0)
%     matlabpool;
% end

pegFilePath = '../JavaBPForMatlab/PCHK/';
initializeJava('../JavaBPForMatlab/JavaBPForMatlab.jar');
initializeJava('./JavaPBPForMatlab.jar');

%% initialize Arrays
particleValues = [];
trueSampleNoise = [];
mseEstimateVsSample = [];
mseEstimateVsTrue = [];
mseInitialVsTrue = [];
decodingError = [];

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% set peg File name here. for differen channel noise or differet profile, one should use different peg File, which can be generated by peg.exe.
pegFileName = 'Reg3_6_10000.dat';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Automatically load PCHK and other parameter settings for BP. user does not need to change the following settings.
randSeed = 1; % set random seed from here.
RandStream.setDefaultStream(RandStream('mt19937ar','seed',randSeed)); % this seed for generate random noise
fprintf('Initializing and loading pchk data\n');
[m, n, pchk, pchk_java] = loadPCHK(strcat(pegFilePath, pegFileName));
codeRate = m/n;
numThreads = feature('numCores'); % automatically detect number of cpu cores; 12;

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% parameters for BP
max_num_codeword = 10;
max_BP_iters = 260;
%%%%%%%%% %%%%%%  %%%%%%  %%%%%%% %%%%%%  %%%%%%%%%% %%%%%  %%%%%%%  %%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%%  %%%%%%%%  %%%%%%%%%%%% 

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% noise channel parameters
channelType_java = setChannelType(javabpformatlab.ChannelType.AWGN); % one can only use AWGN or BSC
% channelType_java = setChannelType(javabpformatlab.ChannelType.BSC); % one can only use AWGN or BSC
if(channelType_java == javabpformatlab.ChannelType.AWGN)
    trueMeanValue = 1.3; % Eb/N0 for AWGN Channel; crossover probability for BSC
    timeVaringSigma = 0.1; % timeVaringSigma = 0; for constant channel Noise
    noise_mismatch = 0;  % 
    fprintf('Using AWGN Channel; ');
elseif(channelType_java == javabpformatlab.ChannelType.BSC)
    trueMeanValue = 0.08; % Eb/N0 for AWGN Channel; crossover probability for BSC
    timeVaringSigma = 0.05; % timeVaringSigma = 0; for constant channel Noise
    noise_mismatch = 0.0;  % 
    fprintf('Using BSC; ');
end
inputTrueChannelNoise = false;  % if true, the channel input will be equal to the true channel Noise
%%%%%%%%% %%%%%%  %%%%%%  %%%%%%% %%%%%%  %%%%%%%%%% %%%%%  %%%%%%%  %%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%%  %%%%%%%%  %%%%%%%%%%%%  

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% parameters for Parameter Estimation using PBP or NONE
useParameterEstimation = setParameterEstimator('PBP');  % the options can be PBP, and NONE, where NONE means no estimator will be used.
connectionRatio = 50;
estimationIterStart = 100;
estimationIterStep = 20;
estimatorLength = floor((n - 1)/connectionRatio) + 1;
fprintf('Using Parameter Esimation: %s; \n', useParameterEstimation);
%%%%%%%%% %%%%%%  %%%%%%  %%%%%%% %%%%%%  %%%%%%%%%% %%%%%  %%%%%%%  %%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%%  %%%%%%%%  %%%%%%%%%%%%  

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% parameters for PBP estimator
if(channelType_java == javabpformatlab.ChannelType.AWGN)
    randomWalkSigma = 0.01; % for random Walk in MH
    numParticles = 12;      % number of particles for
    maxPBP_iter = 1000;     % for MH;
elseif(channelType_java == javabpformatlab.ChannelType.BSC)
    randomWalkSigma = 0.005; % for random Walk in MH
    numParticles = 12;      % number of particles for
    maxPBP_iter = 1000;     % for MH;
end

%%%%%%%%% %%%%%%  %%%%%%  %%%%%%% %%%%%%  %%%%%%%%%% %%%%%  %%%%%%%  %%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%%  %%%%%%%%  %%%%%%%%%%%%  
%%%%%%%%% %%%%%%  %%%%%%  %%%%%%% %%%%%%  %%%%%%%%%% %%%%%  %%%%%%%  %%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%%  %%%%%%%%  %%%%%%%%%%%%  
%%%%%%%%% %%%%%%  %%%%%%  %%%%%%% %%%%%%  %%%%%%%%%% %%%%%  %%%%%%%  %%%%%%%% %%%%%%%%%%%%%%%%%% %%%%%%%%%% %%%%%%%%%%%%  %%%%%%%%  %%%%%%%%%%%%  

%% generate channel Noise
[trueChannelNoise, inputChannelNoise] =generateChannelNoise(channelType_java, estimatorLength, trueMeanValue, timeVaringSigma, noise_mismatch, codeRate, inputTrueChannelNoise);

%% generate connection Map
[connectionMap, connectionMap_java] = getConnectionMap(numThreads, n, connectionRatio, numParticles);

%% create BP decoder from Java side
fprintf('Start Decoding\n');
bpDecoder_java = javabpformatlab.BPDecoder(numThreads, pchk_java, channelType_java);

%% factory injection; select channel noise model
channelNoiseModel_java = getChannelNoiseModel(channelType_java);

%% create pbp estimator from java
if(strcmp(useParameterEstimation, 'PBP'))
    [priorDistributionParams_java, messageFromPrior_java] = getPriorDistributionParams(inputChannelNoise, channelType_java);
    pbpEstimator_java = javapbpformatlab.PBPEstimator(randomWalkSigma, randSeed, maxPBP_iter, n, connectionMap_java, estimatorLength, numParticles, numThreads, channelNoiseModel_java);        
end
%% perform decoding for many different codewords 
randSeed = randSeed + 10; % reset random seed from here, so that each codeword are generated from the same rand sequence.
RandStream.setDefaultStream(RandStream('mt19937ar','seed',randSeed));
for iter_codeword = 1: max_num_codeword       
    fprintf('%05d/%d  ',iter_codeword, max_num_codeword);
    %% generate source according to the channelType and channelNoise
    [source_x, source_y] =generateSource(channelType_java, trueChannelNoise, connectionRatio, n);
    %% initialize bp Decoder
    bpDecoder_java.initialization(source_y, Arrays.expandArray(inputChannelNoise, connectionRatio, n));    
    %% initialize PBP estimator
    if(strcmp(useParameterEstimation, 'PBP'))
        % initialize pbp estimator
        pbpEstimator_java.initialization(repmat(inputChannelNoise, 1, numParticles), priorDistributionParams_java, source_y); 
    end
    %%
    h = progressBar([], [], max_BP_iters);
    estimation_count = 0;
    %% start BP decoding
    for iter_BP = 1:max_BP_iters
        h = progressBar(h, iter_BP);
        %% BP decoding
        bpDecoder_java.decode(iter_BP, bpDecoder_java);
        decodingError = sum(bitxor(bpDecoder_java.decodedWord, source_x)); 
        %% perform PBP estimation is possible
        if(~strcmp(useParameterEstimation, 'NONE') && iter_BP >= estimationIterStart && (mod(iter_BP, estimationIterStep) == 0 || iter_BP == max_BP_iters))
            if(strcmp(useParameterEstimation, 'PBP'))
                estimation_count = estimation_count + 1;
                pbpEstimator_java.estimate(bpDecoder_java.bprbWithoutPrior);
                bpDecoder_java.updateLratio(pbpEstimator_java.getLratio());           
            end
        end  
        %% check decoding errors or max_BP_iters
        if((decodingError == 0) || iter_BP == max_BP_iters)  
            %% finish the final EP/BP estimation and get estimation result
             if(~strcmp(useParameterEstimation, 'NONE'))
                 if(estimation_count == 0) % when decoding is successful, but iter_bp < estimationIterStart, we perform estimation here.
                     estimation_count = estimation_count + 1;
                     if(strcmp(useParameterEstimation, 'PBP'))
                        pbpEstimator_java.estimate(bpDecoder_java.bprbWithoutPrior);                     
                     end
                 end  
                 if(strcmp(useParameterEstimation, 'PBP'))
                    [particleValues(iter_codeword,:), trueSampleNoise(iter_codeword,:),...
                    mseEstimateVsSample(iter_codeword), mseEstimateVsTrue(iter_codeword),...
                    mseInitialVsTrue(iter_codeword)] = getEstimationResultsPBP(pbpEstimator_java, channelType_java, trueChannelNoise, inputChannelNoise, source_y, source_x, connectionRatio); %#ok<SAGROW>                 
                 end
             end;
            %% calculate decoding errors 
            h = progressBar(h, max_BP_iters);
            decodingError_codeword(iter_codeword) = decodingError; %#ok<SAGROW>
            iter_BP_codeword(iter_codeword) = iter_BP; %#ok<SAGROW>
            fprintf('err: %5d; iter: %3d; Ave err: %e; Ave Iter %4.2f; mseEstimateVsSample: %4.5f, mseEstimateVsTrue: %4.5f, mseInitialVsTrue: %4.5f;\n',...
                    decodingError, iter_BP, mean(decodingError_codeword)/n, round(mean(iter_BP_codeword)), mean(mseEstimateVsSample), mean(mseEstimateVsTrue), mean(mseInitialVsTrue));
            break;
        end
    end
end

plotEstimation();

Contact us