아래 포스팅에서 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')

 


  1. 2014.04.17 15:33

    비밀댓글입니다

    • 남성 2014.04.17 18:23 신고

      네~ 방문해 주셔서 감사합니다. ^^ 열공하세요~

  2. 곰비 2014.05.12 00:34

    QPSK의 BER이론치를 구하는 공식을 알수잇을까요.??ㅠㅠㅠ

    • 남성 2014.05.12 08:02 신고

      PSK 의 이론적 BER 이라고 되 있는 부분의 코드를 보시면 수식을 아실 수 있을 겁니다.

  3. yj 2017.09.25 20:31

    BPSK는 해당 이론식에 적용되지 않는건가요?

    • 남성 2017.09.26 04:21 신고

      아니요 위 포스팅에서는 그냥 시뮬레이션 안 한 것뿐입니다. 똑같이 적용 됩니다

  4. 궁금해 2018.11.10 00:54

    이론값과 실제값이 차이나는게 심볼 수 때문이라고 하는데 왜 그런 차이가 나는건가요??
    오차가 생기는 원인이 궁금합니다ㅠㅠ

    • 남성 2018.11.10 21:43 신고

      확률적인 시뮬레이션이니까 시뮬레이션 샘플수가 많을수록 이론적인 값과 근사하게 되겠죠.

이번 포스팅에서는 AWGN (Additive white gaussian noise) 채널 환경에서의 QPSK Bits Error Rate(BER)에 성능 검증 simulation 에 대해 설명한다.

 

QPSK 및 잡음의 분산(σ2) 과 No(noise power spectral density) 의 관계에 대한 내용은 다음 포스팅을 참조 하기 바란다.




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


2011/03/08 - [통신] - [디지털 통신] BPSK, QPSK

 

QPSK 에서는 아래 식의 M =4 이다.

 

 

따라서 다음 식과 같이 정리 된다.

 

 

QPSK 시뮬레이션 시의 주의할 점신호의 파워를 1 로 정규화 해 줘야 하며 잡음 역시 real , imaginary 에 대해 각각 독립적이어야 한다.

 

이제 QPSK 시뮬레이션을 해 보자.

 

  

Command 창의 결과는 다음과 같다.

Eb/No= 0 BER: 0.078801

Eb/No= 2 BER: 0.0374294

Eb/No= 4 BER: 0.0124703

Eb/No= 6 BER: 0.00238721

Eb/No= 8 BER: 0.0001905

Eb/No= 10 BER: 3.07143e-006

 

 

