MATLAB BER simulation 관련 검색을 하다가 아래 주소의 글을 발견했습니다.

 

http://www.mathworks.com/matlabcentral/fileexchange/22316-communication-systems-reference-curves

 

위 파일에서 설명하는 내용은 아래와 같습니다.

 

- PSK and QAM over AWGN Channel (BER and SER)

- BPSK over Rayleigh fading channel (BER)

- Convolutional Coded BPSK over AWGN (BER)

 

일단 BER 이나 SER 시뮬레이션을 어떻게 해야 할 지 모르는 분들이 보기에 참 좋은 코드라는 생각이 들더군요.

 

위 주소에 들어가서 아래 그림과 같이 Download Submission 버튼을 누르면 파일을 다운로드 할 수 있습니다.

 

위 코드가 실행되기 위해서는 MATLAB 뿐만 아니라 Simulink, Communications Toolbox, Communications Blockset 이 필요 합니다.

 

다운로드 받은 Submit_v2.zip 파일을 압축해제 하고~ RunMe.m 파일을 실행시키면 아래 그림과 같이 웹 화면이 나옵니다.

 

웹 화면에서 뒤로 가기를 누르면 정상적으로 실행 안 되는 경우가 있으므로 이전 페이지의 내용을 다시 돌려 보고 싶다는 분들은 RunMe.m 파일을 다시 실행하고 들어가서 실행시키면 됩니다.

 

다양한 내용들이 있는데~ 저는 그 중에서 Rayleigh Flat Fading Channel Bit Error Rate curves 만 한번 실행해 봤습니다.

 

아래 그림과 같이 Simulink 파일이 실행되면서 시뮬레이션이 진행 됩니다.

 

시뮬레이션이 진행 될 때마다~ 아래 그림과 같이 BER 그래프에 결과가 하나씩 추가되며~

 

맨 마지막에는 다음과 같은 결과를 확인 할 수 있습니다.

 

BER 또는 SER 시뮬레이션을 해 보려고 하시는 엔지니어 분들이나 디지털 통신을 공부하시는 학생 분들에게 대단히 좋은 예가 될 것 같네요~


