2011년 11월 9일 수요일

채널코딩

채널코딩은 잡음에 오염된 채널을 통해 신뢰할 만한 통신을 성취하기 위하여, 소스 데이터를 변환하는 것을 말한다. 신뢰도를 높이는 대신 전송률과 대역폭, 복잡도는 절충해야 한다.

전송 오류제어하는 방법은,,,
수신측에서 독자적으로 오류를 일으킨 비트를 정정(FEC, 순방향 오류 정정),
수신단에서 오류 검출을 하고, 오류 검출될 때마다 역방향 링크를 통해 송신측에 보고하여 데이터를 재전송하도록 하는 방식(ARQ, 자동 반복 요구)

FEC 채널코딩을 분류하자면,,,
wave coding , block coding, convolution coding 이 있다.
블록 코딩은 현재의 k비트 메시지어는 과거의 k비트 메시지어와 관계 없이 독립적으로 출력을 만드어 내는데 반하여,
컨볼루션 코딩은 현재의 부호어를 만드는데 있어서 현재의 k비트 메시지어만 필요한 것이 아니라 과거의 메시지어도 필요하다. 컨볼루션 코딩은 처리 지연이 작게 되어 실시간 통신이 중요시 되는 경우에 많이 사용된다.

인터리빙(interleaving)이 지향하는 효과는 연속된 오류를 분산시켜서 한 부호어 내에는 오류정정 능력 이하의 비트 오류가 발생되도록 하는 것이다.

채널용량, 정보용량
엔트로피 H(x) 채널입력이 관찰되기 전에 채널입력의 불확실성을 나타내고, 조건부 엔트로피 H(x|y)는 채널출력이 관찰된 후에 채널 불확실성을 표현하므로 H(x)-H(x|y)의 차이는 채널출력을 관찰함으로써 분석되는 채널입력에 대한 불확실성을 표현한다.
그리고 이것을 상호정보량 I(x;y)이라고 한다.


