디지털통신 Explorer

2012년 6월 5일 화요일

매트랩을 사용한 부가백색잡음환경에서의 QAM 시스템 BER 계산

개인적으로 외적 환경 변화가 많았던 기간을 지냈기도 했고, 잠시 게임 폐인이 되기도 했다.

한동안 이동통신공부에 탄력받았었는데, 지금은 갈피를 잘 못 잡겠고 정체기인것 같다.

아~ 정말 공부할 양도 많고, 어렵기도 하다~ 또 어떤방향으로 해야하나~

그래도 강박증을 버리되, 호기심으로 조금씩 정리해 보기로 한다!

오늘은 주제는 'Compute BER for a QAM System with AWGN Using MATLAB' 이다

매트랩사이트의 예제문서의 내용을 요약해고, 나름대로 부연설명을 추가해서

전반적인 통신의 절차를 확인해 볼 생각이다.

먼저,,, 원문 링크주소는 http://www.mathworks.co.kr/help/toolbox/comm/gs/a1067963807.html

참고로, 매트랩버전의 차이로 일부 명령어가 다르기도 한데, 대략의 의미는 통할것으로 본다.

문제1. 기저대역의 변조, 채널, 복조로 구성된 통신시스템을 이용한 이진데이터 stream을 처리한다. 시스템의 bit error rate (BER)를 구하고, 또한 전송 및 수신한 신호를 scatter plot에 표시하라!

문제2 송신기와 수신기에 펄스파형성형 및 정합필터를 수행하는 한 쌍의 square root raised cosine filters 를 사용하는 그레이코드된 변조로 수정하라!

문제3 주어진 구속장(constraint lengths)과 생성다항식(constraint lengths and generator polynomials)가진 convolutional coding 과 decoding을 포함할 수 있게 수정하라!

%% Setup : Define parameters

% Modified Plot, With Gray Coding
M = 16; % Number of points in constellation
k = log2(M);           % Number of bits per symbol
n = 5e5;                % Number of bits to process
nsamp = 4;            % Oversampling rate
hMod = modem.qammod('M',M,'SymbolOrder','Gray');   % Modulator object
mapping = hMod.SymbolMapping;                            % Symbol mapping vector
pt = hMod.Constellation;                                         % Vector of all points in constellation
scatterplot(pt);                                                       % Plot the constellation.

% Include text annotations that number the points.
text(real(pt)+0.1,imag(pt),dec2bin(mapping));
axis([-4 4 -4 4]); % Change axis so all labels fit in plot.

%% Signal Source : Create a binary data stream as a column vector.

x = randint(n,1); % Random binary data stream

%% Encoder : Define a convolutional coding trellis and use it to encode the binary data.

t = poly2trellis([5 4],[23 35 0; 0 5 13]); % Trellis

code = convenc(x,t); % Encode.

coderate = 2/3;

%% Bit-to-Symbol Mapping : Convert the bits in x into k-bit symbols.

% B. Do ordinary binary-to-decimal mapping.