아래 포스팅에서 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 신고

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

  5. 2020.05.02 23:25

    혹시 qpsk가 awgn 채널말고 rayleigh, rician, exponentail channel 일 경우 에는 매트랩이 어떻게 변형되는지 알 수 있을까요?

  6. 2020.05.02 23:45


    BER_BPSK_Rayleigh.m


    % MATLAB Script for computing the BER for BPSK modulation in a
    % Rayleigh fading channel and compred to AWGN Channel

    % Clear all the previously used variables
    clear all;

    format long;

    % Frame Length
    bit_count = 10000;

    %Range of SNR over which to simulate
    SNR = 0: 1: 40;

    % Start the main calculation loop
    for aa = 1: 1: length(SNR)

    % Initiate variables
    T_Errors = 0;
    T_bits = 0;

    % Keep going until you get 100 errors
    while T_Errors < 100

    % Generate some random bits
    uncoded_bits = round(rand(1,bit_count));

    % BPSK modulator
    tx = -2*(uncoded_bits-0.5);

    % Noise variance
    N0 = 1/10^(SNR(aa)/10);

    % Rayleigh channel fading
    h = 1/sqrt(2)*[randn(1,length(tx)) + j*randn(1,length(tx))];

    % Send over Gaussian Link to the receiver
    rx = h.*tx + sqrt(N0/2)*(randn(1,length(tx))+i*randn(1,length(tx)));

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

    % Equalization to remove fading effects. Ideal Equalization
    % Considered
    rx = rx./h;

    % BPSK demodulator at the Receiver
    rx2 = rx < 0;

    % Calculate Bit Errors
    diff = uncoded_bits - rx2;
    T_Errors = T_Errors + sum(abs(diff));
    T_bits = T_bits + length(uncoded_bits);

    end
    % Calculate Bit Error Rate
    BER(aa) = T_Errors / T_bits;
    disp(sprintf('bit error probability = %f',BER(aa)));

    end

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

    % Finally plot the BER Vs. SNR(dB) Curve on logarithmic scale
    % Calculate BER through Simulation

    % Rayleigh Theoretical BER
    SNRLin = 10.^(SNR/10);
    theoryBer = 0.5.*(1-sqrt(SNRLin./(SNRLin+1)));

    % Start Plotting
    % Rayleigh Theoretical BER
    figure(1);
    semilogy(SNR,theoryBer,'-','LineWidth',2);
    hold on;

    % Simulated BER
    figure(1);
    semilogy(SNR,BER,'or','LineWidth',2);
    hold on;
    xlabel('SNR (dB)');
    ylabel('BER');
    title('SNR Vs BER plot for BPSK Modualtion in Rayleigh Channel');


    % Theoretical BER
    figure(1);
    theoryBerAWGN = 0.5*erfc(sqrt(10.^(SNR/10)));
    semilogy(SNR,theoryBerAWGN,'blad-','LineWidth',2);
    legend('Rayleigh Theoretical','Rayleigh Simulated', 'AWGN Theoretical');
    axis([0 40 10^-5 0.5]);
    grid on;

    계속 질문드려서 죄송한데 위의 매트랩에서 어떻게 고쳐야 QPSK 가 될까요

    • 남성 2020.05.02 23:59 신고

      위 코드중에서는 아래 주석 부분의 변복조 부분을 바꿔야 될 것 같네요.

      % BPSK modulator
      % BPSK demodulator at the Receiver

      현재 코드는 BPSK 로 되어 있으니가 QPSK 로 바꾸기 위해서는 복소수 형태로 만들어 줘야 되겠죠


      아래 주소에도 코드가 있는것 같은데 확인해 보시길~

      https://www.mathworks.com/matlabcentral/fileexchange/35964-qpsk-over-rayleigh-fading-channel

  7. 여울 2020.05.06 22:35

    안녕하세요. 바쁘실텐데 매번 답변해 주셔서 감사합니다. 기록된 이름은 다르지만 매번 질문드렸었어요 ㅎㅎ

    MATLAB을 배운지 얼마안되서 지금 채널에 따른 FUNCTION에 대해 학습하고 있는데,

    일전에 답변해주셔서 AWGN Channel과 Rayleigh Channel에 대해서 알 수 있었습니다.

    이번에는 Exponential channel에 대해서 문의드릴려고 하는데
    [송신단에서 신호를 쏘았을 때, 최단거리를 통과한 신호를 a, 벽에 한번 튕겨온 신호를 b, 벽에 세번 튕겨 온 신호를 c라고 가정하고
    이 때 a를 delay = 0라고 하고, b는 a보다 5ms 더 늦게 도착, c는 a보다 20ms 더 늦게 수신단에 도착했다고 하면,
    이 때의 tap개수는 모두 3개이고, tap delay = 0, 5ms, 20ms로 설정하여 exponential channel을 모델링
    (이렇게 하나의 신호를 송신단에서 보냈을 때 여러 갈래로 수신단에 도착하는 채널 상태를 multipath 채널)]

    위에 기재해주신 FUNCTION에서 어디에 어떻게 아래의 FUNCTION을 추가하였을때 Exponential channel FUNCTION이 되는 지 궁금합니다.
    BER of QPSK(With AWGH Channel)의 채널 생성부분에 입력을 하면 될 것 같기는 한데.. ㅠㅠ 아직 입문단계라 그런지 잘 모르겠습니다. ㅠㅠ

    그리고 아래의 함수는 수정없이 바로 사용하여도 될지 궁금합니다. 바쁘실텐데 매번 댓글로 답변 달아주셔서 감사합니다!

    function [ sd kmax ] = ExpModel( a )

    %---Frequency Selective Channel (Exponential Model)-----%

    kmax= ceil( 10*a);

    if kmax ==0
    sigma0 = 1;
    k = 0;
    VarK = 1;
    else
    beta = exp(-1/a);
    Var0 = ( 1 - beta ) / ( 1 - beta^(kmax+1) );
    k = (0:kmax)';
    VarK = Var0*beta.^k;
    end

    %standard deviation of exponential distribution
    sd = sqrt(VarK);

    %------< original codes >-----------%
    %
    % function [taps] = ExpChanTaps(sampRateMHz, delaySprdNsec)
    %
    % sampTimeNsec = 1000 / sampRateMHz;
    % if delaySprdNsec == 0
    % Kmax = 0;
    % vark = 1;
    % else
    % Kmax = ceil( 10 * delaySprdNsec/sampTimeNsec );
    % var0 = (1 - exp( -1*(sampTimeNsec)/delaySprdNsec ))/ ...
    % (1 - exp( -1*((Kmax+1)* sampTimeNsec)/delaySprdNsec ));
    %
    %
    % k = (0:Kmax)';
    % vark = var0 * exp( - k*sampTimeNsec/delaySprdNsec );
    % end

    %
    % stdDevReOrIm = sqrt(vark/2);
    % taps = stdDevReOrIm .* (randn(Kmax+1,1) + j*randn(Kmax+1,1));
    %

    • 남성 2020.05.07 01:58 신고

      저도 Frequency Selective Channel 은 MATLAB 통신 툴박스에 있는걸 사용했지, 만들어서 사용해 본적이 없어서 그닥 도움은 못 될 것 같네요.

      그래서 저도 잘은 모르지만...
      질문하신 Frequency Selective Channel 은 아시는 바와 같이 multipath 의 형태이므로 impulse 의 형태로 모델링이 될 것이고 신호의 sampling rate 와 delay 시간에 따라 계산해서 모델링이 되어야 할 것으로 보입니다.

      즉 이전 질문하신 코드 중에 h 라고 되어 있는 부분이 현재는 스칼라값이지만 벡터값으로 모델링이 되어야 하지 않을까 생각합니다. 그리고 h .*tx 의 형태에서 conv(h,tx) 또는 filter(h,1,tx) 의 형태로 채널이 신호에 반영되게 될것 같네요. 또 simulation 이 될 때마다 채널 h 의 delay 는 그대로 지만 각 impulse 의 파워는 변화하는 형태가 되도록 만들어야 할 것 같습니다.

      그리고 ExpModel, ExpChanTaps 함수를 사용해도 될지 질문을 주셨는데 음... 잘 모르겠지만 ExpChanTaps 함수에서 argument 로 sampRateMHz, delaySprdNsec 를 사용하는 걸 보면 뭔가 multipath 채널을 만드는것 같기는 하네요.

      제 답변은 잘 모르는 상태에서 그냥 도움이 될까해서 첨언한 것 뿐이니 참조만 하시길

+ Recent posts