Probabilistic Error Correction

%% Lab 9 — XOR 패리티로 사라진 가닥 되살리기 (툴박스 불필요)
A = double('SONG');  B = double('KIM!');
C = bitxor(A, B);                 % 패리티 가닥 C = A xor B
fprintf('C = %s (hex)\n', strjoin(string(dec2hex(C,2)), ' '));
B2 = bitxor(A, C);                % B 가 사라짐 → A xor C
fprintf('복원된 B = %s\n', char(B2));

%% 가닥 4개 + 패리티 1개 (RAID 5 와 같은 원리)
D = randi([0 255], 4, 8);         % 가닥 4개 x 8 byte
P = bitxor(bitxor(D(1,:), D(2,:)), bitxor(D(3,:), D(4,:)));
lost = 3;                         % 3번 가닥이 안 읽혔다
R = P;
for i = setdiff(1:4, lost)
    R = bitxor(R, D(i,:));        % 남은 것을 모두 XOR
end
fprintf('%d번 가닥 복원: %d\n', lost, isequal(R, D(lost,:)));

%% Lab 10 — LT(fountain) 부호: droplet 몇 개면 되나? (툴박스 불필요)
% 값 대신 "그 조각을 알아냈는가"만 추적해 필요한 droplet 수를 센다
for K = [20 100]                  % 원본 조각 수
    c = 0.1;  delta = 0.5;        % robust soliton 모수
    d = 2:K;
    rho = [1/K, 1./(d.*(d-1))];   % ideal soliton
    R = c*log(K/delta)*sqrt(K);  s = floor(K/R);
    tau = zeros(1, K);
    tau(1:s-1) = R./((1:s-1)*K);
    tau(s) = R*log(R/delta)/K;
    cdf = cumsum((rho+tau) / sum(rho+tau));
    need = zeros(1, 30);
    for trial = 1:30
        known = false(1, K);  eqs = {};
        while ~all(known)
            deg = find(rand <= cdf, 1);
            eqs{end+1} = randperm(K, deg);   % droplet 하나
            changed = true;                  % peeling 복호
            while changed
                changed = false;
                for e = 1:numel(eqs)
                    u = eqs{e}(~known(eqs{e}));
                    if numel(u) == 1
                        known(u) = true;  changed = true;
                    end
                end
            end
        end
        need(trial) = numel(eqs);
    end
    fprintf('K=%3d: 평균 %.1f개 필요 (여분 %.0f%%)\n', ...
            K, mean(need), 100*(mean(need)/K - 1));
end

%% Lab 11 — LDPC bit-flipping (툴박스 불필요)
H = [1 1 0 1 0 0;                 % 검사식 3개, bit 6개
     0 1 1 0 1 0;
     1 0 1 0 0 1];
