diff --git a/ChadCole_LDPC_ErrorFloor_Archive/H2048_8023an.mat b/ChadCole_LDPC_ErrorFloor_Archive/H2048_8023an.mat new file mode 100644 index 0000000..7feab3c Binary files /dev/null and b/ChadCole_LDPC_ErrorFloor_Archive/H2048_8023an.mat differ diff --git a/ChadCole_LDPC_ErrorFloor_Archive/ISldpc.m b/ChadCole_LDPC_ErrorFloor_Archive/ISldpc.m new file mode 100644 index 0000000..9e8b218 --- /dev/null +++ b/ChadCole_LDPC_ErrorFloor_Archive/ISldpc.m @@ -0,0 +1,434 @@ +%%%%%%%%%%%%% +%% ISldpc_ms_new_doc.m - Requires input H_sparse matrix. +%% +%% This program implements the last step of the 3-step program outlined +%% in "A General Method for Finding Low Error Rates of LDPC Codes." The +%% last step uses the technique of importance sampling (IS) to get a +%% better estimate than Monte Carlo (MC) can provide of the FER/BER at +%% higher SNR's. +%% +%% The core decoder implementation +%% is the same as in the program 'irregularSPA.m,' another piece of +%% software in this suite. +%%%%%%%%%%%% + +clear all +randn('state',sum(100*clock)); +format compact +%function ISldpc_ms_new_doc() + +% load H1008_nox2 +% load H1200_4_8_no6cycle +% load H2048_8023an +% load H504_no82 +H = [1 1 1 1 0 0 0 0 0 0; + 1 0 0 0 1 1 1 0 0 0; + 0 1 0 0 1 0 0 1 1 0; + 0 0 1 0 0 1 0 1 0 1; + 0 0 0 1 0 0 1 0 1 1] +H_sparse = sparse(H); + + +n = length(H_sparse(1,:)); +n_k = length(H_sparse(:,1)); +k = n - n_k; +rate = k/n; + +H_col_sum = sum(H_sparse,1); +H_row_sum = sum(H_sparse,2); +%B will be sorted in ascending order +B = unique(H_col_sum); +num_v_deg = length(B); +R = unique(H_row_sum); +num_c_deg = length(R); +deg_c_prof = zeros(length(R), 2); +deg_v_prof = zeros(length(B), 2); + +%%%% +% Set up degree profile info for V-nodes & C-nodes and make sure H is +% ordered from highest to lowest degree in columns and rows. +% +% Example - n-k = 500, half parities have deg 6, half 8 +% always list deg profiles in DESCENDING ORDER +% : deg_c_prof = [250, 8; +% 250, 6] +H_new = spalloc(n-k,n,length(find(H_sparse))); +running_index = 0; +for ii = 1:num_v_deg + change_loc = find(H_col_sum == B(num_v_deg + 1 - ii)); + deg_v_prof(ii, 1) = length(change_loc); + deg_v_prof(ii, 2) = B(num_v_deg + 1 - ii); +%This step orders columns from highest degree to lowest + for jj = 1:length(change_loc) + running_index = running_index + 1; + H_new(:,running_index) = H_sparse(:, change_loc(jj)); + end +end +H_sparse = H_new; +clear H_new; + +H_new = spalloc(n-k,n,length(find(H_sparse))); +running_index = 0; +for ii = 1:num_c_deg + change_loc = find(H_row_sum == R(num_c_deg + 1 - ii)); + deg_c_prof(ii, 1) = length(change_loc); + deg_c_prof(ii, 2) = R(num_c_deg + 1 - ii); +%This step orders rows from highest degree to lowest + for jj = 1:length(change_loc) + running_index = running_index + 1; + H_new(running_index,:) = H_sparse(change_loc(jj),:); + end +end +H_sparse = H_new; +clear H_new; + +Vlist=zeros(n-k,max(sum(H_sparse,2))+1); % list of V-nodes +Clist=zeros(n,max(sum(H_sparse,1))+1); % list of C-nodes +for jj=1:n-k + Vlist(:,1)=sum(H_sparse,2); + icnt=0; + for ii=1:n + if H_sparse(jj,ii)==1 + icnt=icnt+1; + Vlist(jj,icnt+1)=ii; + end + end +end + +for ii=1:n + Clist(:,1)=(sum(H_sparse,1))'; + jcnt=0; + for jj=1:n-k + if H_sparse(jj,ii)==1 + jcnt=jcnt+1; + Clist(ii,jcnt+1)=jj; + end + end +end + +num_edges_v = sum(Clist(:,1)); +num_edges_c = sum(Vlist(:,1)); +if (num_edges_c ~= num_edges_v) + text = ['error - row/col mismatch in H'] +end +num_edges = num_edges_c; +clear num_edges_c num_edges_v + +Lr_ind = zeros(1, num_edges); +Lq_ind = zeros(1, num_edges); + +tot_count = 0; +for kk = 2:length(Clist(1,:)) + for jj=1:n + if(Clist(jj,kk) ~= 0) + tot_count = tot_count + 1; + Lq_ind(tot_count) = (kk-2)*n + jj; + end + end +end + +row_index = zeros(n-k, 1); +tot_count = 0; +for kk=2:length(Clist(1,:)) + for jj=1:n + if(Clist(jj,kk) ~= 0) + tot_count = tot_count + 1; + row=Clist(jj, kk); + row_index(row) = row_index(row) + 1; + Lr_ind(tot_count) = row + (row_index(row)-1)*(n-k); + end + end +end + +Lq=zeros(n,max(Clist(:,1))); % messages up +Lr=zeros(n-k,max(Vlist(:,1))); % messages down +LQ=zeros(1,n); +Lc=zeros(1,n); + + +small_num = eps; +%%% Set SNR's for a desired IS performance estimate. +% ebnodb = [4]; +% ebnodb = [4 4.5 5 5.5 6]; +ebnodb = 4:1:24; + +%% The following variables collect some information which can help +%% determine how accurate the IS estimate is. +num_new_error_event = 0; +%% The error_event variable keeps track of new errors that occur throughout +%% the simulation that were not included in the initial list obtained from +%% the TSearch.m program. +error_event = spalloc(1000000,n,10000*30); +%% This counts how many occurrences of each new error events occurs. +error_event_count = zeros(1,1000); + +itenum=50; % maximum # of iterations + + +% load d_eH504_no82_del36_1 +% load H504_no82_del36_op8_ebno5_TS +% load H2048_8023an_cw_TS +load H10_test + +%% Only use the TS with d_e^2 < some threshold +% TS = TS(find(d_e_2 < 20),:); + +if (exist('TS') == 0) || isempty(TS) + num_TS = 0; + TS = []; +else + num_TS = length(TS(:,1)); + temp = zeros(num_TS,2); + temp(:,1) = sum(TS(1:num_TS,:),2); + temp(:,2) = sum(mod(TS(1:num_TS,:)*H_sparse',2),2); +end + +if exist('cw') == 0 + cw = []; +end + + +%% Another way to pick certain TS based on their (a,b) value, for example +%% in the following we just use (10,2) TS. +% TS = TS(intersect(find(temp(:,1)==10), find(temp(:,2)==2)),:); + +%% mu is the list of mean-shift candidates in the f^* used in IS +mu = [TS; cw]; +mu = sparse(mu); +num_mu = length(mu(:,1)); + +%% Keep track of any new codewords found throughout simulation. +num_miscor = 0; +if isempty(cw) + num_cw = 0; +else + num_cw = length(cw(:,1)); +end + +%% mu_coef specifies the amount we mean shift in each bit of a TS. The +%% optimal value is not `1' unless we're shifting towards a valid codeword, +%% but for most TS, `1' will suffice. +mu_coef = 1.0; + +%% We typically generate noise realizations shifted towards each of the mu +%% bias points an equal amount of time, for example 5000 in this case. +num_trials = 400*num_mu*ones(1,length(ebnodb)); + +%% The FER/BER estimates for each SNR. +Pf_vec = zeros(1, length(ebnodb)); +Pb_vec = zeros(1, length(ebnodb)); + +%% The hit_rate counts the number of `hits' (errors) for each shift point at each SNR. +%% We can use this info as an indicator of how good a certain shift point +%% is. If no hits were made in 5000 realizations of noise shifted towards +%% a certain TS, then that TS is probably not a significant contributor to +%% the error rate. `intended hits' occur when the error is the exact same one +%% as that which we shifted towards. These are the most common types of +%% hits and as SNR increases, the number of `intended hits' should approach +%% the number of total hits. +hit_rate = zeros(length(ebnodb), num_mu); +hit_rate_intended = zeros(length(ebnodb), num_mu); + +%% ite_hist keeps a history of the number of iterations needed for each +%% decoding. +ite_hist = zeros(1,num_trials(1)); +%% ite_hist_mean contains average number of iterations required at each SNR +ite_hist_mean = zeros(1, length(ebnodb)); + +for ebnoloop = 1:length(ebnodb) + ebno = 10^((ebnodb(ebnoloop))/10) + esno = ebno*(rate);%*3; %rate 1/2, 3 bits/sym + sigma_2 = 1/(2*esno); + sigma = sqrt(sigma_2); + Lc_coef = 4*esno; + a = ones(1,n); + variance = 0; + num_frame_errors = 0; %don't forget to reset these to 0 for each ebno + num_bit_errors = 0; + t = 0; + coef = -1/(2*sigma_2); + %% s_f is the scale factor necessary to prevent exp(M) evaluating to + %% zero (in Matlab) when we calculate the weight function. + s_f = (n/2); + +while (t < num_trials(ebnoloop)) + + t = t + 1; + %% For easy encoding, send all 0's (s0) message. + %% The initial point in decoding space is biased towards one of the + %% vectors in the mu matrix. + y = a - mu_coef*mu(mod(t,num_mu)+1,:) + sigma*randn(1,n); + %% This is different than the typical Monte Carlo simulation which would + %% send instead: + %% y = a + sigma*randn(1,n); + +% MPA +% step 1: initialization + +Lc=Lc_coef*y; +num_total = 0; +for ii = 1:num_v_deg + Lq(num_total+1:num_total+deg_v_prof(ii, 1), 1:deg_v_prof(ii, 2)) = Lc(num_total+1:num_total + deg_v_prof(ii, 1))'*ones(1,deg_v_prof(ii, 2)); + num_total = num_total + deg_v_prof(ii, 1); +end + + +stopsig=0; % check whether a valid codeword has been found +itestep=0; + +while (stopsig==0) && (itestep= deg_current) + num_combos = 1; + indices = [1:deg_current]'; + else + indices = zeros(v_num, nchoosek(deg_v_prof(length(deg_v_prof(:,1))+1-iii,2), v_num)); + for pp = 1:2^deg_current + binary = int2vec(pp, deg_current); + if (sum(binary) == v_num) + num_combos = num_combos + 1; + indices(:, num_combos) = find(binary); + end + end + end + num_combos_save(iii) = num_combos; + + for jj = 1:deg_v_prof(length(deg_v_prof(:,1))+1-iii,1) + v_count = v_count - 1 + for kk = 1:num_combos + +% The c_nodes structure contains the check-nodes involved in the tree rooted +% at the v_count^{th} variable node + c_nodes = Clist(v_count, indices(:,kk)+1); +% The v_tier_1 structure contains all possible variable nodes that are in +% variable tier 1 of the tree rooted at the v_count^{th} variable node. + v_tier_1 = zeros(length(c_nodes), max(Vlist(c_nodes, 1))); + for ppp = 1:length(c_nodes) + v_tier_1(ppp, 1:Vlist(c_nodes(ppp),1)-1) = setxor(Vlist(c_nodes(ppp),2:Vlist(c_nodes(ppp),1)+1),v_count); + end + % v_tier_1_ind contains indices of nonzero elements of v_tier_1 + v_tier_1_ind = find(v_tier_1); + %store candidates for epsilon_2 impulses in this matrix which only has + %to be calculated once for each of our n trees + epsilon_mat = zeros(length(v_tier_1_ind), (max(Vlist(:,1))-1)*(max(Clist(v_tier_1(v_tier_1_ind),1))-1)); + for ppp = 1:length(v_tier_1_ind) + tot = 0; + for qqq = 1:Clist(v_tier_1(v_tier_1_ind(ppp)),1) + if sum(Clist(v_tier_1(v_tier_1_ind(ppp)),qqq+1)*ones(1,length(c_nodes)) == c_nodes)==0 + epsilon_mat(ppp, tot+1:tot+Vlist(Clist(v_tier_1(v_tier_1_ind(ppp)),qqq+1),1)-1) ... + = setxor(Vlist(Clist(v_tier_1(v_tier_1_ind(ppp)),qqq+1), 2: ... + Vlist(Clist(v_tier_1(v_tier_1_ind(ppp)),qqq+1),1)+1), ... + v_tier_1(v_tier_1_ind(ppp))); + tot = tot + Vlist(Clist(v_tier_1(v_tier_1_ind(ppp)),qqq+1),1)-1; + end + end + end + + + %tot_perms counts the number of permutations of impulse bits for each + % tree + tot_perms = 1; + for ppp = 1:length(c_nodes) + tot_perms = tot_perms*(Vlist(c_nodes(ppp),1)-1); + end + + %bit_combos enumerates all tot_perms length(c_nodes)-bit error impulse + % candidate sets + bit_combos = zeros(length(c_nodes),tot_perms); + count = zeros(length(c_nodes),1); + for ppp = 1:tot_perms + v_prod = 1; + for zzz = 1:length(c_nodes) + if (mod(ppp, v_prod)==0) + count(zzz) = mod(count(zzz)+1,Vlist(c_nodes(zzz),1)-1); + end + v_prod = v_prod*(Vlist(c_nodes(zzz),1)-1); + end + for zzz = 1:length(indices(:,1)) + bit_combos(zzz, ppp) = v_tier_1(zzz, count(zzz)+1); + end + end + + for ll = 1:tot_perms + +%%%%The following code can be +%%%%uncommented for finding TS in codes where searching the first layer +%%%%of variable nodes falls short, typically for larger variable degree +%%%%nodes, such as {4,8} codes with n > 2000 or so. Don't forget to also +%%%%uncomment the end of the for loop below this block of code. + +%%%% Find c_nodes connected to leftmost v_node at tier 1. + +% v_root=intersect(Vlist(c_nodes(1),2:Vlist(c_nodes(1),1)+1), ... +% bit_combos(:,ll)); +% v_root_cnodes = setxor(Clist(v_root,2:Clist(v_root,1)+1), c_nodes(1)); +% v_tier_2 = zeros(length(v_root_cnodes), max(Vlist(v_root_cnodes, 1))); +% for ppp = 1:length(v_root_cnodes) +% v_tier_2(ppp, 1:Vlist(v_root_cnodes(ppp),1)-1) = setxor(Vlist(v_root_cnodes(ppp),2:Vlist(v_root_cnodes(ppp),1)+1),v_root); +% end +% +% tot_perms_2 = 1; +% for ppp = 1:length(v_root_cnodes) +% tot_perms_2 = tot_perms_2*(Vlist(v_root_cnodes(ppp),1)-1); +% end +% +% bit_combos_2 = zeros(length(v_root_cnodes),tot_perms_2); +% count_2 = zeros(length(v_root_cnodes),1); +% for ppp = 1:tot_perms_2 +% v_prod = 1; +% for zzz = 1:length(v_root_cnodes) +% if (mod(ppp, v_prod)==0) +% count_2(zzz) = mod(count_2(zzz)+1,Vlist(v_root_cnodes(zzz),1)-1); +% end +% v_prod = v_prod*(Vlist(v_root_cnodes(zzz),1)-1); +% end +% for zzz = 1:length(v_root_cnodes) +% bit_combos_2(zzz, ppp) = v_tier_2(zzz, count_2(zzz)+1); +% end +% end +% +% for mm = 1:tot_perms_2 +% + num_checks = num_checks + 1; + d = zeros(1,n); + d(v_count) = 1; + d(bit_combos(:, ll)) = 1; +% d(bit_combos_2(:, mm)) = 1; + v_node_hit_hist(bit_combos(:, ll)) = v_node_hit_hist(bit_combos(:, ll)) + 1; + v_node_hit_hist(v_count) = v_node_hit_hist(v_count) + 1; + syn_sum = sum(mod(d*H_sparse',2)); + if (syn_sum <= 14) + syn_save_hist(syn_sum+1) = syn_save_hist(syn_sum+1) + 1; + + y = gamma(iii)*a - d*epsilon_1; + %Only want nonzero values of epsilon_mat to index into y, + % this not problem for (check) regular codes. also a problem + % here if girth < 8 + [garbage,garbage2,useful] = intersect(bit_combos(:, ll),v_tier_1(v_tier_1_ind)); +% y(epsilon_mat(useful,:)) = epsilon_2; + +% MPA +% step 1: initialization + +Lc=Lc_coef*y; + +num_total = 0; +for ii = 1:num_v_deg + Lq(num_total+1:num_total+deg_v_prof(ii, 1), 1:deg_v_prof(ii, 2)) = Lc(num_total+1:num_total + deg_v_prof(ii, 1))'*ones(1,deg_v_prof(ii, 2)); + num_total = num_total + deg_v_prof(ii, 1); +end + +stopsig=0; % check whether a valid codeword has been found +%itenum=20; % maximum # of iterations +itestep=0; + +while (stopsig==0) && (itestep val) + while (lll < num_TS) && not_done_outer + lll = lll + 1; + if (sum(TS(lll,:) == TS_temp) == n) + not_done_outer = 0; + TS_unique_count(lll) = TS_unique_count(lll) + 1; + end + end + if (not_done_outer) %add new TS + num_TS = num_TS + 1 + TS(num_TS,:) = TS_temp; + TS_unique_count(num_TS) = 1; + end + end +end + +end + +end %if we decode + +% weight = 1; +% num_frame_errors = num_frame_errors + weight; +% variance = variance + weight^2; + +% end % tot_perms_2 loop - uncomment along with above stuff +end %tot_perms loop + + + end % num_combos loop + end %v_count loop +end % num deg prof while +toc + +TS = TS(1:num_TS,:); +cw = cw(1:num_cw,:); +%save Hfilename_TS TS cw diff --git a/ChadCole_LDPC_ErrorFloor_Archive/determ_boundary_finder.m b/ChadCole_LDPC_ErrorFloor_Archive/determ_boundary_finder.m new file mode 100644 index 0000000..fd18a64 --- /dev/null +++ b/ChadCole_LDPC_ErrorFloor_Archive/determ_boundary_finder.m @@ -0,0 +1,346 @@ +%%%%%%%%% +% This file contains the implementation of Step-2 of the 3-Step error floor +% analysis procedure. The input is an H matrix in sparse format called +% H_sparse. The output variable is d_e_2 which should be saved and used +% in step 3, which is implemented in ISldpc.m. + +clear all + +% load H4000_wesel_test + +% load H504_4_8_no6_fallback4_best +% load test_1200H_3_6_no6 + +% load H1000_4_8_gallager +% load H1008_nox2 +% load QR_H8192 +% load Margulis11qHCV +% load H1200_4_8_no6cycle +% load H504goodcode_no82 +% load H1008_peg +load H96_firstmackay + +n = length(H_sparse(1,:)); +n_k = length(H_sparse(:,1)); +k = n - n_k; +rate = k/n; + +%%%%%% +% These are MPA messages: +% Lq = Variable->Check +% Lr = Check->Variable +% LQ = Marginal LLR used in making a hard decision. +Lq=zeros(n,max(sum(H_sparse,1))); +Lr=zeros(n-k,max(sum(H_sparse,2))); +LQ=zeros(1,n); + +%%%% Organize H columns from highest degree -> lowest. this is necessary +%%%% for the decoding algorithm to make use of Matlab vectorized commands +%%%% which significantly speed up program run-time. +H_col_sum = sum(H_sparse,1); +H_row_sum = sum(H_sparse,2); +%B will be sorted in ascending order +B = unique(H_col_sum); +num_v_deg = length(B); +R = unique(H_row_sum); +num_c_deg = length(R); +deg_c_prof = zeros(length(R), 2); +deg_v_prof = zeros(length(B), 2); + +%%%% +% Set up degree profile info for V-nodes & C-nodes and make sure H is +% ordered from highest to lowest degree in columns and rows. +% +% Example - n-k = 500, half parities have deg 6, half 8 +% always list deg profiles in DESCENDING ORDER +% : deg_c_prof = [250, 8; +% 250, 6] +H_new = spalloc(n-k,n,length(find(H_sparse))); +running_index = 0; +for ii = 1:num_v_deg + change_loc = find(H_col_sum == B(num_v_deg + 1 - ii)); + deg_v_prof(ii, 1) = length(change_loc); + deg_v_prof(ii, 2) = B(num_v_deg + 1 - ii); +%This step orders columns from highest degree to lowest + for jj = 1:length(change_loc) + running_index = running_index + 1; + H_new(:,running_index) = H_sparse(:, change_loc(jj)); + end +end +H_sparse = H_new; +clear H_new; + +H_new = spalloc(n-k,n,length(find(H_sparse))); +running_index = 0; +for ii = 1:num_c_deg + change_loc = find(H_row_sum == R(num_c_deg + 1 - ii)); + deg_c_prof(ii, 1) = length(change_loc); + deg_c_prof(ii, 2) = R(num_c_deg + 1 - ii); +%This step orders rows from highest degree to lowest + for jj = 1:length(change_loc) + running_index = running_index + 1; + H_new(running_index,:) = H_sparse(change_loc(jj),:); + end +end +H_sparse = H_new; +clear H_new; + +%%%%%%% The Clist Vlist Data Structures +% Clist is an n by (max(deg_v)+1) matrix with the first column having the +% variable node degree (deg_v) of each column of H. The next deg_v entries +% for that row are the check node numbers which that variable node is +% connected to. Any extra columns following are set to zero. +% Vlist is an (n-k) by (max(deg_c)+1) matrix with the first column having the +% check node degree (deg_c) of each row of H. The next deg_c entries +% for that row are the variable node numbers which that check node is +% connected to. Any extra columns following are set to zero. + +Vlist=zeros(n-k,max(sum(H_sparse,2))+1); +Clist=zeros(n,max(sum(H_sparse,1))+1); + +for jj=1:n-k + Vlist(jj,1)=sum(H_sparse(jj,:)); + icnt=0; + for ii=1:n + if H_sparse(jj,ii)==1 + icnt=icnt+1; + Vlist(jj,icnt+1)=ii; + end + end +end + +for ii=1:n + Clist(ii,1)=sum(H_sparse(:,ii)); + jcnt=0; + for jj=1:n-k + if H_sparse(jj,ii)==1 + jcnt=jcnt+1; + Clist(ii,jcnt+1)=jj; + end + end +end + +%Maximum number of MPA iterations before we declare a failure. +itenum = 50; + +small_num = eps; %how expensive is eps()? + +% Simple H matrix consistency check. +num_edges_v = sum(Clist(:,1)); +num_edges_c = sum(Vlist(:,1)); +if (num_edges_c ~= num_edges_v) + text = ['error - row/col mismatch in H'] +end +num_edges = num_edges_c; +clear num_edges_c num_edges_v + +Lr_ind = zeros(1, num_edges); +Lq_ind = zeros(1, num_edges); + +%%%%% +% Set up the edge connections between Lr (C->V messages) and Lq (V->C). +% Remember, Matlab matrices are indexed by columns, so if Lr is of +% dimension (n-k)X max(d_c)+1, then Lr(n-k+1) is the same as Lr(1,2). +tot_count = 0; +for kk = 2:length(Clist(1,:)) + for jj=1:n + if(Clist(jj,kk) ~= 0) + tot_count = tot_count + 1; + Lq_ind(tot_count) = (kk-2)*n + jj; + end + end +end + +row_index = zeros(n-k, 1); +tot_count = 0; +for kk=2:length(Clist(1,:)) + for jj=1:n + if(Clist(jj,kk) ~= 0) + tot_count = tot_count + 1; + row=Clist(jj, kk); + row_index(row) = row_index(row) + 1; + Lr_ind(tot_count) = row + (row_index(row)-1)*(n-k); + end + end +end + + + +%load H1008_peg_del4_op6_ebno6_TS +%load H504_peg_del3_op8_ebno6_TS +%load H504_no82_del36_op8_ebno5_TS +%load H1200_4_8_no6cycles_5bitimpulse_TS +% load H1008_nox2_4delta_07op_7db_TS +% load QR_8192_4bit_TS +% load Margulis_TS_tot + +% load H504_4_8_no6_fallback4_best_TS +% load test_1200H_3_6_no6_TS +% load dinoi_TS_save +% load H204_TS +% load H1200_4_8_no6cycles_5bitimpulse_TS_extra_2 +% load H1200_4_8_no6cycles_5bitimpulse_TS +% load H504_no82_TS_temp +% load H4000_wesel_test_TS_cw +% load H1008_peg_del3_op6_ebno6_TS +load H96-964_5_1_TS + +num_TS = length(TS(:,1)); +temp = zeros(num_TS,2); +temp(:,1) = sum(TS(1:num_TS,:),2); +temp(:,2) = sum(mod(TS(1:num_TS,:)*H_sparse',2),2) + +% TS = TS(find(d_e < 30),:); +TS = TS(intersect(find(temp(:,1)==5), find(temp(:,2)==1)),:); + + +e_max = 3.5; +e_min = 1; +num_intervals = 10; + +mu_coef = 1; + +num_bits = zeros(length(ebnodb), length(TS(:,1))); + +ebnodb = 6; +% ebnodb = [3 4 5 6 7]; +a = ones(1,n); +e_vec = zeros(length(TS(:,1)), length(ebnodb), 2); + + +for outer = 1:length(TS(:,1)) + + mu = TS(outer,:); + + for ebnoloop = 1:length(ebnodb) + + upper = e_max + lower = e_min; + + ebno = 10^((ebnodb(ebnoloop))/10); + esno = ebno*(rate);%*3; %rate 1/2, 3 bits/sym + sigma_2 = 1/(2*esno); + sigma = sqrt(sigma_2); + Lc_coef = 4*esno; + + variance = 0; + num_frame_errors = 0; %don't forget to reset these to 0 for each ebno + + for e_loop = 1:num_intervals + mu_coef = lower + (upper - lower)/2; + + y = a - mu_coef*mu; + + % MPA +% step 1: initialization +Lc=Lc_coef*y; +num_total = 0; +for ii = 1:num_v_deg + Lq(num_total+1:num_total+deg_v_prof(ii, 1), 1:deg_v_prof(ii, 2)) = Lc(num_total+1:num_total + deg_v_prof(ii, 1))'*ones(1,deg_v_prof(ii, 2)); + num_total = num_total + deg_v_prof(ii, 1); +end + + +stopsig=0; % check whether a valid codeword has been found +itestep=0; + +while (stopsig==0) && (itestep= 2^i) + val = val - 2^i; + binvec(i+1) = 1; + else + binvec(i+1) = 0; + end + end +% return binvec; + + + \ No newline at end of file diff --git a/ChadCole_LDPC_ErrorFloor_Archive/irregularSPA.m b/ChadCole_LDPC_ErrorFloor_Archive/irregularSPA.m new file mode 100644 index 0000000..009576b --- /dev/null +++ b/ChadCole_LDPC_ErrorFloor_Archive/irregularSPA.m @@ -0,0 +1,385 @@ +%%%%%%%%%%%%% +% This program is a general Matlab LDPC decoder simulator for an AWGN channel. +% Both regular and irregular LDPC code performance can be simulated. The +% full belief propagation and min-sum approximation, using full double +% precision messages, are the two decoding algorithms currently +% implemented. +% +% The full belief propagation version of this code can simulate roughly 100 +% Kb/s (coded bits). Although this is not near as much throughput as most +% hardware simulations, the dominant TS making up the error floor are +% collected is this software implementation. +% +% The codes to be simulated should be in H_sparse form + + +clear all +% function irregularSPA() +randn('state',sum(100*clock)); + +%%%%% +% These would be some example files containing H_sparse +% load H1000_4_8_gallager +% load H504goodcode_no82 +% load QCH2331 +% load QCH8160 +% load ramanujanq17p5 +% load H504_4_8_no4cycle +% load H96_4_8_no4cycle +% load 128x256regular.mat H +% load H1008_nox2 +% load H96_4_8_70cw +% load H504_4_8_no6_fallback4_best +% load H26_reallysmall_no4cycle +% load test_1200H_3_6_no6 +% load H1200_4_8_no6cycle +load H2048_8023an +% load H4000_wesel_test +% load H128_highrate_test +% load H500_rate08_3_15 +% load H500_rate08_3_15_best +% load H1000_rate08_4_20 + +% H_sparse = sparse(H); + +n = length(H_sparse(1,:)); +n_k = length(H_sparse(:,1)); +k = n - n_k; +rate = k/n; + +%%%%%%%% To generate G matrix for encoding uncomment the following. +% This is generally not necessary because sending all zeros codeword ok for +% linear codes. Requires files rearrange_cols, inv_GF2 + +% H = full(H_sparse); +% H = rearrange_cols(H); +% %swap full rank cols to last part +% temp = H(:, 1:n-k); +% H(:, 1:n-k) = H(:, n-k+1:n); +% H(:, n-k+1:n) = temp; +% C2=H(:, n-k+1:n); +% %only need to do this N^3 operation once +% C2_inv = inv_GF2(C2); +% H_system = mod(C2_inv*H, 2); +% P_trans = H_system(:, 1:n-k); +% G = cat(2, eye(k), P_trans'); + +%Max number of noisy messages sent per SNR +num_trials = 20000000; + +%%%%%% +% These are MPA messages: +% Lq = Variable->Check +% Lr = Check->Variable +% LQ = Marginal LLR used in making a hard decision. +Lq=zeros(n,max(sum(H_sparse,1))); +Lr=zeros(n-k,max(sum(H_sparse,2))); +LQ=zeros(1,n); + +num_frame_errors = 0; +num_bit_err = 0; + +%SNR values in DeciBels +% ebnodb = [1.5:0.5:3.0]; +ebnodb = [3:.5:6]; + +%%%% Organize H columns from highest degree -> lowest. this is necessary +%%%% for the decoding algorithm to make use of Matlab vectorized commands +%%%% which significantly speed up program run-time. +H_col_sum = sum(H_sparse,1); +H_row_sum = sum(H_sparse,2); +%B will be sorted in ascending order +B = unique(H_col_sum); +num_v_deg = length(B); +R = unique(H_row_sum); +num_c_deg = length(R); +deg_c_prof = zeros(length(R), 2); +deg_v_prof = zeros(length(B), 2); + +%%%% +% Set up degree profile info for V-nodes & C-nodes and make sure H is +% ordered from highest to lowest degree in columns and rows. +% +% Example - n-k = 500, half parities have deg 6, half 8 +% always list deg profiles in DESCENDING ORDER +% : deg_c_prof = [250, 8; +% 250, 6] +H_new = spalloc(n-k,n,length(find(H_sparse))); +running_index = 0; +for ii = 1:num_v_deg + change_loc = find(H_col_sum == B(num_v_deg + 1 - ii)); + deg_v_prof(ii, 1) = length(change_loc); + deg_v_prof(ii, 2) = B(num_v_deg + 1 - ii); +%This step orders columns from highest degree to lowest + for jj = 1:length(change_loc) + running_index = running_index + 1; + H_new(:,running_index) = H_sparse(:, change_loc(jj)); + end +end +H_sparse = H_new; +clear H_new; + +H_new = spalloc(n-k,n,length(find(H_sparse))); +running_index = 0; +for ii = 1:num_c_deg + change_loc = find(H_row_sum == R(num_c_deg + 1 - ii)); + deg_c_prof(ii, 1) = length(change_loc); + deg_c_prof(ii, 2) = R(num_c_deg + 1 - ii); +%This step orders rows from highest degree to lowest + for jj = 1:length(change_loc) + running_index = running_index + 1; + H_new(running_index,:) = H_sparse(change_loc(jj),:); + end +end +H_sparse = H_new; +clear H_new; + +%%%%%%% The Clist Vlist Data Structures +% Clist is an n by (max(deg_v)+1) matrix with the first column having the +% variable node degree (deg_v) of each column of H. The next deg_v entries +% for that row are the check node numbers which that variable node is +% connected to. Any extra columns following are set to zero. +% Vlist is an (n-k) by (max(deg_c)+1) matrix with the first column having the +% check node degree (deg_c) of each row of H. The next deg_c entries +% for that row are the variable node numbers which that check node is +% connected to. Any extra columns following are set to zero. + +Vlist=zeros(n-k,max(sum(H_sparse,2))+1); +Clist=zeros(n,max(sum(H_sparse,1))+1); + +for jj=1:n-k + Vlist(jj,1)=sum(H_sparse(jj,:)); + icnt=0; + for ii=1:n + if H_sparse(jj,ii)==1 + icnt=icnt+1; + Vlist(jj,icnt+1)=ii; + end + end +end + +for ii=1:n + Clist(ii,1)=sum(H_sparse(:,ii)); + jcnt=0; + for jj=1:n-k + if H_sparse(jj,ii)==1 + jcnt=jcnt+1; + Clist(ii,jcnt+1)=jj; + end + end +end + +%Maximum number of MPA iterations before we declare a failure. +itenum = 50; +%Keep a history of number of iterations necessary for each decoding. +ite_hist = zeros(1,num_trials); +%Number of UnSatisfied Checks at each iteration. +num_USC = zeros(1, itenum); +%Store which checks unsatisfied. +USC_vec = zeros(1, n-k); +c_hat_hist_save = cell(1,itenum); +USC_vec_save = cell(1,itenum); +c_hat_hist = zeros(itenum, n); + +%Store Trapping Set, noise, and codeword information +num_TS = 0; +TS = zeros(100, n); +noise_save = zeros(100, n); +cw = zeros(10,n); + +%small_num defines the magnitude of largest messages from C->V +small_num = eps; + +%Number of miscorrections (MPA converges to incorrect valid codeword) +num_miscor = 0; + +% Simple H matrix consistency check. +num_edges_v = sum(Clist(:,1)); +num_edges_c = sum(Vlist(:,1)); +if (num_edges_c ~= num_edges_v) + text = ['error - row/col mismatch in H'] +end +num_edges = num_edges_c; +clear num_edges_c num_edges_v + +Lr_ind = zeros(1, num_edges); +Lq_ind = zeros(1, num_edges); + +%%%%% +% Set up the edge connections between Lr (C->V messages) and Lq (V->C). +% Remember, Matlab matrices are indexed by columns, so if Lr is of +% dimension (n-k)X max(d_c)+1, then Lr(n-k+1) is the same as Lr(1,2). +tot_count = 0; +for kk = 2:length(Clist(1,:)) + for jj=1:n + if(Clist(jj,kk) ~= 0) + tot_count = tot_count + 1; + Lq_ind(tot_count) = (kk-2)*n + jj; + end + end +end + +row_index = zeros(n-k, 1); +tot_count = 0; +for kk=2:length(Clist(1,:)) + for jj=1:n + if(Clist(jj,kk) ~= 0) + tot_count = tot_count + 1; + row=Clist(jj, kk); + row_index(row) = row_index(row) + 1; + Lr_ind(tot_count) = row + (row_index(row)-1)*(n-k); + end + end +end + +temp1 = zeros(size(Lr)); +temp2 = zeros(size(Lq)); +tic +%%%%% +% Main SNR loop +for ebnoloop = 1:length(ebnodb) + ebno = 10^((ebnodb(ebnoloop))/10) + esno = ebno*(rate); + sigma_2 = 1/(2*esno); + sigma = sqrt(sigma_2); + Lc_coef = 4*esno; + num_frame_errors = 0; + num_bit_errors = 0; + t = 0; + %Simulate till we collect 20 errors or num_trials reached w/o getting + % 20 errors. This can be adjusted to affect quality of estimate. +while ((t < num_trials) && (num_frame_errors < 20)) + t = t + 1; +%%% If encoding necessary +% u = (1-sign(rand(1,k)-0.5))/2; +% a = mod(u*G, 2); +% a_mod = -1*(a*2 - 1); + + noise = randn(1,n); + %for easy encoding, send all 0's message. + a_mod = ones(1,n); + + y = a_mod + sigma*noise; %regular MC + +% MPA +% Initialization +Lc=Lc_coef*y; +num_total = 0; +for ii = 1:num_v_deg + Lq(num_total+1:num_total+deg_v_prof(ii, 1), 1:deg_v_prof(ii, 2)) = Lc(num_total+1:num_total + deg_v_prof(ii, 1))'*ones(1,deg_v_prof(ii, 2)); + num_total = num_total + deg_v_prof(ii, 1); +end + + +stopsig=0; % If a valid codeword has been found, set this to `1' +itestep=0; % Iteration number + + +while (stopsig==0) && (itestep