MATLAB PSK(Phase-shift keying) BER simulation

아래 포스팅에서 BPSK, QPSK BER(Bit error rate) simulation in AWGN channel 에 대해 설명 드렸었는데~

   

2011/03/27 - [programming language/MATLAB] - MATLAB QPSK BER simulation in AWGN channel

2011/03/19 - [programming language/MATLAB] - MATLAB BPSK BER simulation in AWGN channel

 

8 PSK(Phase-shift keying) 시뮬레이션에 대해 질문하신 분이 있어서 답변 드립니다.

 

먼저 제 경험을 말씀 드리면, 8 PSK 나 16 PSK 의 경우 책에서 공부한적은 있지만 실제 통신 시스템에서 사용하는 것을 본적은 없는 것 같습니다. 사용하지 않는 이유에 대해서는 아래 글을 참조해 보시기 바랍니다. 아래 글은 16 QAM 과 16 PSK 의 SER (Symbol Error Rate)성능 비교에 대한 내용인데 16 PSK 가 16 QAM 대비하여 확실히 SER 성능이 안 좋다는 것을 확인 할 수 있습니다.

 

http://www.dsplog.com/2008/03/29/comparing-16psk-vs-16qam-for-symbol-error-rate/

 

일단 QPSK, BPSK 의 경우 각 constellation 지점들이 90 도 또는 180 도 위상차이를 보이기 때문에 심볼에 대해 위상을 구할 필요가 없이 BPSK 의 경우 간단하게 x=0 지점을, QPSK 의 경우 x=0, y=0 지점을 기준으로 심볼 demapping 을 할 수 있었습니다.

 

하지만 8 PSK 의 constellation 은 아래 그림과 같이 8개 지점(45도 간격)으로 구성되어 있습니다. 이런 경우에는 당연히 45/2 = 22.5 도 간격으로 symbol demapping 을 해야 하므로 수신 심볼의 위상을 구해야 할 것입니다.

 

이미지 출처 : http://commons.wikimedia.org/wiki/File:8PSK_Gray_Coded.svg

 

PSK 시뮬레이션을 위해서는 PSK 의 Gray code 에 대해 공부를 하셔야 하고~ 아래 주소에 굉장히 좋은 설명이 있더군요. MATLAB/OCTAVE 코드도 있으니 실습해 보시면서 천천히 공부해 보시기 바랍니다.

 

http://www.dsplog.com/2008/05/12/gray-code-to-binary-conversion-for-psk-pam/

 

http://www.dsplog.com/2008/05/11/binary-to-gray-code-conversion-psk-pam/

 

다음으로 질문하신 8 PSK BER 시뮬레이션에 대해 소개하겠습니다.

 

아래 주소에 16 PSK BER 시뮬레이션에 대해 소개하는데~

 

http://www.dsplog.com/2008/05/18/bit-error-rate-for-16psk-modulation-using-gray-mapping/

 

코드는 아래 주소에서 다운로드 받을 수 있습니다.

 

http://www.dsplog.com/db-install/wp-content/uploads/2008/05/script_16psk_gray_mapping_bit_error_rate.m

 

위 코드를 실행해보면~ 16 PSK 뿐만 아니라 8 PSK, QPSK도 시뮬레이션 해 볼 수 있습니다.

 

위 주소에서 설명한 16 PSK 부터 시뮬레이션을 해보죠~ 아래 코드의 모든 권리는 http://www.dsplog.com 사이트의 Krishna Pillai 에게 있습니다. 재 배포를 하실 경우 원저자를 명시 하시기 바랍니다.

 

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

% All rights reserved by Krishna Pillai, http://www.dsplog.com

% The file may not be re-distributed without explicit authorization

% from Krishna Pillai.

% Checked for proper operation with Octave Version 3.0.0

% Author    : Krishna Pillai

% Email        : krishna@dsplog.com

% Version    : 1.0

% Date        : 17 May 2008

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

 

% Bit Error Rate for 16-PSK modulation using Gray modulation mapping

 

clear

N = 10^5; % number of symbols

M = 16; % constellation size

k = log2(M); % bits per symbol

 

 

thetaMpsk = [0:M-1]*2*pi/M; % reference phase values

 

Eb_N0_dB = [0:25]; % multiple Es/N0 values

Es_N0_dB = Eb_N0_dB + 10*log10(k);

 