c = [0 0 1 0 1 1];                % codeword (H*c' = 0)
r = c;  r(3) = ~r(3);             % 3번 bit 오류
for it = 1:5
    s = mod(H*r', 2);             % 깨진 검사식 = 1
    if ~any(s), break, end
    cnt = s' * H;                 % bit 마다 깨진 검사식 수
    [~, b] = max(cnt);
    fprintf('반복 %d: 깨진 식 %s, bit %d 뒤집기\n', ...
            it, num2str(s','%d'), b);
    r(b) = ~r(b);
end
fprintf('복원 %s, 정답과 같음 %d\n', num2str(r,'%d'), isequal(r,c));

%% (Communications Toolbox) 긴 LDPC + soft 정보(LLR)
H2  = dvbsLDPCPCM('1/2');        % 32400 x 64800 희소 행렬
enc = ldpcEncoderConfig(H2);
dec = ldpcDecoderConfig(H2);
msg = randi([0 1], enc.NumInformationBits, 1);
cw  = ldpcEncode(msg, enc);
sg  = 0.8;                        % 잡음 크기
y   = (1 - 2*cw) + sg*randn(size(cw));   % 0 → +1, 1 → -1
llr = 2*y / sg^2;                 % 양수면 0 쪽
nerr = sum((llr < 0) ~= cw);     % hard 판정 오류
fprintf('복호 전 bit 오류 %d / %d\n', nerr, numel(cw));
out = ldpcDecode(llr, dec, 50);
fprintf('LDPC 후 bit 오류 %d\n', sum(out ~= msg));

%% Lab 12 — convolutional 부호 (7,5) 와 Viterbi
u = [1 0 1 1 0 0];                % message 1011 + 꼬리 00
s = [0 0];                        % 상태 = 직전 두 bit
c = [];
for t = 1:numel(u)
    o1 = mod(u(t) + s(1) + s(2), 2);   % 생성기 7 = 111
    o2 = mod(u(t) + s(2), 2);          % 생성기 5 = 101
    c  = [c o1 o2];
    s  = [u(t) s(1)];             % 기억을 한 칸 민다
end
fprintf('손 부호기: %s\n', num2str(c, '%d'));

%% (Communications Toolbox) 같은 부호와 Viterbi 복호
tr = poly2trellis(3, [7 5]);      % 기억 3칸, 생성기 7, 5
c2 = convenc(u, tr);
fprintf('convenc 와 같은가: %d\n', isequal(c, c2));
r  = c2;  r([2 9]) = ~r([2 9]);   % bit 2개 오류
d  = vitdec(r, tr, 5, 'term', 'hard');
fprintf('받은 %s\n', num2str(r, '%d'));
fprintf('복호 %s\n', num2str(d(1:4), '%d'));

%% Lab 13 — 자리가 밀리는 오류와 VT 부호 (툴박스 불필요)
a = 'ACGTACGTAC';
b = a;  b(4) = [];                % 4번 염기 삭제
b(end+1) = 'A';                   % 길이만 맞춰 비교
fprintf('Hamming %d칸 vs edit distance %d\n', ...
        sum(a ~= b), editdist(a, b));

%% VT0(6): sum(i * x_i) mod 7 == 0 인 단어들
n = 6;
X = dec2bin(0:2^n-1) - '0';
C = X(mod(X*(1:n)', n+1) == 0, :);
fprintf('VT0(%d) codeword %d개\n', n, size(C, 1));
ok = 0;  tot = 0;
for i = 1:size(C, 1)
    for p = 1:n                   % 모든 자리를 하나씩 삭제
        y = C(i,:);  y(p) = [];
        ok = ok + isequal(vtdecode(y, n), C(i,:));
        tot = tot + 1;
    end
end
fprintf('단일 삭제 복원 %d / %d\n', ok, tot);

function d = editdist(a, b)       % Levenshtein 거리 (DP 표)
D = zeros(numel(a)+1, numel(b)+1);
D(:,1) = 0:numel(a);  D(1,:) = 0:numel(b);
for i = 2:numel(a)+1
    for j = 2:numel(b)+1
        D(i,j) = min([D(i-1,j)+1, D(i,j-1)+1, ...
                      D(i-1,j-1) + (a(i-1) ~= b(j-1))]);
    end
end
d = D(end, end);
end

function x = vtdecode(y, n)       % VT0(n) 단일 삭제 복호
w = sum(y);                       % 남은 1의 수
D = mod(-sum((1:n-1) .* y), n+1); % 모자란 checksum
if D <= w                         % 지워진 것은 0
    r = [cumsum(y, 'reverse') 0]; % 그 자리 오른쪽 1의 수
    p = find(r == D, 1);  ins = 0;
else                              % 지워진 것은 1
    l = [0 cumsum(1 - y)];        % 그 자리 왼쪽 0의 수
    p = find(l == D - w - 1, 1);  ins = 1;
end
x = [y(1:p-1) ins y(p:end)];
end

%% Lab 14 — homopolymer 제약의 용량 (툴박스 불필요)
for r = 1:5                       % 같은 염기 최대 r 개 연속
    T = zeros(r);                 % 상태 = 지금 run 길이
    T(:, 1) = 3;                  % 다른 염기 3개 → run 1
    for k = 1:r-1
        T(k, k+1) = 1;            % 같은 염기 → run + 1
    end
    lam = max(abs(eig(T)));       % 가짓수가 늘어나는 비율
    fprintf('run <= %d : %.3f bit/nt\n', r, log2(lam));
end

%% 무작위 150 nt 가 제약을 통과할 확률
N = 1e4;  ok = 0;
for i = 1:N
    s = randi(4, 1, 150);         % 1 A, 2 T, 3 G, 4 C
    runs = diff([0 find(diff(s) ~= 0) 150]);
    gc = mean(s >= 3);
    ok = ok + (max(runs) <= 3 && gc >= 0.4 && gc <= 0.6);
end
fprintf('run<=3, GC 40-60%% 통과율 %.1f%%\n', 100*ok/N);