I(x;y)=H(x)-H(x|y)
결합엔트로피 I(x;y)=H(x)+H(y)-H(x,y) , H(x,y)=∑∑p(x,y)*log2(1/(p(x,y))

이산 무기억 채널의 채널용량을 채널의 단일 사용에 있어 최대 평균 상호정보로 정의한다.
C=max[I(x;y)]

채널용량 C는 채널사용당 비트, 혹은 전송당 비트로서 측정된다.
 




2진대칭채널(BSC)

채널은 두 개의 입력심벌과 두 개의 출력심벌을 갖는다. 채널은 0이 보내졌을 때 1로 받을 확률은 1이 보내졌을 때 0으로 받을 확률과 같기 때문에 대칭적이다. 오류의 조건부 확률은 p로 정의되어진다.

H(x)는 채널입력 확률이 1/2일 때, 최대(1)로 된다.

C=1+p*log2(p)+(1-p)*log2(1-p) =1-H(p)

① 채널은 잡음이 없을 때, p=0이라 하면 채널용량 C는 채널사용당 한 비트의 최대값이 되는데 그것은 각 채널입력에서의 정확한 정보가 된다. 위의 p값에서 엔트로피 함수 H(p)는 최소값 0에 도달한다.
② 잡음에 의한 조건부 오류확률 p=1/2일 때 채널용량 C는 최소치 0에 도달하고 반면에 엔트로피 함수 H(p)는 최대값 1에 도달하는데 그러한 경우의 채널은 쓸모 없다.

p=linspace(0,1,100); C=1+p.*log2(p)+(1-p).*log2(1-p); plot(p,C)
 





Shannon's 2st theorem : 2진 대칭채널에서 채널용량 C와 같거나 이하인 어떤 부호율 r(code rate)에 대해서도 원하는 평균 오류확률보다 작은 오류확률을 가지는 부호가 존재하한다.

(예시) 반복부호에서 코드 rate에 대한 평균 오류확률
r=1일때, Pe는 1/100 , r=1/3일때, Pe=∑i=2~3 (3C2)*p^i*(1-p)^(3-i) ∴ Pe=3/10000



가우시안 확률분포함수(pdf) f_X (x)=1/(σ√2π) e^([-(x-m)^2/2σ^2 ])

평균 0, 분산 1로 정규화된 가우시안 누적분포함수(=Q함수)

Q(x)=∫_x~∞[1/√2π e^([-u^2/2]) du]

erfc(x)=2Q(√2x) ※ Q(0)=0.5 , Q(1)=0.1587 , erfc(x)=1-erf(x)

비트오류확률 Pb=Q(sqrt(Eb((1-p)/N0)) , p는 상관계수 ∴ BPSK의 Pb=Q(sqrt(SNR))

Xk가 전송된 신호의 표본이라고 한다. 채널의 출력은 평균이 영이고 전력스펙트럼밀도가 N0/2의 부가백색 가우시안 잡음(AWGN)에 의해 방해를 받는다. 잡음도 B Hz로 대역제한된다. 연속랜덤변수 Yk=Xk+Nk 가 수신된 신호의 표본이라고 하자

잡음표본 Nk는 평균이 영이고, 분산이 σ^2=N0*B , E[Xk^2]=P

C=max[I(Xk;Yk):E[Xk^2]=P] =H(Yk)-H(Yk|Xk)

C=H(Yk)-H(Nk) =1/2*log2(2πe(P+σ^2)-1/2*log2(2πeσ^2) =1/2*log2(1+P/(σ^2))

단위시간당 정보용량 C=B*log2(1+P/(N0*B))
 


Shannon's 3st theorem : 시스템이 전력제한되고, 대역제한된 가우시안 채널에 대해 오류 없는 전송의 속도에 대한 기본적인 제한을 정의
C=B*log2(1+P/(N0*B)) bit/sec

(예시) 대역폭이 19.2[kHz]인 RF채널의 용량이 9.6[kbps]가 되도록 하기 위해 필요한 신호의 SNR은? 9600=19200*log2(1+SNR)[bps]  ∴ SNR=-3.83[dB]


Linear Block Code
Hamming 거리 d는 임의의 두 부호어 벡터의 같은 위치에 있는 원소들 중 서로 다른 원소의 개수로 정의된다.
오류검출능력 q=d_min-1 , 오류정정능력 t=(d_min-1)/2

블록코딩은 k 비트의 message vector를 n 비트의 code vector로 매핑시키는 것을 말하는데, 길이 n의 codeword 2^k개로 구성된 code 를 (n, k) block code라고 한다.
c=mG , G_k*n=[I_k*k|P_k*(n-k)] , H=[P'_(n-k)*k|I_(n-k)*(n-k)] , G*H'=P+P=0
syndrome vector s=r*H'=(c+e)H'=e*H'
c=r-e
디코더는 s에 해당하는 에러유형 e를 찾아서 그것을 수신벡터 r로부터 뺌으로써 원래의 올바른 code vector c를 구한다.




% (7, 4) 블록코딩

SNRbdB=10; MaxIter=1e6;

n=3; N=2^n-1; K=2^n-1-n; % Codeword (Block) length and Message size

Rc=K/N; % Code rate

SNRb=10.^(SNRbdB/10);

SNRbc=SNRb*Rc; % SNRbc*n=SNR*k

sqrtSNRbc=sqrt(SNRbc);

pemb_uncoded=Q(sqrt(SNRb)); % Uncoded msg BER with BPSK

et=Q(sqrt(SNRbc)); % Crossover probability

L=2^K;

for i=1:L
M(i,:)=deci2bin1(i-1,K); % All message vectors

end

P=[1 1 0; 0 1 1; 1 1 1; 1 0 1];

G=[P eye(K)];, H=[P' eye(N-K)]; % [H,G]=Hammgen(n);

Hamming_code=rem(M*G,2);

Min_distance =min(sum(Hamming_code(2:L,:)')); % weight

No_of_correctable_error_bits=floor((Min_distance-1)/2); % 오류정정비트수

E= eye(N); % Error patterns 에러유형

S= rem(E*H',2); NS=size(S,1); % The syndrome matrix

nombe=0;

for iter=1:MaxIter
msg=randint(1,K); % Message vector
coded= rem(msg*G,2); % Coded vector
modulated= 2*coded-1; % BPSK-modulated vector
r= modulated +randn(1,N)/sqrtSNRbc; % Received vector with noise
r_sliced= r>0; % Sliced

r_c=r_sliced; % To be corrected only if it has a syndrome

s= rem(r_sliced*H',2); % Syndrome

for m=1:NS % Error correction depending on the syndrome

if s==S(m,:), r_c=rem(r_sliced+E(m,:),2); break; end

end
notbe1=sum(coded~=r_sliced); notbec1=sum(coded~=r_c);

% 2개이상 비트 오차가 발생했을 때 정정은 오히려 오차개수를 증가시킴

nombe=nombe+sum(msg~=r_c(N-K+1:N)); if nombe>100, break; end

end

pemb=nombe/(K*iter); % Message bit error probability

pemb_t=prob_err_msg_bit(et,N,No_of_correctable_error_bits);

fprintf('\n Uncoded Messsage BER=%8.6f',pemb_uncoded)

fprintf('\nMessage BER=%5.4f(Theoretical BER=%5.4f)\n',pemb,pemb_t)

G, Hamming_code, H', [E S]

 

G =
     1     1     0     1     0     0     0
     0     1     1     0     1     0     0
     1     1     1     0     0     1     0
     1     0     1     0     0     0     1

Hamming_code =
     0     0     0     0     0     0     0
     1     0     1     0     0     0     1
     1     1     1     0     0     1     0
     0     1     0     0     0     1     1
     0     1     1     0     1     0     0
     1     1     0     0     1     0     1
     1     0     0     0     1     1     0
     0     0     1     0     1     1     1
     1     1     0     1     0     0     0
     0     1     1     1     0     0     1
     0     0     1     1     0     1     0
     1     0     0     1     0     1     1
     1     0     1     1     1     0     0
     0     0     0     1     1     0     1
     0     1     0     1     1     1     0
     1     1     1     1     1     1     1

H' =
     1     1     0
     0     1     1
     1     1     1
     1     0     1
     1     0     0
     0     1     0
     0     0     1

[E | S] =
     1     0     0     0     0     0     0     1     1     0
     0     1     0     0     0     0     0     0     1     1
     0     0     1     0     0     0     0     1     1     1
     0     0     0     1     0     0     0     1     0     1
     0     0     0     0     1     0     0     1     0     0
     0     0     0     0     0     1     0     0     1     0
     0     0     0     0     0     0     1     0     0     1


 참고로 BER은 berawgn 과 비교해 보면 큰 향상이 있음.


% BER of a block code

n = 23; k = 12; % Lengths of codewords and messages
dmin = 7; % Minimum distance
EbNo = 1:10;
ber_block = bercoding(EbNo,'block','hard',n,k,dmin);
berfit(EbNo,ber_block) % Plot BER points and fitted curve.
ylabel('Bit Error Probability');
title('BER Upper Bound vs. Eb/No, with Best Curve Fit');




Convolutional Code와 Vitebi Decodeing
컨볼루션 코딩은 정보 비트열이 직렬로 채널 부호화기에 입력되고 전송량이 증가한 부어열이 직렬로 출력되어 전송되는 방식이다. (n, k , K)코드로 표현하는데, K는 n비트의 부호어를 생성하기 위하여 필요한 메시지어의 개수를 나타내고, 이를 구속장이라 한다.

(2, 1, 3) 컨볼루션 부호화기를 예로 살펴보자!
부호화기에 연결벡터는 g1=(111), g2=(101)





trel = poly2trellis(3, [7 5]) % constraint length(구속장) K=3 , 1 1 1 ->7 , 1 0 1 ->5

msg = [1 1 0 1 1 0 0]

[code state]=convenc(msg,trel), tblen = 1; % traceback length

out = vitdec(code, trel, tblen, 'cont', 'hard')


code =

1 1 0 1 0 1 0 0 0 1 0 1 1 1


state =
0

out =
0 1 1 0 1 1 0

비터비 복호 알고리즘은 최우복호를 위한 연산량을 줄일 수 있는 효율적인 방법을 제시한다. 비터비 알고리즘은 각 순간마다 각 상태로 들어오는 트렐리스 경로들과 수신 코드열의 거리를 비교하여 최우경로의 가능성이 없는 경로는 고려 대상에서 제거함으로써 연산량을 줄일 수 있다.

이상,,, 상당히 이해하기 난해한 부분이 많이 있지만, 가능한 감이라도 오도록 참고내용을 요약해서 그대로 적는다. 디지털 데이터전송을 잘하기 위해 비트열을 가지고 이렇게 줄였다 늘렸다 곱했다 나눴다 찌지고 뽂는 가공(프로세싱)을 하는 것을 보면, 정보에는 공장에서 찍어만드는 공산품보다도 더 눈물겨운 수고로움이 축척되어 있음을 알게된다.
마지막으로 비터비란 통신의 대가는 퀄컴의 창업주라고 한다. 역시 자신의 알고리즘을 CDMA 휴대폰에 전송에러를 획기적으로 줄이데 활용했으리라. 뿐만아니라 여러 공학적인 응용에서 비터비알고리즘이 많이 쓰이는 듯 하다. 책에서 보던 공학자가 세계의 몇백번째 부자로 사회적으로 성공하고, 아직도 실존하고 있다는 것이 나에게는 약간은 충격적이고, 부럽다.ㅎ

2011년 11월 7일 월요일

소스코딩

소스코딩은 어떤규칙에 의해서 정보를 압축시켜 전송량을 줄이고자하는 절차이다.

엔트로피는 심벌당 평균 정보량 이다. (불확실성의 측정량)

H(x)=∑p∙log_2 (1/p)

엔트로피는 사건의 확률이 같을때, 최대가 된다.

 (예시)  00, 01, 10, 11 네개의 심볼이 p=1/4 확률로 발생하면, 엔트로피는 2 bits

Shannon's 1st theorem : 정보를 코딩하는데 필요한 bit의 개수는 적어도 그 엔트로피 이상 필요



% 허프만 코딩

p=[0.2 0.15 0.13 0.12 0.1 0.09 0.08 0.07 0.06];    % 어떤 사건이 일어날 확률들

h={'1' '0'};

M=length(p);  N=M-1; p=p(:); % make p a column vector

if M>2

  pp(:,1)=p;

  for n=1:N

     [pp(1:M-n+1,n),o(1:M-n+1,n)]=sort(pp(1:M-n+1,n),1,'descend'); % in descending order

     if n==1, ord0=o; end  % Original descending order

     if M-n>1

       pp(1:M-n,n+1)=[pp(1:M-1-n,n); sum(pp(M-n:M-n+1,n))];

     end

  end

  for n=N:-1:2

     tmp=N-n+2; oi=o(1:tmp,n);

     for i=1:tmp, ht{oi(i)}=h{i}; end

     h=ht; h{tmp+1}=h{tmp};

     h{tmp}=[h{tmp} '1'];

     h{tmp+1}=[h{tmp+1} '0'];

  end

  for i=1:length(ord0), ht{ord0(i)}=h{i}; end

  h=ht;

end

L=0;

for n=1:M

  L=L+p(n)*length(h{n});% Average codeword length

end

H=-sum(p.*log2(p)); % Entropy


>> L,H,pp,o,h
L =    3.1000
H =     3.0731

pp =
0.2000 0.2000 0.2000 0.2200 0.2600 0.3200 0.4200 0.5800
0.1500 0.1500 0.1700 0.2000 0.2200 0.2600 0.3200 0.4200
0.1300 0.1300 0.1500 0.1700 0.2000 0.2200 0.2600 0
0.1200 0.1300 0.1300 0.1500 0.1700 0.2000 0 0
0.1000 0.1200 0.1300 0.1300 0.1500 0 0 0
0.0900 0.1000 0.1200 0.1300 0 0 0 0
0.0800 0.0900 0.1000 0 0 0 0 0
0.0700 0.0800 0 0 0 0 0 0
0.0600 0 0 0 0 0 0 0

o =
1 1 1 6 5 4 3 2
2 2 7 1 1 1 1 1
3 3 2 2 2 2 2 0
4 8 3 3 3 3 0 0
5 4 4 4 4 0 0 0
6 5 5 5 0 0 0 0
7 6 6 0 0 0 0 0
8 7 0 0 0 0 0 0
9 0 0 0 0 0 0 0

h =
Columns 1 through 8
'00' '110' '101' '011' '010' '1111' '1110' '1001'
Column 9
'1000'

부연설명을 하자면, 1 ~9까지의 심벌의 발생확률이 p라고 했을때,
100심벌을 전송한다고 가정해 보자!
메시지를 2진코드로 4bit씩 그대로 코딩하면, 400bit 가 필요한 반면,
허프만코딩은 20*2+(15+13+12+10)*3+(9+8+7+6)*4=310bit 로 정보를 압축할 수 있다.

이와달리, 심벌들의 발생과 확률분포가 필요없이 주어진 소스 시퀀스에 대한 codetable을 만들어 사용하는 Lempel-Ziv-Welch 코딩도 있다.
Lempel-Ziv 알고리즘은 허프만을 대신해 파일압축을 위한 기본으로 널리 이용되고 있다.

지금까지는 통계적 중복성 제거 를 통한 정보압축이었다.


시각적/청각적 특징 을 기반한 정보의 압축으로는
영상의 디지털화 형식(원래 RGB=8:8:8 → Ycbcr=8:2:2, 8:4:4), 음성의 마스킹 효과가 있다.


시간적 중복성 제거
영상에서 움직임 예측 및 움직임 보상 이론


공간적 중복성 제거
예로는, JPEG의 압축기술인 DCT를 들 수 있다.

DCT(Discrete Cosine Transform)은 직교 변환을 이용한 것으로, 공간상의 영상 데이터를 주파수상의 에너지 분포로 변환시키는 것이다.

F(u,v)=(1/4)∙C(u)C(v)∑_(i=0~7)∑_(j=0~7)f(i,j)cos((2i+1)uπ/16) cos((2j+1)vπ/16)
여기서 C(x)={(1/√2 ,x=0 @ 1 ,otherwise)}



Xr=[]; X=[];
for u=0:7
           for v=0:7
                     for i=0:7
                                for j=0:7
                    x(i+1,j+1)=cos((2*i+1)*u*pi/16).*cos((2*j+1)*v*pi/16);
                                end
                     end
                     Xr=[Xr x];
                subplot(8,8,u*8+v+1), imshow(x,[-1,1])
           end
           X=[X; Xr]; Xr=[];
end
 



64개의 DCT 기본함수


용량큰 JPEG를 화면에 뿌려줄때, 흐릿하다가 좀 뒤에 또렷해지는 것은...
위에 저주파의 신호부터 지그재그로 읽어들이는 도중에 먼저 화면에 뿌리고, 고주파의 신호까지 완벽히 읽은 뒤에 마지막으로 다시 화면에 뿌려주기 때문일 것이다.

부연설명을 하자면, 영상을 DCT변환하면 수직·수평 저주파 계수들의 값은 크고, 고주파로 갈 수록 값이 작아지는데, 이를 양자화하는데서 압축을 하면 정보량을 많이 줄일 수 있다.
 

RGB = imread('autumn.tif');
        I = rgb2gray(RGB);
        J = dct2(I);
        imshow(log(abs(J)),[]), colormap(jet), colorbar

        J(abs(J)<10) = 0;
        K = idct2(J);
        figure, imshow(I)

        figure, imshow(K,[0 255])
 
 


웨이블릿 변환

 정지영상압축 표준인 JPEG2000에서 매우 중요한 부분을 차지하고 있으며, MPEG-4, H.264등에도 사용되고 있다. 또한 잡음제거, 에지 검출과 같은 다양한 영상신호 처리에 응용된다.
웨이블릿 변환에 사용되는 기저 함수의 집합은 하나의 기본 웨이블릿 기저 함수에 대한 시간축 방향으로의 확대 및 축소 그리고 평행 이동을 통해 얻어진다. 기본 웨이블릿 기저 함수는 특별한 형태의 band-pass 필터로 생각할 수 있다. 웨이블릿 변환에서는 주파수 대역이라는 용어 대신에 스케일(scale)이라는 용어를 주로 사용한다.


Ψ_j,k(t) : 웨이블릿 크기 압축계수(scale, j), 시간축으로의 이동 전이계수(translattion, k)

DWT : f(t)=∑_(j⊂Z)C_0*Φ_0,j(t)+∑_(k≥0))∑_(j⊂Z)d_k,j*Ψ_k,j(t)



load wbarb;
wavelet = 'haar'; % Define wavelet of your choice
level = 2; % Define wavelet decomposition level
[C S] = wavedec2(X,level,wavelet); % Compute multilevel 2D wavelet decomposition
% decomposition vector C = [ A(N) | H(N) | V(N) | D(N) | ... | H(1) | V(1) | D(1) ]
% the corresponding bookkeeping matrix S


% Define colormap and set rescale value
colormap(map); rv = length(map);
% Plot wavelet decomposition using square mode

A = cell(1,level); H = A; V = A; D = A;
for k = 1:level
A{k} = appcoef2(C,S,wavelet,k); % Extract 2-D approximation coefficients
[H{k} V{k} D{k}] = detcoef2('a',C,S,k); % Extract 2-D detail coefficients
A{k} = wcodemat(A{k},rv); % Extended pseudocolor matrix scaling % approximation coefficients
H{k} = wcodemat(H{k},rv); % hori. detail coefficients
V{k} = wcodemat(V{k},rv); % vert. detail coefficients
D{k} = wcodemat(D{k},rv); % diag. detail coefficients
end
dec = cell(1,level);
dec{level} = [A{level} H{level} ; V{level} D{level}];

for k = level-1:-1:1
dec{k} = [imresize(dec{k+1},size(H{k})) H{k} ; V{k} D{k}];
end
figure(1);
image(dec{1});
title(['Decomposition at level ',num2str(level)]);

% Peak Signal To Noise Rate
X0 = waverec2(C,S, wavelet);
error=X-X0;
s_error=error^2;
ms_error=sum(sum(s_error),2)/(256*256);
Peak_SNR=10*log(255*255/ms_error)

% Analyze Frequency Characteristic of above wavelet filters
[D_L, D_H, R_L, R_H] = wfilters(wavelet);
den = [1];
figure(2);
freqz(D_L, den);
figure(3);
freqz(D_H, den);
figure(4);
subplot(221); stem(D_L);
title('Decomposition low-pass filter');
subplot(222); stem(D_H);
title('Decomposition high-pass filter');
subplot(223); stem(R_L);
title('Reconstruction low-pass filter');
subplot(224); stem(R_H);
title('Reconstruction high-pass filter');
xlabel(['The four filters for ' wavelet])
 

 




마찬가지로,,, 대부분의 정보는 왼쪽 위에 모서리에 집중되어 있다.

참고로, wfilters의 종류에 따른 웨이블릿 필터의 주파수응답(Fig2, 3)과 임펄스응답(Fig4)을 확인 할 수도 있다.