% Mapping for binary <--> Gray code conversion

ref = [0:M-1];

map = bitxor(ref,floor(ref/2));

[tt ind] = sort(map);

 

ipPhaseHat = zeros(1,N);

for ii = 1:length(Eb_N0_dB)

 

% symbol generation

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

ipBit = rand(1,N*k,1)>0.5; % random 1's and 0's

bin2DecMatrix = ones(N,1)*(2.^[(k-1):-1:0]) ; % conversion from binary to decimal

ipBitReshape = reshape(ipBit,k,N).'; % grouping to N symbols having k bits each

ipGray = [sum(ipBitReshape.*bin2DecMatrix,2)].'; % decimal to binary

 

% Gray coded constellation mapping

ipDec = ind(ipGray+1)-1; % bit group to constellation point

ipPhase = ipDec*2*pi/M; % conversion to phase

ip = exp(j*ipPhase); % modulation

s = ip;

 

% noise

% -----

n = 1/sqrt(2)*[randn(1,N) + j*randn(1,N)]; % white guassian noise, 0dB variance

 

y = s + 10^(-Es_N0_dB(ii)/20)*n; % additive white gaussian noise

 

% demodulation

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

% finding the phase from [-pi to +pi]

opPhase = angle(y);

% unwrapping the phase i.e. phase less than 0 are

% added 2pi

opPhase(find(opPhase<0)) = opPhase(find(opPhase<0)) + 2*pi;

 

% rounding the received phase to the closest constellation

ipPhaseHat = 2*pi/M*round(opPhase/(2*pi/M))    ;

% as there is phase ambiguity for phase = 0 and 2*pi,

% changing all phases reported as 2*pi to 0.

% this is to enable comparison with the transmitted phase

ipPhaseHat(find(ipPhaseHat==2*pi)) = 0;

ipDecHat = round(ipPhaseHat*M/(2*pi));

 

% Decimal to Gray code conversion

ipGrayHat = map(ipDecHat+1); % converting to decimal

ipBinHat = dec2bin(ipGrayHat,k) ; % decimal to binary

 

% converting binary string to number

ipBinHat = ipBinHat.';

ipBinHat = ipBinHat(1:end).';

ipBinHat = str2num(ipBinHat).' ;

 

% counting errors

nBitErr(ii) = size(find([ipBit- ipBinHat]),2); % couting the number of errors

 

end

simBer = nBitErr/(N*k);

theoryBer = (1/k)*erfc(sqrt(k*10.^(Eb_N0_dB/10))*sin(pi/M)); % PSK 의 이론적 BER

 

close all; figure

semilogy(Eb_N0_dB,theoryBer,'bs-','LineWidth',2);

hold on

semilogy(Eb_N0_dB,simBer,'mx-','LineWidth',2);

axis([0 20 10^-5 1])

grid on

legend('theory', 'simulation');

xlabel('Eb/No, dB')

ylabel('Bit Error Rate')

title('Bit error probability curve for 16-PSK modulation')

 

위 코드를 돌려보면~ 아래 그림과 같이 16 PSK 결과가 나옵니다. 아래 BER 그래프에서 16 dB 이상에서 이론 값과 약간 달라 보이는 것은 시뮬레이션 횟수가 부족해서 BER 이 0 으로 나오는 경우가 생겨서 그런 것입니다. 고 dB 에서도 이론값과 같은 결과를 보이고 싶다면 심볼수인 N = 10^5 값을 좀더 늘려 보시기 바랍니다. 

 

위 코드를 약간 수정해서~질문하신 8 PSK 시뮬레이션을 해보죠~ 위 코드에서 빨간색 표시한 부분들만 아래와 같이 변경합니다.

 

M = 8; % constellation size

Eb_N0_dB = [0:15]; % multiple Es/N0 values

axis([0 15 10^-5 1])

title('Bit error probability curve for 8-PSK modulation')

 

시뮬레이션을 해보시면~ 다음과 같이 8 PSK 결과를 확인 할 수 있습니다.

 

다음과 같이 변경하면 QPSK 시뮬레이션이 되는거죠~

 

M = 4; % constellation size

Eb_N0_dB = [0:10]; % multiple Es/N0 values

axis([0 10 10^-5 1])

title('Bit error probability curve for QPSK modulation')