xsym = bi2de(reshape(code,k,length(code)/k).','left-msb');

%% Modulation

y = modulate(modem.qammod(M),xsym); % Modulate using 16-QAM.

%% Filter Definition : Define filter-related parameters.

filtorder = 40; % Filter order

delay = filtorder/(nsamp*2); % Group delay (# of input samples)

rolloff = 0.25; % Rolloff factor of filter

% Create a square root raised cosine filter.

rrcfilter = rcosine(1,nsamp,'fir/sqrt',rolloff,delay);

% Plot impulse response.

figure; impz(rrcfilter,1);

%% Transmitted Signal : Upsample and apply square root raised cosine filter.

ytx = rcosflt(y,1,nsamp,'filter',rrcfilter);

% Create eye diagram for part of filtered signal.

eyediagram(ytx(1:2000),nsamp*2);

%% Channel : Send signal over an AWGN channel.

EbNo = 10; % In dB

snr = EbNo + 10*log10(k*coderate)-10*log10(nsamp);

ynoisy = awgn(ytx,snr,'measured');

%% Received Signal : Filter received signal using square root raised cosine filter.

yrx = rcosflt(ynoisy,1,nsamp,'Fs/filter',rrcfilter);

yrx = downsample(yrx,nsamp); % Downsample.

yrx = yrx(2*delay+1:end-2*delay); % Account for delay.

%% Scatter Plot : Create scatter plot of received signal before and after filtering.

h = scatterplot(sqrt(nsamp)*ynoisy(1:nsamp*5e3),nsamp,0,'g.');

hold on;

scatterplot(yrx(1:5e3),1,0,'kx',h);

title('Received Signal, Before and After Filtering');

legend('Before Filtering','After Filtering');

axis([-5 5 -5 5]); % Set axis ranges.

%% Demodulation : Demodulate signal using 16-QAM.

zsym = demodulate(modem.qamdemod(M),yrx);

%% Symbol-to-Bit Mapping : Undo the bit-to-symbol mapping performed earlier.

z = de2bi(zsym,'left-msb'); % Convert integers to bits.

% Convert z from a matrix to a vector.

z = reshape(z.',numel(z),1);

%% Decoder : Decode the convolutional code.

tb = 16; % Traceback length for decoding

z = vitdec(z,t,tb,'cont','hard'); % Decode.

%% BER Computation : Compare x and z to obtain the number of errors and the bit error rate. Take the decoding delay into account.

decdelay = 2*tb; % Decoder delay, in bits

[number_of_errors,bit_error_rate] = biterr(x(1:end-decdelay),z(decdelay+1:end))

awgn 의 함수의 검증은 아래와 같다.

SNR_input=30;

         % To specify the power of X to be 0 dBW and add noise to produce
         % an SNR of 10dB, use:
         N=2^10;
         X = sqrt(2)*sin(0:2*pi/N:2*pi-2*pi/N);
         Y = awgn(X,SNR_input,0);

X_pow=X*X'/N                 % =var(X)
Y_pow=Y*Y'/N                 % =var(Y)
plot(1:N, X, 1:N, Y)
noise_pow=sum((X-Y).^2)/N    % =var(X-Y)
SNR_dB=10*log10(X_pow/noise_pow)

먼저, 컨벌루션 코딩를 수행하는 하드웨어와 격자구조를 비교해 보자.

2개 입력에 3개 출력이 나오는 구조이다.

A Rate-2/3 Feedforward Encoder

(trel = poly2trellis([5 4],[23 35 0;0 5 13])

(The ratio of bit energy to noise power spectral density ) k * Eb/N0

= Es/N0 (the ratio of symbol energy to noise power spectral density)

(the signal-to-noise ratio) SNR = (Eb * k * coderate / nsamp) / N0

또한 각 단계별 정보의 변형(예시로 8bit )을 추적해 보자!

위 정보 가공의 과정을 요약하자면, 랜덤이진데이터x에서 채널코딩후 16-QAM변조 시키고

noise를 s/n비로 오염 후 업샘플하여 필터처리하여 송신한다.

수신단에서는 송신신호를 동일한 정합필터를 통하여 심볼을 검출하고 디코딩해서

원 랜덤이진데이터x와 비교하여 BER을 계산하면 끝~

덧붙여 Phase Noise(위상잡음)에 대해서도 알아보자!

'Compute BER for a QAM System with AWGN and Phase Noise Using Simulink' 이다.

원문 링크주소는,,, http://www.mathworks.co.kr/help/toolbox/comm/gs/fp8022.html

Carrier Frequency : 100MHz

Offset Frequency : 1KHz

Carrier/Offset Power : 60dBc 이면,

위상잡음은 -60dBc/KHz = -60-10log(1000) = -90dBc/Hz

중심주파수와 오프셋주파수에서의 전력차이가 클수록(Phase noise ↓)

BER이 우수하다는 것을 알 수 있다.

위상잡음의 영향은 Contellation 이 둥근패턴이 되어 검출의 정확성을 떨어뜨린다.

좀 더 자세한 실험은 원문링크주소를 참고하기 바란다.

이상,,, matlab의 '시작하기' 였다... 헐,,,

그래도 시작이 반이다.^^

이번 주제에는 RF 주파수천이와 동기화 방법들이 생략되었다.

동기화방법은 짧게 정리하자면,,,

Carrier Recovery - 수신측에서 복조에 사용되는 반송파(local carrier)의 주파수는 물론 그 위상도 수신신호의 carrier 위상에 맞춰져야 한다는 것

Symbol Synchronization(Timeing Recovery) - 수신측의 복조기출력을 한 심볼구간(Ts)마다 샘플링해서 detector로 보내는데, 이 샘플링시점도 신호가 송신기로부터 수신기까지의 전송시간 만큼 지연시켜 심볼구간에 맞도록 동기화

Frame Synchronization - 사전에 약속된 일정 길이의 접두부호(Prefix)를 매 프레임이 시작할 때마다 전송함으로써, 수신측에서 어느 비트가 MSB(Most Significant Bit)이고, 어느 비트가 LSB(Least Significant Bit)인지 구별

다음에는 이런것들과 측정방식에 대해서도 정리를 해야겠다.

일단,,, 기말고사 준비나 해야겠다...

2012년 2월 23일 목요일

모스부호(Morse code)와 Tip

간만에 블로그를 올린다.

이번에는 직접 제작한 매트랩 프로그램과 알면 유용한 매트랩 명령어에 대하여 알아보자!

그러기에 앞서, 항상 느끼는 것이지만,,,

매트랩은 공대생들에게 쉽고 유용한 필수도구라고 생각이든다.

짧은코딩으로, 쉽게 목적을 달성할 수 있는지 미리 살필 수 있고, 원리를 파악하기도 쉽다.

하드웨어랑 연결하고 조작도 쉬운것 같은데,,, 아직은 잘 모르겠다.

분명한것은 공부해야할 것들이 무궁무진하게 많고,

쉽게 코딩해서 강력한 기능을 맘대로 사용할 수 있다는 점이다.

수학과 공학 지식의 축적으로 이루어진 거대한 탑...^^

너무 높고, 방대해서 한 사람의 노력으로는 다 이해하기 어려운,,,

아무튼, 멋진 공학의 기초 도구이가 될것이다!ㅎ

그럼 이제, K5 자동차 광고에 이용해서 고급스러움과 신비스러움을 더 해줬던

모스부호의 세계로 빠져봅시다~

아래 광고의 '뚜~ 뚜 뚜~ 뚜 뚜 뚜 뚜 뚜' 의 세련된 음은

모스부호로써 'K5' 를 연속해서 가르키고 있다.

아래코드는 자작한 모스부호생성, 모스부호분석하여 변환하는 프로그램이다.</><>

아직 wave 파일 분석은 일반적이지 못해 좀더 수정해 볼 생각이다.

table={'a' '.-';
'b' '-...';

'c' '-.-.';

'd' '-..';

'e' '.';

'f' '..-.';

'g' '--.';

'h' '....';

'i' '..';

'j' '.---';

'k' '-.-';

'l' '.-..';

'm' '--';

'n' '-.';

'o' '---';

'p' '.--.';

'q' '--.-';

'r' '.-.';

's' '...';

't' '-';

'u' '..-';

'v' '...-';

'w' '.--';

'x' '-..-';

'y' '-.--';

'z' '--..';

'1' '.----';

'2' '..---';

'3' '...--';

'4' '....-';

'5' '.....';

'6' '-....';

'7' '--...';

'8' '---..';

'9' '----.';

'0' '-----'};

code=[];

c=input('문자를 입력해주세요!\n', 's'); c=lower(c);

for i=1:length(c)

for j=1:length(table)

if strcmp(c(i), table{j,1}), code=[code table{j,2}]; end

end

code=[code ' '];

end

fs=8000;

l=0.1;

dot_t=0:1/fs:l; %

dash_t=0:1/fs:3*l; % -

dot=sin(2*pi*1812.5*dot_t);

dash=sin(2*pi*1812.5*dash_t);

space=zeros(1,length(dot));

code

morse=[];
for i=1:length(code)

if strcmp(code(i),'.')

morse=[morse dot];

elseif strcmp(code(i),'-')

morse=[morse dash];

else

morse=[morse space space];

end

morse=[morse space];

end

sound(morse,fs)

% wavwrite(morse,fs,16, c)

% morese=wavread(c)

ev=(abs(hilbert(morse))>0.5);
% s(t)=Re[s_hat(t)·exp(jwt)]=si(t)cos(wt)-sq(t)sin(wt) , s_hat(t)=si(t)+j·sq(t)

evs=downsample(ev,fs*l/8);

re=[]; k=1;

while k~=length(evs)

if length(evs)>=k+23 & sum(evs(k:k+23))==24, re=[re '-']; k=k+24;

elseif length(evs)>=k+23 & sum(evs(k:k+23))==0, re=[re ' ']; k=k+24;

elseif length(evs)>=k+7 & sum(evs(k:k+7))==8, re=[re '.']; k=k+8;

end

k=k+1;

end

id=[1 strfind(re, ' ')+1];

rec=[]; i=1; j=1;

for i=1:length(id)-1

for j=1:length(table)

if strcmp(re(id(i):id(i+1)-2), table{j,2}), rec=[rec table{j,1}]; end

end

upper(rec)

별건아니지만, hello world의 순수함으로 만들고, 뿌듯함을 느낀다.ㅎㅎ

다음은 잘 몰랐지만, 알고나면 굉장히 유용한 matlab의 주요한 기능을 살펴보자!

%-----------------------------------------------

% Create a new AVI file
%-----------------------------------------------
t = linspace(0,2.5*pi,40);

fact = 10*sin(t);

fig=figure;

aviobj = avifile('example.avi')

[x,y,z] = peaks;

for k=1:length(fact)

h = surf(x,y,fact(k)*z);

axis([-3 3 -3 3 -80 80])

axis off

caxis([-90 90])

F = getframe(fig);

aviobj = addframe(aviobj,F);

end

close(fig)

aviobj = close(aviobj);

%-----------------------------------------------

% Display wait bar

%-----------------------------------------------

h = waitbar(0,'Please wait...');

for i=1:1000,

% computation here %

waitbar(i/1000,h)

end

%-----------------------------------------------

% Display legend

%-----------------------------------------------

x = 0:.2:12;

plot(x,bessel(1,x),x,bessel(2,x),x,bessel(3,x));

legend('First','Second','Third');

legend('First','Second','Third','Location','NorthEastOutside')

%-----------------------------------------------

% meshgrid, surf, contour

%-----------------------------------------------

[X,Y] = meshgrid(-1:.1:1);

Z = 1-(X.^2 + Y.^2);

surf(X,Y,Z)

surfc(X,Y,Z)

[C, h]=contour(X, Y, Z);

clabel(C, h)

save('f.mat', 'X', 'Y', 'Z');

load('f.mat')
%-----------------------------------------------

%-----------------------------------------------

% Record and play back 2 seconds of single audio sampled at 8 kHz

%-----------------------------------------------

Fs = 8000;

y = wavrecord(2*Fs, Fs, 'Single');

wavplay(y, Fs);
%-----------------------------------------------

%-----------------------------------------------

% Calendar

%-----------------------------------------------

calendar(date)
%-----------------------------------------------

%-----------------------------------------------

% Capture

%-----------------------------------------------

vid = videoinput('winvideo', 1, 'YUY2_160x120');

src = getselectedsource(vid);

vid.FramesPerTrigger = 1;

start(vid);

vid.ReturnedColorspace = 'rgb';

nicelcw = getdata(vid);

imshow(nicelcw);
%-----------------------------------------------

%-----------------------------------------------

% Font color

%-----------------------------------------------

a='test';

text(1,1,a,'color','r')
%-----------------------------------------------

이상,,, 블로그에 업로드하기까지 오랜기간이 걸렸다. 그동안 항상 맘에 있었지만,,,

게을렀다고 할까?!ㅎ

한편으로는 허접한 블로그를 업로드하기위해, 나는 지난날을 그렇게 끙끙댔나보다~

올해,,, 만학도로 대학원을 입학했다.

몇일해본 대학원은 해주는 것을 받아먹는게 아니라,

노력하는 만큼 얻어가는 곳이라는 느낌이 확 든다.

공부해서,,, 졸업 쯔음엔 떳떳하고 멋있게 내 논문 발표할 수 있으면 좋겠다^^

그리고 좋은 블로그를 게시하여,

디지털통신 전공생의 가려운 부분을 박박 긁어주고 싶고, 지식도 나누고 싶다.

나 자신을 위하여 파이팅~!!ㅎ

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

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

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)을 확인 할 수도 있다.