위 그래프를 보면 이론적인 결과와 근사적으로 일치 함을 확인 할 수 있다.


  1. 남성 2012.11.25 18:47 신고

    이렇게 바꿔보세요.

    Symbols=((bits_v(:,1)*(-2)+1) + i*(bits_v(:,2)*(-2)+1))./sqrt(2);

  2. 과객 2013.06.18 23:52

    BPSK QPSK 정말 감사하게 잘 실행되었습니다.
    다만 QPSK에서는 비트 2개를 하나는 리얼로 하나는 이미지로 생각해서 처리하셨는데 8PSK에서는 비트 3개라 어떻게 해야 할지 모르겠네요 8PSK에서는 어떻게 해야 할까요?

    • 남성 2013.06.19 06:45 신고

      8 psk 에서도 3 비트를 한 점의 symbol 로 mapping 하시면 됩니다.

      아래 주소의 그림처럼 말이죠.

      http://commons.wikimedia.org/wiki/File:8PSK_Gray_Coded.svg


      방문해 주셔서 감사합니다. ^^

    • 과객 2013.06.20 20:11

      감사합니다.

  3. 답변부탁드릴게요 2013.11.17 23:06

    도움 감사합니다 하지만 실행이 되지 않아서 질문드립니다
    답변해주시면 정말 감사하겠습니다.

    clc
    clear
    close all

    N_bits=70000000; %생성할 비트 수 / 2
    Eb_No_dB=0:2:10; %Eb_No dB scale
    BER_buffer=zeros(size(Eb_No_dB)); %BER 저장할 버퍼


    for n=1:length(Eb_No_dB)
    Eb_No_ral_scale=10^(Eb_No_dB(n)/10); %Eb/No real scale 값 반환
    sigma_v=sqrt(1./(4*Eb_No_ral_scale)); %잡음의 표준 편차 값 계산

    bits_v=randi([0 1],N_bits,2); %비트 생성

    Symbols=(bits_v(:,1)*(-2)+1 + i *(bits_v(:,2)*(-2)+1))./sqrt(2); %symbol mapping & normalizat

    noise_v=(randn(N_bits,1)*sigma_v)+ i * (randn(N_bits,1)*sigma_v); % complex AWGN 생성

    tx_signal = Symbols * noise_v; % 신호*잡음

    demapped_bits=[real(tx_signal) imag(tx_signal)] < 0; %수신단 demapping

    BER_buffer(n)=sum(bits_v(:) ~= demapped_bits(:)) ./ (N_bits*2); %BER calculation

    fprintf('Eb/No=%g BER: %g\n' ,Eb_No_dB(n) ,BER_buffer(n))
    end

    %%이론적인 QPSK BER 계산 Using communication toolbox

    Eb_No_theory = 0:0.05:10;
    Ber_value_theory = berawgn(Eb_No_theory,'psk',4,'nondiff');
    %% 그래프그리기

    figure(1),semilogy(Eb_No_dB,BER_buffer,'b+',Eb_No_theory,Ber_value_theory,'r--'),x;xlabel('Eb/No [dB]'),ylabel('BER')
    grid.axis([0 10 0.000001 1]), legend('Simulated result','Theoretical BPSK result')


    제시해주신 코드와 똑같이 작성하려고 노력하였습니다,
    하지만 아래와 같은 오류가 나와서 질문드립니다.

    ??? Error using ==> randn
    Out of memory. Type HELP MEMORY for your options.

    Error in ==> QPSK_BER at 18
    noise_v=(randn(N_bits,1)*sigma_v)+ i * (randn(N_bits,1)*sigma_v); % complex
    AWGN 생성

    어떤 문제인지 알려주시면 정말 감사하겠습니다!

    • 남성 2013.11.17 23:59 신고

      randn() 에서 메모리 에러 나는것은 N_bits 를 너무 크게 줘서 그렇습니다.

      아래 세 부분이 문제더군요.

      다음과 같이 바꾸세요.

      N_bits=7000000; %생성할 비트 수 / 2, 너무 크게 주면 컴퓨터뻗어요.

      tx_signal = Symbols + noise_v; % 신호 + 잡음, 곱하기가 아닙니다.


      % 아래 부분에서는 단순 오타가 좀 있더군요.
      figure(1),semilogy(Eb_No_dB,BER_buffer,'b+',Eb_No_theory,Ber_value_theory,'r--');xlabel('Eb/No [dB]'),ylabel('BER')
      grid on , axis([0 10 0.000001 1]), legend('Simulated result','Theoretical BPSK result')

    • 답변부탁드릴게요 2013.11.18 01:06

      감사합니다! 정상작동됩니다!

    • 남성 2013.11.18 01:39 신고

      네~ 열심히 공부하세욧 !! 방문해 주셔서 감사합니다. ^^

  4. 답변좀 부탁드려요 2013.11.20 17:06

    안녕하세요~ 과제 때문에 들렀다가 오류가 떠서 질문드립니다..ㅠ
    Eb_N0_theory=0:0.05:10 ; 이 부분에서요
    다음과 같은 오류가 뜹니다..
    The expression to the left of the equals sign is not a valid target for an assignment.
    Eb_N0 대신 혹시나 해서 Eb_No로 바꿔서 입력해보아도 같은 오류 메시지가 뜨는데 어떻게 해야 하나요?

    • 남성 2013.11.20 18:30 신고

      Eb_N0_theory=0:0.05:10 ; 부분에 코드상의 문제는 없습니다.
      다른 부분에 문제가 있는것 같은데요.

  5. 119 2014.04.01 17:59

    음.. ㅠㅠ 저기 tx_signal 차원이 [N_bits,1]맞죠?? symbol하구 noise 왜 더하는데 둘이 자꾸 차원이 안맞다구 나오는지 ㅠㅠ

  6. 119 2014.04.01 19:58

    저는 sigma(n)으로 고쳐서 곱해주니까다시됏어요.. 그리구 demappping이랑 BER계산 어떻게 하신건가요ㅠㅠ??

    • 남성 2014.04.01 21:23 신고

      demapping은 0 보다 작으면 1 그렇지 않으면 0 으로 demapping 한 거구요. 이론적 ber 게산은 berawgn 이라는 통신 툴박스의 함수를 사용했습니다. 위에 보시다 시피 코드는 저게 다에요.

  7. 2014.04.16 23:14

    비밀댓글입니다

    • 남성 2014.04.17 01:06 신고

      아래 포스팅에 답변 합니다.

      http://iamaman.tistory.com/1251

  8. 공돌공돌 2014.06.06 22:06

    clc
    clear all
    close all

    N_bits=7000000; %생성할 비트 수/2
    Eb_No_dB=0:2:10; %Eb/No dB scale
    BER_buffer=zeros(size(Eb_No_dB)); %BER 저장할 버퍼

    for n=1:length(Eb_No_dB)
    Eb_No_ral_scale=10^(Eb_No_dB(n)/10); %Eb/No real scale 값 변환
    sigma_v=sqrt(1./(4*Eb_No_ral_scale)); %잡음의 표준 편차 값 계산

    bits_v=randi([0 1],N_bits,2); %비트 생성

    Symbols=(bits_v(:,1)*(-2)+1 +i*(bits_v(:,2)*(-2)+1))./sqrt(2); %symbol mapping & normalization

    noise_v=(randn(N_bits,1)*sigma_v) +i*(randn(N_bits,1)*sigma_v); %complex AWGN 생성

    tx_signal=Symbols+noise_v; %신호 + 잡음

    demapped_bits=[real(tx_signal) imag(tx_signal)]<0; %수신단 demapping

    BER_buffer(n)=sum(bits_v(:)~=demapped_bits(:))./(N_bits*2); %BER calculation

    fprint('Eb/No=%g BER:%g\n',Eb_No-dB(n), BER_buffer(n))
    end

    % 이론적인 QPSK BER 계산
    Eb_No_theory=0:0.05:10;
    Ber_value_theory=berawgn(Eb_No_theory,'psk',4,'nondiff');

    % 그래프 그리기
    figure(1),semilogy(Eb_No_dB,BER_buffer,'b+',Eb_No_theory,Ber_value_theory,'r--');xlabel('Eb/No[dB]'), ylabel('BER')
    grid on, axis([0 10 0.000001 1]), legend('Simulated result','Theoritical QPSK result')

    이렇게 했는데 왜 BERTool에서 run 누르면 invalid simulation M-file or model 이라는 에러가 나오면서 그림이 안그려지죠? ㅜㅜ 제발 알려주세요 ㅜㅜㅜㅜ

  9. 남성 2014.06.07 06:07 신고

    위 코드에서는 25 번째 줄에 오타가 몇개 있더군요. 아래와 같이 하면 정상적으로 동작 합니다.

    fprintf('Eb/No=%g BER:%g\n',Eb_No_dB(n), BER_buffer(n))

    위 질문에서 bertool 이라고 얘기하시는데... 위 코드를 m 파일로 작성하시고 bertool 에서 monte carlo simulation 하셨다는 건가요?

    • 남성 2014.06.07 10:34 신고

      먼가 잘 못 이해하고 계신거 같은데... 위 코드는 그냥 m 파일로 돌리면 되는건데요? bertool 안 써도 되는겁니다.

    • 공돌공돌 2014.06.07 10:35

      25번째줄 말씀하신대로 수정하고 저장한 후에 bertool 켜서 monte carlo simulation 했는데 또 invalid simulation M-file or model 이라는 에러가 뜨네요 ㅜㅜ 대체 왜 이러는지 모르겠습니다 ㅜ 위에 님이 포스팅하신 것처럼 이론적인 QPSK 그래프와 코드 시뮬레이션 한 값과 비교하는 그래프를 봐야하는데 그게 안뜹니다 ㅜㅜ

    • 2014.06.07 10:38

      비밀댓글입니다

    • 남성 2014.06.07 11:51 신고

      bertool을 사용하신다면 변수명 등을 bertool 에 맞게 써 줘야 하는것 같더군요.

  10. 공돌이 2014.06.08 15:36

    안녕하세요 자료를 찾아보다가 오게 되어서 질문드립니다 ㅠ

    제가 지금 설계프로젝트를 하고 있는데요 32QAM을 써서 통신시스템 설계를 진행하고 있는데 awgn만 첨가했을 때

    bertool이랑은 전혀 다른 결과가 나와서 문의드립니다.

    들어가는 비트 Eb - 채널코딩 한 후 Ec - 변조 한 후 Es로 해서

    N0=10.^(EcN0dB(snrindex)/10);
    gain=sqrt(N0/2);
    noise=randn(200,1)*gain

    으로 했는데 bertool보다 훨씬 더 좋게 나오네요,,, 이러면 안되는데,,, awgn 채널을 잘못구성한건가요??

    • 남성 2014.06.08 17:08 신고

      채널 코딩도 들어갔으면 파워 계산 과정에서 채널 코드 레이트에 대한 수치가 반영 되어야 합니다. 채널 코딩시 encoder 의 입력 비트에 비해 출력비트가 더 많이 나오므로 출력의 각 비트당 파워는 줄어들게 되겠죠. 이를 정확히 반영하셔야 정확한 시뮬레이션이 될 겁니다. 이를 제대로 반영안했다면 이론 값보다 당연히 좋게 나오겠죠.

  11. 멍청이 2015.07.22 20:32

    질문이 있는데 qpsk 에서 위에 코드에서 Symbols=(bits_v(:,1)*(-2)+1 + i *(bits_v(:,2)*(-2)+1))./sqrt(2); 루트 2 이걸로 나눠주는데 이건 왜 나눠 주는거죠???

    • 남성 2015.07.23 04:18 신고

      위의 포스팅에서 적어 놨지만 신호의 파워, 여기서는 심볼의 파워를 1로 맞춰주는겁니다.
      Symbols = ±1/sqrt(2) ± j · 1/sqrt(2) 이므로 절대값이 1 이라는것을 확인 할 수 있습니다. 심볼 파워를 1로 정규화한 상태에서 잡음의 파워 즉 분산을 계산해서 넣어주는 겁니다.

+ Recent posts