오늘은 AWGN (Additive white gaussian noise) 채널 환경에서의 BPSK Bits Error Rate(BER)에 대해 설명한다. 일단 BPSK 가 뭔지 잘 모르는 사람은 다음 페이지를 참조 하기 바란다.


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

 

위 페이지를 참조하면 bit 를 BPSK 심볼로 mapping 하는 방식에 대해 알았을 테고, 이제 잡음을 생성하는 방법을 알아 보자. 잡음을 생성하기 위해서는 몇 지 통신 관련 수식에 대해 알아야 한다. 보통 BER 그래프를 그릴 때는 SNR 축에 따른 그래프가 아닌 Eb/No 에 따른 그래프를 그리곤 한다. 따라서 BER 시뮬레이션을 정확히 하려면 Eb/No 에 따른 잡음을 생성할 수 있어야 하겠다.

 

잡음의 분산(σ2) 과 No(noise power spectral density) 는 다음과 같은 관계가 있다.

 

 

주의 : 본 식에서는 oversample 이나 channel coding 에 대해 고려 하지 않은 것이다.

 

 

 

위 식에서 각 기호의 의미는 아래와 같다.

 

  • Es : Signal energy

 

  • Eb : Bit energy

 

  • N0 : Noise power spectral density

 

  • M : PSK modulation order

 

 

시뮬레이션의 Es 는 1 로 정규화 한다. 그럼 잡음의 분산 값은 다음과 같은 식이 된다.

 

 

BPSK 의 경우 성좌도의 경우의 수가 2 이므로 M = 2 이다. 그럼 이제 최종적으로 BPSK 시뮬레이션시의 잡음 분산은 다음과 같이 된다.

 

 

이제 BPSK 의 BER 시뮬레이션을 해 보자.

 

MATLAB 코드 내용은 다음과 같다.


2014-07-03 추가 내용 

캡쳐 화면을 올렸더니 오타로 고생하시는 분들이 너무 많더군요. Github 에 파일 올립니다. 

파일 다운로드 링크 입니다. 

https://github.com/ssgkd/MATLAB-BPSK-AWGNSimulation/archive/master.zip

 

위 코드를 실행 했을때의 command 창 결과는

 

Eb/No= 0 BER: 0.0784773

Eb/No= 2 BER: 0.0374897

Eb/No= 4 BER: 0.0125513

Eb/No= 6 BER: 0.00238729

Eb/No= 8 BER: 0.000199714

Eb/No= 10 BER: 4.85714e-006

 

이며 그래프는 다음과 같이 나온다.

 

 


 

이론값과 시뮬레이션 값이 거의 일치하는 것을 확인 할 수 있고, 10 dB 일때 약간 값이 다른 것은 시뮬레이션 bits 수를 늘리면 해결 된다. 그리고 위 코드에서 이론적인 BER 값을 계산할 때는 communication toolbox berawgn() 함수를 이용하였다.



  1. 이전 댓글 더보기
  2. 흐미 2014.04.13 14:56

    좋은자료 감사합니다. 그런데 왜 저는 simulated된 BER값이 일정하게 나올까요? 어떤게 잘못됬을까요?

    • 남성 2014.04.13 15:11 신고

      방문해 주셔서 감사합니다. ^^
      일정하게 나오다는건 0.5 정도로 나온다는 거죠? 0.5 로 나온다는건 시뮬레이션이 잘못 됐다는겁니다. 오타일 확률이 클것 같네요.

    • 흐미 2014.04.13 20:57

      계속 찾아보려고 해도 오타가 보이지 않네요ㅜ
      한번 봐주실 수 있으세요?
      clc
      clear
      close all

      N_bits=7000000;
      Eb_No_dB=0:1:10;
      BER_buffer=zeros(size(Eb_No_dB));


      for n=1:length(Eb_No_dB)
      Eb_No_real_scale =10^(Eb_No_dB(n)/10);
      sigma_v=sqrt(1./(2*Eb_No_real_scale));

      bits_v=randi([0,1],N_bits,1);
      Symbols=bits_v*(-2)*1;

      noise_v=randn(N_bits,1)*sigma_v;

      tx_signal=Symbols + noise_v;

      demapped_bits= tx_signal < 0;

      BER_buffer(n)=sum(bits_v~= demapped_bits)./N_bits;

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

      %% Theory BPSK with BER
      Eb_NO_theory = 0:0.05:10;
      Ber_value_theory=berawgn(Eb_NO_theory,'psk',2,'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','Theoretical BPSK result')

      덕분에 좋은 자료로 디지털 통신 이해가 되고 있습니다 감사합니다

    • 남성 2014.04.13 22:32 신고

      다 하셨는데 딱 한군데 틀리셨네요.
      Symbols=bits_v*(-2)*1; 를 Symbols=bits_v*(-2)+1;
      로 바꾸세요. 곱하기 1을 더하기 1로 바꾸시면 됩니다.
      바꾸고 돌려보니 잘 되네요.

    • 흐미 2014.04.13 22:34

      정말 감사합니다 ^^ 진짜 많이 배우고 갑니다.
      좋은 자료 정말 감사합니다!

    • 남성 2014.04.13 22:37 신고

      학교 숙젠가봐여? 도움 되셨다니 다행이네요.

    • 2014.04.13 23:07

      비밀댓글입니다

  3. 급해요 2014.04.15 13:29

    close all

    N_bits=200; %비트생성
    Eb_No_dB=0:1:10; %0부터 10까지 1씩지점을 체크하겠다.dB스케일로
    BER_buffer=zeros(size(Eb_No_dB)); %저장할 버퍼

    for n=1:length(Eb_No_dB)
    Eb_No_ral_scale=10^(Eb_No_dB(n)/10); %Eb/No값스케일변환
    sigma_v=sqrt(1./(2+Eb_No_ral_scale)); %잡음분산하려고 표준편차값 계산
    bits_v=randi([0 1],N_bits,1); %비트생성
    Symbols=bits_v*(-2)+1; %심볼 매핑
    noise_v=randn(N_bits,1)*sigma_v; %AWGN생성

    tx_signal=Symbols + noise_v; %원신호에 잡음더해주기
    res_bits=tx_signal < 0 %수신된것
    BER_buffer(n)=sum(bits_v ~= res_bits) ./ N_bits; %버퍼값 계산
    fprintf('Eb/No= %g BER: %g\n',Eb_no_dB(n), BER_buffer(n))
    end

    Eb_No_theory=0:0.05:10;
    Ber_value_theory=berawgn(Eb_No_theory,'psk',2,'mpmdoff');

    figure(1),semilogy(Eb_No_dB,BER_buffer,'b+',Eb_No_theory,Ber_value_theory,'r--'),xlabel('Eb/No [dB]'),ylabel('BER')
    grid,axis([0 10 0.000001 1]), legend('계산결과','이론적인 BPSK결과') %이론적값은 berawgn()함수사용

  4. 급해요 2014.04.15 13:30

    위에처럼 코딩했는데요
    Run하면 커맨드창에

    Undefined function 'Eb_no_dB' for input arguments of type 'double'.

    Error in Untitled2 (line 17)
    fprintf('Eb/No= %g BER: %g\n',Eb_no_dB(n), BER_buffer(n))

    나옵니다

  5. 급해요 2014.04.15 13:30

    그리고 그래프나오게 어떻게하죠??

  6. 부탁드려요 2014.05.08 11:09

    암만봐도어느부분이틀린건지알수가없네요 ㅠㅠ 한번만 확인해주시면 안될까요
    clc
    clear
    close all


    N_bits=7000000; % 생성할 비트 수
    Eb_No_dB=0:2:10; % Eb/N0 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/N0 real scale 값 변환
    sigma_v=sqrt(1./(2+Eb_No_ral_scale)); %잡음 표준편차 값 계산

    bits_v=randi([0 1],N_bits,1); % 비트 생성
    Symbols = bits_v*(-2)+1; % symbol mapping

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

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

    demapped_bits=tx_signal < 0; % 수신단 demapping

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

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

    %% 이론적인 BPSK BER 계산 Using Communication Toolbox

    Eb_No_theory=0:0.05:10;
    Ber_value_theory=berawgn(Eb_No_theory,'psk',2,'nondiff');

    %% 그래프 그리기

    figure(1), semilogy(Eb_No_dB,BER_buffer,'b+',Eb_N0_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')

    에러가 이렇게 뜨는데 봐도봐도 안보이네요 ㅠㅠ
    ??? Undefined function or method 'berawgn' for input
    arguments of type 'double'.

    Error in ==> BPSK at 32
    Ber_value_theory=berawgn(Eb_No_theory,'psk',2,'nondiff');

    공부하는데 정말 큰 도움이 되고 있습니다 고맙습니다~

    • 남성 2014.05.08 20:01 신고

      두군데 고치시면 되겠네요. 두군데 다 오타네요.


      아래 부분에서는 더하기를 곱하기로 고침

      sigma_v=sqrt(1./(2 * Eb_No_ral_scale)); %잡음 표준편차 값 계산



      아래 부분에서는 Eb_N0_theory 를 Eb_No_theory 로 고침, No의 O 를 소문자로

      figure(1), semilogy(Eb_No_dB,BER_buffer,'b+',Eb_No_theory,Ber_value_theory,'r--'),xlabel('Eb/No [dB]'),ylabel('BER')



      위 코드들은 수정한 코드이니 그대로 사용하시면 됩니다.

  7. ㅠㅠ 주무시겟죠 2014.05.12 01:16

    demapped_bits=tx_signal<0; % 수신된 비트신호
    이렇게 생각했는데 이게 왜 수신 신호인지 모르겟어요 ㅠㅠ

    • 남성 2014.05.12 08:01 신고

      mapping 을 할때 비트 0 을 1로 비트 1을 -1 로 mapping 을 하기 때문에 수신 신호가 0 보다 크면 비트 0으로 0보다 작으면 비트 1로 mapping 하는 겁니다.

  8. 2014.11.11 01:46

    비밀댓글입니다

    • 남성 2014.11.11 13:11 신고

      방문해 주셔서 감사합니다. ^^ 수정했습니다.

  9. 디시설 2014.11.16 23:32

    방금 위에서 물어 봣는데 그럼 혹시 저코드를 사진을 이용해 영상처리를 해서 잡음이 많은 거랑 적은걸 표현을 어떻게 하나요 ㅜㅜ

    • 남성 2014.11.17 01:25 신고

      사용하는 영상이 몇비트로 구성되어 있는지는 모르겟지만 영상의 값들을 비트로 표현한 후에 잡음을 섞어서 imshow() 와 같은 함수로 이미지를 보여 주면 잡음이 섞인 정도에 따라 비교가 될것 같은데요.

  10. 급한질문 2014.11.29 12:16

    혹시 이 코드에서 엔코더와 디코더를 넣을려면 어떻게해야하는지 알수잇을까요?

    • 남성 2014.11.29 18:22 신고

      어떤 인코더와 디코더를 말씀 하시는지... 위 코드의 시작과 끝이 비트 니까 각각의 외부에 넣으면 되겠네요.

  11. 서영화 2014.12.03 18:42

    송신에 Raised cosine 수신에 matched filter(square root RC로)추가하려면 어디에 너야할까요??

    • 남성 2014.12.03 21:29 신고

      송신 필터는 송신단 끝단에 넣으면 되고 수신 필터는 채널 통과후 수신 첫단에 넣으면 됩니다.

  12. 2014.12.10 18:17

    비밀댓글입니다

    • 남성 2014.12.11 18:12 신고

      아래 주소 참조하시면 쉽게 할 수 있을것 같네요.

      http://www.mathworks.com/matlabcentral/fileexchange/37455-bpsk-modulatian-demodulation-by-using-simulink

      Communication system toolbox 가 있어야 할 겁니다.

    • 2014.12.12 00:58

      비밀댓글입니다

    • 남성 2014.12.12 08:49 신고

      위에 그래프 그리는 코드도 들어가 있는데요.

    • 2014.12.13 01:40

      비밀댓글입니다

    • 남성 2014.12.13 02:30 신고

      아래 글에 작성한 바와 같이 simulink 에는 to workspace 라는 블록이 있습니다.

      http://iamaman.tistory.com/485

      simulink 에서 BER 계산 결과와 Eb/N0 값을 to workspace 블록을 이용하여 MATLAB workspace 에 저장한 후에 그래프를 그리면 될것 같네요.

    • 도와주세요 ㅠㅠ 2014.12.16 12:57

      감사합니다. 별거 아닌데 너무 어렵게 생각했네요. 많이배워갑니다

    • 남성 2014.12.16 13:21 신고

      도움되었다니 다행이네요. 방문해 주셔서 감사합니다. ^^

  13. 디시설 2014.12.20 15:37

    clc
    clear
    close all


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


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

    bits_v=randi([0 1],N_bits,1); % 비트 생성
    Symbols=bits_v*(-2)+1; % symbol mapping

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

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

    demapped_bits= tx_signal < 0 ; % 수신단 demapping

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

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

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

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

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

    'berawgn'은(는) 유형 'double'의 입력 인수에 대해 정의되지 않은 함수입니다.

    오류 발생: Untitled (line 33)
    Ber_value_theory=berawgn(Eb_N0_theory,'psk',2,'nondiff');

    이렇게오류가 나는데 무슨 문제인가요 ㅜㅜ

    • 남성 2014.11.16 20:02 신고

      Communications System Toolbox 가 있는건지 확인해 보세요.

  14. dduss 2015.01.14 13:42

    디지털공부하는데 정말 많은 도움이 되고 있습니다! 이글 보고 BER 코드를 짜고있는데 송수신단에 각각 다른 필터를 사용해 성능비교를 하려고 하는데 raised cosine filter를 추가할때 oversample을 해서인지 demapped bit 와 bits_v가 사이즈가 맞지않아 연산이 수행되지 않는것 같습니다.. 어떻게 해결해야할까요? ㅠㅠ

    • 남성 2015.01.14 23:03 신고

      송신단에서 oversample 을 했으면 수신단 필터링 후에 downsample 을 해야 겠죠. 그리고 필터 delay 때문에 송신단에서는 0 을 delay 만큼 추가를 해 줘야 송신한 모든 데이터를 수신 할 수 있을 겁니다.

    • dduss 2015.01.15 13:29

      그럼 oversampling만큼 upsampling 하고 필터링 후에 downsampling하면 되는 건가요?? 필터를 추가하면 noise를 더할 때도 비트수가 안맞는 것 같아요.. 공부를 더 많이 해야겠네요 ㅠㅠ

    • 남성 2015.01.15 20:58 신고

      upsampling 하고 필터링 하는게 oversampling 이라고 하고요. oversampling 한 다음에 잡음을 섞어야 하고 잡음을 섞을때는 sample 에 섞어야 하니까 sample 에 따른 잡음 분산을 다시 계산해야겠죠. 그리고 수신단에서 필터링 한다음에 downsampling 을 하는 겁니다.

    • dduss 2015.01.15 22:18

      개념을 제가 잘못이해하고있었네요! 정말 감사합니다. 본문도 댓글도 정말많은 도움이 되었어요~

    • 남성 2015.01.16 08:19 신고

      도움되었다니 다행이네요. ^^

  15. 온지 2016.04.19 20:33

    통신이론에 대해 공부하는 학생인데요 BER그래프를 그릴때 X축을 SNR이 아닌 Eb/No으로 해야되는건지 정확한 이유좀 알수있을까요?ㅠㅠㅠ

    • 남성 2016.04.19 21:30 신고

      SNR 축으로 그릴 수는 있습니다. 어떤걸 확인하는게 목적이냐에 따라 다르겠죠. 보통 BER 은 비트의 에너지에 따른 BER 성능을 구하려고 하는거죠. BER 을 구한다 함은 변조에 따른 심볼당 비트수 차에 의한 신호의 파워 차이를 없애고 순수 비트 에너지에 따른 BER 성능을 구하는게 타당 할 겁니다. BPSK 와 QPSK 를 Eb/No 를 x 축으로 해서 그리면 성능이 같지만 SNR 을 축으로 해서 BER 을 구한다면 성능이 다르겠죠? 이때 Eb/No 를 기준으로 한 BER 을 보고 BPSK 와 QPSK 는 BER 성능이 같다고 해석하는게 맞을까요? 아니면 SNR 축에 따라 그려진 BER 을 보고 BPSK 가 QPSK 보다 BER 성능이 좋다라고 해석하는게 맞을까요?

  16. 임명호 2016.11.28 20:22

    안녕하세요. 통신을 전혀 모르는 직장인인데, 아는 지인 분께 통신 관련 논문으로 BPSK랑 QPSK BER을 시뮬레이션 해달라고 하셔서 이리저리 헤메다가 질문 올려봅니다. (통신에 대한 지식이 전무합니다 ㅠㅠ...엉뚱한 소리라도 이해좀 부탁드립니다)

    1. Eb/N0 가 왜 Matlab 소스에서 0에서 10까지로 2간격으로 되어있는지요? 단순히 축을 설정하기 위한 것인가요?
    단순히 축 설정을 위한 것이라면 X축이 같으므로 Y 값이 다 똑같은 결과가 나오는 것이 아닙니까?
    제가 참조하는 논문에는 Eb/N0 = (S/Rb) / (N/W) 라고 되어있습니다. (S: 신호전력 Rb:전송속도 N:잡음전력 W:대역폭)

    밑에는 제가 올린 지식인 질문입니다. 해당 내용에 논문에 대한 참조 데이터가 있습니다. (원하신다면 논문을 보내드리겠습니다.)
    http://kin.naver.com/qna/detail.nhn?d1id=1&dirId=108&docId=265206138

    • 남성 2016.11.28 21:21 신고

      ber 시뮬레이션은 eb/n0 에 따른 비트 에러율을 구하는 겁니다 따라서 eb/n0 값을 시뮬레이션 때마다 조절해줘야 합니다. 즉 S값을 달리하면서 시물레이션 하는 거죠

    • 2016.11.28 23:45

      비밀댓글입니다

    • 남성 2016.11.29 00:13 신고

      제가 질문자님의 의도를 제대로 이해한건지 모르겠는데... 비콘과 발리스의 비트레이트가 다르다고해서 비콘은 bpsk로 변조하고 발리스는 qpsk로 변조해서 에러 레이트 비교한 다음에 둘중에 어느하나가 더 좋다라고 결론내려는건가요?

  17. 임명호 2016.11.30 16:02

    안녕하세요 우선 답변주셔서 감사드립니다. 제가 글을 명확히 안 써서 혼돈이 온 점 죄송스럽게 생각합니다.

    비콘
    발리스에
    대해 각각
    Bpsk ber 시뮬레이션
    Qpsk ber 시뮬레이션
    을 하여 성능을 비교해보려합니다.

    1. Bpsk ber
    비콘 vs 발리스

    2. Qpsk ber
    비콘 vs 발리스

    이런 시뮬레이션을 하려합니다.

    • 남성 2016.11.30 16:13 신고

      해당 내용에 대해 위 포스팅 대로 시뮬레이션을 하면 성능차이가 없을거 같네요

  18. 임명호 2016.12.28 10:19

    남성님 안녕하세요. 1달만에 다시 질문 드립니다.
    현재 BPSK와 QPSK 시뮬레이션을 하려하는데
    현재 SNR 값을 구하였습니다. 제가 조사한 바에 의하면 SNR은 아날로그 통신에서 주로 사용하고
    디지털통신에서는 SNR[DB]을 정규화한 Eb/N0[DB] 로 변환하는걸로 알고있습니다.
    최종적으로 정리한 공식을 찾아보니 아래와 같습니다.

    SNR=10log(S/N) [DB]
    Eb/N0 = (S/N) * (W/R)

    여기서,
    S: 신호전력 [V/m]
    N: 노이즈전력[V/m]
    W : 대역폭 [HZ]
    R : 비트전송속도 [bits/s]

    현재 저는 S,N,R 값을 알고있는데 W 값을 모르겠습니다.
    제가 가지고 있는 정보는

    샘플링 시간 : 4.4325 * 10^(-8)
    심볼 시간 :1.773 * 10^(-6)
    전송 주파수 4.23 MHZ 입니다.

    그리고 BPSK와 QPSK로 시뮬레이션을 하려하는데 제가 W 값을 위 정보를 바탕으로 구할 수 있는지요?

    • 남성 2016.12.28 10:32 신고

      아래 주소 보시면 심볼레이트와 대역폭 사이의 관계를 알 수 있습니다.

      http://www.ktword.co.kr/abbr_view.php?m_temp1=4185

  19. zzzxxxccc 2017.11.25 14:21

    BPSK에서 unipolar signal과 bipolar signal 을 동시에 그래프상에서 보여주고 싶은데 어떻게하나요??

    • 남성 2017.11.25 14:31 신고

      두 시그널에 따른 ber 성능을 비교하고 싶다는 건가요? 아니면 그래프 두개를 동시에 그리는 방법을 물어보시는 건가요?

  20. here 2018.11.25 12:46

    안녕하세요! AWGN에서 BER값을 내기 때문에 BER값을 구할때 계속 다른값이 나오잖아요!
    이걸 백번돌려서 나온 백개들의 평균을 구해서 그래프를 그릴려고 합니다! 이 때 이걸 백번돌리는 함수는 어떻게 만들면 될까요?

    • 남성 2018.11.25 15:27 신고

      백번 돌려서 평균을 구하려면 당연히 버퍼를 만들어서 저장을 하는 형태로 만들어야 되겠죠.

  21. 안녕하세요 2018.12.11 16:35

    제가 QPSK BER curve를 plot하려고 하는데요
    QPSK지만 1bit에 대해서만 ber curve를 그리려고 합니다. 그래서 matlab코드를 이렇게 짜봣는데요
    clc;
    clear all;
    %QPSK BER
    %1->√(2Eb/Tb)*cos(2πfct)이고 0->-√(2Eb/Tb)*cos(2πfct)로 mapping 된다.
    Eb=1;
    sample_num=100000;%각 Eb/No당 샘플갯수
    Eb_over_No_dB=-2.5:0.1:7.5;
    Eb_over_No=10.^(Eb_over_No_dB./10);
    No=1./Eb_over_No.*Eb;
    Bit_error=zeros(1,length(Eb_over_No_dB));%각 Eb/No값마다 bit error값을 저장할 행렬
    for k=1:length(Eb_over_No_dB)% Noise variance를 바꿈.즉, Eb/No값을 -2.5dB에서 7.5dB까지 바꿈. 각각의 Noise variance에 대하여 BER을 찾는다.
    A=randi([0,1],1,sample_num);%랜덤한 0,1 sequence를 sample number만큼 생성.즉, binary데이터 열을 만든것.
    No_over_2=No(1,k)./2;%Noise Variance No/2
    %QPSK oddbit에 대한 receiver 출력구하기.평균이 ±√Eb variance가 No/2인
    %Random Varable(=Output)
    QPSK_odd=(A*2-1)*sqrt(Eb);
    Noise=sqrt(No_over_2)*(randn(1,sample_num));
    Output=QPSK_odd+Noise;


    %BER계산
    for i=1:sample_num
    if(A(1,i)==1) %송신된 Binary data가 1일때
    if( Output(1,i)<0)%receiver 출력이 0보다 작으면 -----
    Bit_error(1,k)=Bit_error(1,k)+1;%-->error
    end
    else if(A(1,i)==0) %송신된 Binary data가 1일때
    if( Output(1,i)>0)%receiver 출력이 0보다 크면 -----
    Bit_error(1,k)=Bit_error(1,k)+1;%error
    end
    end
    end
    end
    end

    BER=Bit_error./sample_num.*100;
    BER_Theory=1/2*erfc(sqrt(Eb_over_No));
    semilogy(Eb_over_No_dB,BER);

    방법에 대해 설명해보자면 binary data와 receive된 data를 비교해서 1을 보냈을때 receive signal이 0보다 작으면 error이고
    0을 보냈을때 receive signal이 0보다 크면 error가 나오는 이런 식으로 짰는데 BER이 너무 높게 나옵니다.
    이론값과 100배이상 나는데 어디가 잘못된 걸까요?

    • 남성 2018.12.11 17:42 신고

      위 코드 링크 넣어놨는데 다운로드 받아서 돌려 보시면서 비교해 보는게 좋을것 같네요

이번 포스팅에서는 MATLAB 의 rounding 관련 함수들에 대해 정리한다. 다음과 같은 벡터 값에 대해 위 함수들을 적용해 보고 그 특징을 알아본다.

x=[-4.3 -1.8 0.7 3.6 1.2+2.6i ]

x =

-4.3000 -1.8000 0.7000 3.6000 1.2000 + 2.6000i

 

  • ceil() : + 무한대 쪽으로 rounding 한다.

ceil(x)

ans =

-4.0000 -1.0000 1.0000 4.0000 2.0000 + 3.0000i

 

  • floor(): ceil 과는 반대로 - 무한대 쪽으로 rounding 한다.

floor(x)

ans =

-5.0000 -2.0000 0 3.0000 1.0000 + 2.0000i

  • fix(): 0 방향으로 rounding 한다.

fix(x)

ans =

-4.0000 -1.0000 0 3.0000 1.0000 + 2.0000i

 

  • round(): 가장 가까운 정수 방향으로 rounding 한다.

round(x)

ans =

-4.0000 -2.0000 1.0000 4.0000 1.0000 + 3.0000i

 

입력이 복소수의 경우에는 real, imaginay 값 각각에 대해 각 함수가 적용 된다.

 

위에서 설명한 rounding 함수들은 몫과 나머지를 구하는데 자주 이용되곤 한다. 예를 들면 다음과 같은 경우이다.

126 / 10 의 몫과 나머지를 구하여라.

몫 : floor(126/10) = 12

나머지 : 126-floor(126/10)*10 = 6

 

나머지와 관련된 함수로는 rem() , mod() 같은 함수가 있다.

rem(), mod() 를 이용하여 나머지를 구해 보자.

 

rem(126,10)

ans =

6

 

mod(126,10)

ans =

6

 

위 결과와 같이 값이 양의 값일 때는 둘 다 같은 결과 6이 나오는 것을 알 수 있다. 하지만 아래 예와 같이 음의 값이 입력되는 경우에는 결과가 다르다.

 

rem(-126,10)

ans =

-6

 

mod(-126,10)

ans =

4

 

rem(x,y) 함수는 x - n.*y 값을 리턴 해 주며 이때 n 을 다음과 같이 처리 한다. n = fix(x./y)

mod(x,y) 함수는 x - n.*y 값을 리턴 해 주며 이때 n 을 다음과 같이 처리 한다. n = floor(x./y)

 

몫을 어떻게 처리 하는가에 따라 결과가 다르므로 함수에 대해 정확히 공부를 하고 사용해야 한다.


  1. a 2015.10.23 13:27

    fix round가 바뀐것 같습니다

  2. 박남욱 2016.06.20 01:35

    fix(-1.8) 은 -1 입니다. 틀리셨어요. 수정 바랍니다.

    • 남성 2016.06.20 08:52 신고

      fix() round() 실행한 내용의 위치가 바뀌었군요. 수정 했습니다.

실험 데이터가 너무 적을 경우 그 경향성을 알아보기 위하여 실험 데이터의 중간 값들을 알아내야 하는 경우가 있다. 이런 경우 interpolation 을 수행하여 수학적으로 중간의 값들을 추정 하곤 한다. 이러한 수치적인 interpolation 을 해 주는 함수 중 하나가 interp1 함수 이다.

 

다음과 같은 데이터를 가정해 보자.

 

X=0:10

X =

0 1 2 3 4 5 6 7 8 9 10

 

Y=sin(X)

Y =

0 0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 0.9894 0.4121 -0.5440

 

위 코드의 Y 값은 간단하게 0, 1, 2, … 10 값에 대한 sin() 값을 나타낸다. 위 값들을 이용하여 X 값이 0, 0.01, 0.02, 0.03, 과 같이 0~10 사이의 중간의 값들을 interp1() 함수를 통해 확인 할 수 있다.

 

X2=0:0.01:10; % 0~10 사이를 0.01 간격으로 설정

 

Y2=interp1(X, Y, X2);

 

plot(X,Y,'r>', X2,Y2,'b--'), grid on

 

 

위 결과를 보면 interpolation 이란 중간 값들을 얻기 위해 빨간색 삼각형 들을 이은 것이라 할 수 있다. interp1() 함수는 기본적으로 linear interpolation 방식을 이용한다. 빨간 삼각형 사이를 직선으로 쭉쭉 이은 것이니 만큼 뭔가 sin 그래프 처럼 보이지가 않는다.

 

그럼 spline interpolation 방식을 이용하여 약간 곡선으로 interpolation 을 해보자.

 

 

Y_spline=interp1(X, Y, X2, 'spline');

 

plot(X,Y,'r>', X2,Y_spline,'b--'), grid on

 

 

그림을 보면 일단 번째 그림 보다 훨씬 sin 그래프에 가까운 것을 확인 수가 있다.

 

다음으로 extrapolation 수행해 보자. extrapolation 이란 데이터 범위 밖의 값들을 추정 하는 이다. 예에서는 x 값이 10 초과 또는 0 미만의 값들은 extrapolation 통해 수학적으로 추정하게 되는 것이다.

 

 

X4=-2:0.1:12;

 

Y_spline_extra=interp1(X, Y, X4, 'spline', 'extrap');

 

figure, plot(X,Y,'r>', X4,Y_spline_extra,'b--'), grid on

 

 

위 그림에서는 0~10 사이가 아닌 값들도 extrapolation 을 통해 수학적인 추정을 함을 확인 할 수 있다. 하지만 0~10 밖의 값들은 sin 그래프의 형태가 아닌 것을 알 수 있다.


MATLAB 에서 매트릭스의 인덱스에 대한 처리는 정말 기본 중에 기본이라 할 수 있다. 이번 포스팅에서는 매트릭스 인덱스 처리에 대해 설명한다.

 

1부터 10까지의 정수로 구성된 벡터 X를 발생시켜 보자.

 

  • 콜론 연산자 이용

X=1:10

X =

1 2 3 4 5 6 7 8 9 10

  • linspace() 함수 이용

X=linspace(1,10,10)

X =

1 2 3 4 5 6 7 8 9 10

 

  • 무식하게 다 써주기

X=[1 2 3 4 5 6 7 8 9 10]

X =

1 2 3 4 5 6 7 8 9 10

 

위 세가지 방법 모두 동일한 결과를 나타낸다. linspace(1, 10, 10) 함수의 첫 번째 1 은 시작 값 , 두 번째 10은 마지막값, 그리고 세번째 10 은 전체 개수를 의미 한다.

 

그럼 X 라는 벡터에서 인덱스를 이용하여 각 인자 또는 벡터를 만들어 보자. MATLAB 의 인덱스는 C/C++ 와는 다르게 1 부터 시작한다는데 주의 해야 한다.

 

인덱스 10 , 2 , 5 번째 값을 추려 보자 .

X([10 2 5])

ans =

10 2 5

 

위 코드와 같이 X라는 변수의 인덱스에 접근을 할 때는 괄호 ( ) 를 하고 그 안에 인덱스 값들을 써 주면 된다. 한 개의 인덱스에 접근 하는 경우에는 X(5) 이런 식으로 한 값만 써주면 되지만 여러 인덱스에 접근하는 경우에는 [ ] 을 써서 벡터 또는 매트릭스의 형태로 써 줘야 한다.

 

다음으로 X 의 값 중에 5 보다 큰 값들만 뽑아 보자.

 

  • logical index 를 이용하는 방법

logical_index=X>5

logical_index =

0 0 0 0 0 1 1 1 1 1

 

class(logical_index)

ans =

logical

 

X(logical_index)

ans =

6 7 8 9 10

 

  • index 를 이용하는 방법

index=find(X>5)

index =

6 7 8 9 10

 

X(index)

ans =

6 7 8 9 10

 

MATLAB에서는 위와 같이 두 가지 방식으로 인덱스에 접근 할 수 있다. 1 부분과 같이 logical 값을 이용하는 방법도 있다는 것을 알아두기 바란다. logical index 를 이용할 때는 코드가 더 간단해 지며, 보통은 아래와 같이 쓰는 게 일반 적이다.

X(X>5)

 

ans =

6 7 8 9 10

 

이제 매트릭스의 인덱스에 대해 살펴 보자.

 

X=magic(3)

 

X =

8 1 6

3 5 7

4 9 2

 

위 X 의 3, 2 , 5 번째 값들을 추려 보자.

 

X([3 2 5])

 

ans =

4 3 5

 

매트릭스는 항상 열을 기준으로 인덱싱을 한다는 것을 알아야 한다.

 

따라서 X 의 3, 2, 5 번째 값을 인데싱 할 때는 매트릭스 X 를 다음과 같이 열을 기준으로 나열해 보면 인덱스를 정확히 알 수 있을 것이다.

 

X(:)

 

ans =

8

3

4

1

5

9

6

7

2

 

위 결과를 보면 첫 번째 열 다음에 두 번째 열, 세 번째 열 이렇게 나열하는 것을 확인 할 수 있다.

 

그리고 위 결과에서 3번째 2번째 5번째 값을 뽑으면 4, 3, 5 라는 값이 출력된다.

 

 

이제 행 과 열에 따른 인덱싱을 알아 보자.

 

위의 X 매트릭스에 대하여 [1,3] 행 2 열 의 값들을 추려 보자.

 

X([1 3], 2)

 

ans =

1

9

 

위와 같이 매트릭스의 행과 열 인덱스를 통해 간단하게 인자 값들에 접근 할 수 있다.


오늘은 MATLAB subplot()에 대해 알아본다. subplot()은 하나의 figure 창에 여러 그래프를 표시 할 때 이용한다. 다음 예를 통해 subplot() 에 대해 설명한다.

 

x=1:5;

y=1:5;

y2=(1:5)*2;

 

figure, subplot(3, 2, 1), plot(x,y)

subplot(3, 2, 4), plot(x,y2,'r:.')

 

 

subplot(3, 2, 1) 에서 첫 번째 3은 그래프의 행의 개수를 의미하고, 두 번째 2는 그래프의 열의 개수를 의미한다. 세 번째 숫자 1은 위 그림에서 빨간색 글자로 표시한 1 부분에 그림을 넣겠다는 의미이다.

subplot(3, 2, 4)는 에서 세 번째 숫자 4는 위 그림의 4번 부분에 그래프를 넣겠다는 의미이다.

 

그럼 위 그림을 다른 방식으로 그려 보자.

 

figure, subplot(3, 2, [1 2]), plot(x,y)

subplot(3, 2, [4 6]), plot(x,y2,'r:.')

 

 

subplot(3, 2, [1 2]) 에서 [1 2] 라고 나타냈으므로 figure 창의 1번 2번 칸에 그림을 나타내겠다는 뜻이며 subplot(3, 2, [3 4])은 [3 4] 라고 나타냈으므로 figure 창의 3번 4번 칸에 그림을 나타내겠다는 뜻이다.

 

이상으로 subplot() 에 대한 설명을 마친다.


통신 또는 신호 처리에 있어서 delay 는 항상 존재 하기 마련이다. 필터링이나, 공기중의 매질을 통과하는 동안의 시간 지연RF 소자에 의한 지연 등.. 굉장히 다양한 지연 요소가 있다. 이런 delay 값을 확인 하기 위하여 수신신호와 기준 신호의 correlation을 이용하곤 한다. correlation 은 상관도로서 reference 값과 들어오는 값이 얼마나 잘 매치가 되는지를 알아 보는 척도이다. 이런 correlation 을 수행하는 함수로 xcorr() 함수가 있다. 본 함수는 Signal processing toolbox 내에 포함된 함수이다.

 

x=1:5

x =

1 2 3 4 5

 

y=[0 0 0 1:5]

y =

0 0 0 1 2 3 4 5

 

위의 x와 y 값을 비교 해 보면 y 가 x 에 비해 3 sample 만큼 delay 가 있다는 것을 확인 할 수 있다.

x와 y 에 대해 correlation 을 취해 보면

[V, D] = xcorr(x,y);

V 에는 correlation 에 따른 결과 값이 입력 되고, D 에는 delay 값이 입력된다. 이에 대해 그래프를 그려 보면 다음과 같은 결과를 확인 할 수 있다.

figure, plot(D, V,'r:.'), grid on

 

 

-3 위치 즉 y 가 3 만큼 delay 된 위치에서 correlation 이 최대가 됨을 확인 할 수 있다. Communication toolbox 에 있는 finddelay() 함수를 이용하면 delay 값을 좀더 쉽게 알 수 있다.

 

D=finddelay(x,y)

D =

3


MATLAB 의 연산은 기본적으로 매트릭스 연산을 기본으로 한다. 따라서 매트릭스 또는 벡터의 연산에 매우 편리하도록 구성이 되어 있다. 반면에 for 나 while 같은 루프 문의 경우에는 비교적 그 처리 속도가 느린 단점이 있다. 따라서 for 나 while 같은 루프 문은 가급적이면 안 쓰는 방향으로 코딩을 하는 것이 좋다. for 나 while 같은 루프 문을 매트릭스 또는 벡터 연산으로 처리 하는 것을 벡터화 기법이라 한다.

다음 수식에 대해 M 파일 코딩을 해 보자.

 

 

위 수식에 대해 for 문을 이용해서 코딩을 하고 tic, toc 을 이용하여 처리 시간을 측정 해보자.

 

tic

sum_v=0; % sum 초기값 설정

for n = 1:100

sum_v =sum_v + n^4; % 4 승

end

fprintf('%g \n',sum_v); % 화면 출력

toc

 

위 코드를 실행하면 Command 창에 다음과 같은 결과가 나옴을 확인 할 수 있다.

2.05033e+009

Elapsed time is 0.000682 seconds.

 

이제 위 for 루프 문을 벡터화 하여 처리를 하고 그 처리 시간을 확인 해 보자.

tic

Sum_v=sum((1:100).^4); % 각 인자 값의 4승을 하기 위해 .^ 연산을 이용

fprintf('%g \n',Sum_v); % 화면 출력

toc

 

위 코드의 결과는 다음과 같다.

2.05033e+009

Elapsed time is 0.000318 seconds.

 

for 문을 이용하였을 경우 0.000682 초가 걸린 반면에 벡터화 해서 처리하는 경우 0.000318 초가 걸리는 것을 확인 할 수 있다. 위 예제는 매우 간단한 예였으므로 그 차이가 그리 크지 않지만 매우 큰 루프문을 이용할 때는 그 차이가 훨씬 더 확연하다.

다음 페이지에서 



2011/03/14 - [programming language/MATLAB] - MATLAB 적분 int(), quad()



MATLAB 에서 이용 가능한 적분에 대해 학습을 한 적이 있다. 오늘은 추가적으로 수치 적분 함수 중 사다리꼴 기법으로 적분을 수행하는 trapz() 함수에 대해 설명한다.

 

trapz() 함수의 원형은 다음과 같다.

Z = trapz(X,Y)

  • X 값은 함수의 입력이고 Y 값은 함수의 출력이다.

 

trapz() 함수를 이용하여 다음 수식에 대한 적분을 수행해 보자.

 

일단 X 의 범위를 정한다. 적분 구간이 0~3 이므로

X=0:3

X =

0 1 2 3

 

Y= X.^2 + 2.*X + 1

Y =

1 4 9 16

 

integral_value=trapz(X,Y)

integral_value =

21.5000

 

실제 위 식에 대해 적분을 해 보면 21 이 되야 하는데 21.5 가 되는 것을 확인 할 수 있다. 이는 X 값의 간격 1 이 너무 크기 때문이다. 따라서 X 값의 간격을 0.1로 좁게 설정을 하고 다시 계산을 수행해보자.

 

X=0:0.1:3;

Y= X.^2 + 2.*X + 1;

integral_value=trapz(X,Y)

integral_value =

21.005

 

위 결과도 21은 아니다. 그렇다면 간격을 0.01로 더 좁게 하고 계산을 다시 해보자.

X=0:0.01:3;

Y= X.^2 + 2.*X + 1;

integral_value=trapz(X,Y)

integral_value =

21.00005

 

X 값의 간격이 1 à 0.1 à 0.01 로 좁아 질수록 값이 21이라는 값에 근접해 감을 확인 할 수 있다.


MATLAB 을 이용하여 적분을 수행해 보자. MATLAB 을 이용하여 적분을 하는 방법은 크게 두 가지 정도로 구분 할 수 있을 것 같다. 

첫 번째로는 수치적인 적분 방법이고 두 번째는 Symbolic math toolbox 를 이용한 수학적인 접근 방법이다.

수치적인 접근 방법이라 하면 사다리꼴 방법 과 같이 함수를 매우 작은 조각으로 나눠서 부분부분의 면적의 합을 구하는 방법이 될 것이다. 물론 이런 적분을 구하는 알고리즘을 짜서 적분을 할 수도 있겠지만 고맙게도 matlab 은 이런 함수를 제공 해 준다.

quad() 함수를 이용하면 수치 적분이 가능하다. 본 함수는 수치 적분이므로 부정적분을 해 주지는 못한다. 따라서 다음 수식과 같이 구간이 정해진 식에 대해서만 적분 할 수 있다.

 

 

위 수식에 대한 적분은 다음과 같다.

 

F = @(x) x.^2+2.*x+1

F =

@(x)x.^2+2.*x+1

 

quad(F,0,3)

ans =

21

간단하게 해가 나오는 것을 확인 할 수 있다.

하지만 이런 수치적분은 위에서 밝힌 대로 부정적분이 불가능 하다. 부정적분을 할 때는 symbolic math toolbox 를 이용해야 한다. Symbolic math toolbox 는 Mathematica 나 Maple 과 같이 기호로서 계산을 해주는 툴로서 복잡한 수학 계산시 굉장히 편리하게 이용할 수 있다.

 

다음과 같은 수식에 대해 간단하게 부정적분 그리고 임의의 구간에 대한 적분을 수행 해보자.

 

 

부정적분은 굉장히 간단하다.

symbol 로 사용한 변수인 x 를 symbolic 으로 설정해준다.

 

syms x

int(x^2+2*x+1)

그리고 int() 함수 내에 적분할 수식을 써 준다. int() 함수는 integral 의 약자이다. 위 코드의 실행 결과는 다음과 같은 부정적분 결과가 나온다.

 

ans =

(x*(x^2 + 3*x + 3))/3

 

그럼 임의의 구간 a~b 구간에 대해서도 적분을 해 보자. a, b 도 symbolic 으로 설정을 해 준 후에 int() 함수에서 적분 구간을 a, b 라고 명시 해 준다.

 

syms x a b

int(x^2+2*x+1,a,b)

ans =

b*(b*(b/3 + 1) + 1) - a*(a*(a/3 + 1) + 1)

 

다음과 같이 간단하게 a, b 라는 적분구간에 대해 적분이 실행 되는 것을 확인 할 수 있다.


  1. 조민형 2013.08.09 16:46

    적분 질문 있습니다 ㅠㅠ
    바쁘실텐데 답변 주실 수 있으시려나요..

    매트랩 초보라 왜 이런지 모르겠습니다.

    식이 좀 복잡한데,,

    >> int((64*(cos((B*q*cos(phi))/2)^2 - 1)*(cos((A*q*sin(phi))/2)^2 - 1)*(sin(C*q)/(2*C*q) - 1/2))/(A^2*B^2*C^2*q^6*cos(phi)^2*(cos(phi)^2 - 1)),phi,0,pi/2)
    Warning: Explicit integral could not be found.

    ans =

    int(((cos((A*q*sin(phi))/2)^2 - 1)*(64*cos((B*q*cos(phi))/2)^2 - 64)*(sin(C*q)/(2*C*q) - 1/2))/(A^2*B^2*C^2*q^6*cos(phi)^2*(cos(phi)^2 - 1)), phi = 0..pi/2)
    와 같이 적분을 찾을 수 없다고 나옵니다.

    이런 경우 부정적분 값을 도출할 수 있는 스킬이나 노하우 같은게 있나요??

    답변 미리 감사드립니다 !!

    • 남성 2013.08.09 19:02 신고

      저도 몇가지 방법으로 시도는 해봤는데 마찮가지로 Warning: Explicit integral could not be found. 라는 메시지가 뜨네요. 도움 못되 죄송합니다. ^^;

  2. 서형관 2016.12.02 09:22

    임의의 구간 적분에 대해서 질문이 있습니다.
    매트랩을 인터넷에 올라온 글을 보면서 독학하는중 입니다.
    구간a b c는 a=1000 b=8000 c=9000으로 정의하였고 그 나눈 부분을 적분하는 프로그램을 짜고 있는데

    syms n1 0 a
    n1=int((i*n)/(2*(pi)*a^2),n,0,a);
    syms n2 a b
    n2=int(i/(2*(pi)*n),a,b);
    syms n3 b c
    n3=int((i*(c^2-n^2))/((2*(pi)*n)*(c^2-b^2)),b,c);


    이렇게 만들어 봤는데 쉽지 않네요 ㅠㅠ
    어떤부분이 문제인지 알 수 있으신가요??

    • 남성 2016.12.02 21:03 신고

      다음과 같이 해 보세요.

      a=1000;
      b=8000;
      c=9000;

      syms n
      n1=int((1i*n)/(2*(pi)*a^2),n,0,a)
      double(n1)

      n2=int(1i/(2*(pi)*n),a,b)
      double(n2)

      n3=int((1i*(c^2-n^2))/((2*(pi)*n)*(c^2-b^2)),b,c)
      double(n3)


  3. 크레스포 2017.04.15 18:22

    Syms로 적분한 함수의 그래프를 출력하는 방법은 없나요?

MATLAB 에서 다항식 또는 연속 함수의 해를 구하는 방법에 대해 설명한다. fzero() 은 연속 함수의 한 지점에서의 해를 구하는 함수이다. 함수의 원형은 x = fzero(fun,x0) 이며 fun 이라는 함수에 대하여 x0 근처에서의 해를 찾아 준다. fzero() 함수는 y 값의 부호가 바뀌는 지점을 찾아서 해를 구하는 방식이다. 내부 알고리즘으로는 bisection, secant, inverse quadratic interpolation methods 이렇게 세 개 알고리즘이 조합된 방식을 이용한다고 한다. 해를 찾는 방식이 y 의 부호 변경 지점을 찾는 것이므로 y=x2 과 같이 x축에 접하는 함수에 대한 해를 찾을 때는 사용하지 않는 것이 좋다.

 

Command 창에 다음과 같은 명령을 해 보자.

fzero('x.^2',0) command 창에는 0 이라는 결과가 나올 것이다. 하지만 다음과 같이 fzero('x.^2',1) 이라고 한다면 ans 로서 NaN 이 뜨게 된다. 이런 경우 때문에 개인적인 생각으로는 fzero() 를 이용할 때는 그래프를 한번 그려 봐서 개형을 파악하고 fzero() 함수를 적용하는 것이 좋다.

 

sin(x) 그래프에 대해 x 값1 근처의 해를 구해 보자.

일단 그래프를 한번 그려 본다. 다음 과 같음 명령을 하면 그래프가 그려진다.

 

ezplot('sin(x)'), grid on

 

위 그래프를 보면 y 의 값들이 0 을 통과 하므로 fzero()를 이용해도 괜찮다고 판단한 후 아래 명령을 이용하여 1 근처의 해를 구한다.

 

fzero('sin(x)',1)

ans =

1.5485e-024

 

결과로 나온 1.5485e-024 은 0 이라는 것을 알 수 있다.

 

다음으로 다항식의 해를 구해주는 roots 함수에 대해 알아 보자.

 

 

위 수식에 대한 해를 구해 보자. 함수의 이용 방법은 굉장히 간단하다.

r = roots([1 4 2 -1])

r =

-3.3028

-1.0000

0.3028

 

위에서 구한 다항식의 해를 이용하여 다항식 계수를 역으로 계산해 보자.

 

coeff=poly(r)

coeff =

1.0000 4.0000 2.0000 -1.0000

 

간단하게 다항식의 계수 값들이 역으로 계산되는 것을 확인 할 수 있다.

물론 poly() 라는 함수를 이용하지 않고 conv() 함수를 이용하는 방법도 있다. 

이에 대해서는 다음 포스트를 확인 하기 바란다.

2011/03/01 - [MATLAB] - conv(), filter(), 인수분해 전개

 


  1. 1 2011.04.07 11:00

    좋은정보 ㄳ

오늘은 ezplot() 함수를 이용하여 그래프를 그리는 방법에 대해서 설명하겠습니다. ezplot() 함수는 함수에 대해 기본적으로 x 축의 범위가 -2π ~ 2π 범위에 대해 그래프를 그려 주는 함수 입니다. 물론 함수의 범위는 사용자가 설정할 수도 있습니다.

 

몇 가지 예를 통해 사용 방법을 설명 드리겠습니다. y=sin(x) + x2 + 2x+4 라는 함수에 대해 그래프를 그려 보겠습니다. 위 수식에 대한 그래프를 그리는 방법은 다음과 같습니다.

 

matlab command 창에서 ezplot('sin(x) + x.^2 + 2.*x+4 ') 라고 타이핑 하면 다음과 같은 결과가 나옵니다. 다음 결과는 x 축의 범위가 위에 설명 드린 데로 -2π ~ 2π 범위의 값이 나왔습니다.

 

그럼 조금 다른 형태로 함수를 표현하고 ezplot() 을 이용하여 -10~10 사이의 축에 대해 그래프를 그려 보겠습니다. 위 수식을 다음과 같이 function handle 을 이용하여서 표현할 수도 있습니다. 일반적으로 함수는 함수파일로 따로 저장을 해서 이용하곤 하지만 다음과 같음 function handle을 이용하면 코딩 중간에 바로 함수를 만들어서 쓸 수 있다는 점에서 편리합니다.

 

f = @(x) sin(x)+ x.^2 + 2.*x+4

f =

@(x)sin(x)+x.^2+2.*x+4

 

이제 위 f 라는 핸들을 이용하여 그래프를 그려보면

 

ezplot(f,[-10 10]), grid on

 

위 결과와 같이 굉장히 간단하게 그래프를 그릴 수 있습니다. ezplot 함수를 사용하실 때 항상 vectorize 하시는 게 좋습니다. 그렇지 않음 다음과 같은 warning 이 뜨게 됩니다. vectorize 에 대해서는 다음 페이지를 참조 하시기 바랍니다.     


 2010/01/01 - [MATLAB] - MATLAB - vectorize()

 

Warning: Function failed to evaluate on array inputs; vectorizing the function may

speed up its evaluation and avoid the need to loop over array elements.

 

마지막 예로 원을 하나 그려 보겠습니다.

 

보통 반지름 r 이고 중점이 (a, b) 인 원은 다음과 같이 표현 됩니다.

 

 

중점이 (2, 3) 이고 반지름이 3인 원을 한번 그려 보겠습니다.

 

먼저 function handle 로 함수를 지정 해 줍니다.

 

f = @(x,y)(x-2).^2 + (y-3).^2 -3^2

f =

@(x,y)(x-2).^2+(y-3).^2-3^2

 

다음으로 ezplot() 함수를 이용하여 그림을 그려 줍니다.

 

ezplot(f), grid on

 

다음과 같은 결과가 나오는 것을 확인 할 수 있습니다.

 



  1. 와글 2014.12.16 03:05

    필요한 지식이였는데 도움많이됬습니다 감사합니다!

polyfit() 함수는 입력과 출력 값으로부터 다항식의 계수를 찾아 주는 함수 입니다. 예를 들면 라는 식에서 x, y 값을 알고 있으면 다항식의 계수값 a, b, c 값을 찾아 준다는 것입니다. polyfit() 함수의 원형은 다음과 같습니다.

 

p = polyfit(x,y,n)

 

  • x 는 다항식의 입력값, y 는 다항식의 결과값입니다.
  • n 은 차수를 의미 합니다.
  • p 는 차수에 따른 다항식의 계수 값을 의미 합니다.

 

다음과 같은 x, y 값에 대하여 다항식의 계수를 찾아보겠습니다.

x=0:5

x =

0 1 2 3 4 5

 

y=2*x.^3 + 5*x.^2 + 6*x+4

y =

4 17 52 121 236 409

 

위와 같은 값이 있을 때 다음과 같이 3차 다항식의 계수를 찾습니다.

결과 p 의 값을 보면 위에 있는 y 식에서의 계수를 정확히 찾는 것을 볼 수가 있습니다.

 

p=polyfit(x,y,3)

p =

2.0000 5.0000 6.0000 4.0000

 

만약 4차 다항식으로 fitting을 해보면

 

p4=polyfit(x,y,4)

p4 =

0.0000 2.0000 5.0000 6.0000 4.0000

 

다음과 같은 값이 나오는 군요.

위 결과로부터 알 수 있는 것은 차수를 모를 때는 최대한 큰 차수를 넣어 보는 것이 좋습니다.

위 결과에서는 4 차 다항식의 계수 값은 0 이므로 x, y 사이의 관계는 3차 다항식이라는 것을 알 수가 있습니다.

 

다음으로 다항식 계수로부터 값을 계산하는 함수인 polyval 함수에 대해 알아봅니다.

polyfit 함수로 찾은 계수의 값이 제대로 나온 건지 검증 하기 위하여 polyval 함수를 주로 이용합니다.

 

함수의 원형은 아래와 같습니다.

 

y = polyval(p,x)

p 는 다항식 계수 값이고 x 는 입력 값입니다.

 

앞에서 y=2*x.^3 + 5*x.^2 + 6*x+4 라고 표현했던 부분은

polyval 함수를 이용해서 y=polyval([2 5 6 4], x) 라고 표현해도 동일한 결과가 나오는 겁니다.

 


MATLAB 으로 M 파일 작성시 가끔 다른 사람들과 파일을 교환해야 하는 경우가 종종 있습니다. 이런 경우 회사 또는 개인의 보안상 파일의 소스는 숨겨서 보내야 할 경우가 있죠. C/C++ 같은 경우에는 라이브러리 파일처럼, MATLAB 에는 .p 라는 확장자를 갖는 파일을 만들 수가 있습니다.

 

이렇게 m 파일을 p 파일로 바꿔주는 함수가 바로 pcode 함수 입니다.

 

pcode 함수의 사용 방법은 굉장히 간단 합니다.

 

pcode abc.m 과 같은 형태로 써주면 abc.p 라는 pcode 가 만들어 집니다. 외부 또는 다른 사람에게 이 abc.p 파일을 넘겨주면 파일 소스는 볼 수 없지만 실행은 시킬 수가 있습니다.

 

pcode *.m 이라고 써 주게 되면 현재 폴더 안의 모든 m 파일들을 다 p 파일로 바꿔지기도 합니다.

 

같은 이름을 가진 abc.p 파일과 abc.m 파일이 같은 폴더 안에 있을 경우에는 주의를 해서 이용해야 합니다.

abc 라는 파일을 실행 시킬 때 같은 이름의 m 파일이 있어도 p 파일의 실행이 우선 합니다.

  1. 바나나 2018.02.26 23:10

    관리자의 승인을 기다리고 있는 댓글입니다

오늘은 EYE PATTERN 에 대해 알아 보려 합니다. 


학부 때 EYE PATTERN 이란 걸 보고 이게 뭘까 ~ 정말 도무지 감이 안 왔던 기억이 있습니다.


EYE PATTERN 이란 특정 시간 구간 동안 들어오는 신호의 파형을 계속 겹쳐서 나타낸 것뿐입니다.


보통 EYE PATTERN을 그릴 때는 2 심볼 구간 동안 표시를 하곤 하죠.


오늘 EYE PATTERN을 하기 전에 BPSK 에 모르시는 분들은 아래 글을 읽고 오시기 바랍니다.

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

 

그리고 Raised cosine filter 에 대해서는 아래 글을 참조해 주세요.

2011/03/02 - [MATLAB] - [디지털 통신] Raised Cosine Filter

 


그럼 이제 matlab 코딩을 시작 해 보죠. BPSK 신호의 송수신을 예로 들어서 설명해 보겠습니다. 잡음은 고려하지 않겠습니다.


  • 먼저 Square root raised filter 를 만듭니다. 8 oversamplegroup delay 3심볼, roll off 0.2 인 필터를 만들겠습니다.

Sqrt_r_filter= rcosine(1e6, 8e6, 'sqrt', 0.2, 3);



  • 비트를 만듭니다. 여기선 5000 비트를 만들겠습니다.


bits=randi([0 1], 5000, 1);



  • 비트 0은 -1로 비트 1은 1로 Mapping 해줍니다.


Mapped_signal=bits*2-1;



  • 8 oversample 이므로 8배upsample 을 합니다. upsample 이란 중간에 0을 채워 넣는 겁니다.

Upsampled_signal=upsample(Mapped_signal, 8);


  • Tx filtering 을 수행합니다.

Tx_signal=filter(Sqrt_r_filter,1, Upsampled_signal);



  • 수신 단에서는 Rx filtering 을 수행합니다.


Rx_signal=filter(Sqrt_r_filter,1, Tx_signal);

  

이제 신호의 송수신은 다 됐습니다.

근데 위에서도 말했듯이 EYE PATTERN 을 보기 위해서는 특정 구간으로 겹쳐서 그려야 합니다.

따라서 위의 수신 신호 Rx_signal 을 적절히 구성만 해주면 됩니다.


2 심볼 구간 동안의 EYE PATTERN을 보기 위한 2 심볼은 8*2 = 16 샘플 입니다.

그리고 Rx_signal 의 1부터 48번째까지의 값은 delay 값 일명 쓰레기 값입니다. 이는 필터의 delay 가 3 심볼이니깐 8*3에다 Tx , Rx 필터 통과니깐 또 *2 를 해서 48 이라는 숫자가 나오는 겁니다.


이제 다음과 같이 매트릭스를 구성합니다.

 

Rx_signal1=Rx_signal(49);

Rx_signal_cutted=Rx_signal(50:end);

N_column=floor(length(Rx_signal_cutted)/16);

Rx_signalmatrix=reshape(Rx_signal_cutted(1:16*N_column),16,N_column);

 


아래 코드의 빨간색 부분은 이전 신호의 온 타임 지점을 반복시켜 주기 위한 것입니다.


Rx_signalmatrix_all=[[Rx_signal1 Rx_signalmatrix(end,1:end-1)] ;Rx_signalmatrix];

 

이제 그림을 그려줍니다.


figure, plot(Rx_signalmatrix_all,'b'), grid on

 

결과는 다음과 같습니다.

 

그럼 ~ 맨~ 처음 만들었던 필터에서 roll off 를 바꿔가며 위에서 진행한 실험을 똑같이 해 봅니다.


roll off 가 0.5 일 때

 

roll off 가 0.7 일 때

 

roll off 가 1 일 때

 

쭉~ 보니깐 roll off 가 올라갈수록 그림이 깔끔해 지죠? roll off 가 올라갈수록 EYE 는 점점 더 열리게 되 므로 신호 품질은 더 좋아지는 겁니다.


하지만 공학이란 게 어디까지나 트레이드 오프가 있는 거니깐~ roll off가 높아진다고 무조건 좋은 건 아닙니다. roll off 가 높아질수록 BW(대역폭) 을 더 이용하게 되는 겁니다.

   

EYE PATTERN 이란 게 직접 시뮬레이션 해 보면 별거 아닌데 처음에 보면 이해가 어렵습니다. 이 글을 보시는 분들도 위 코드대로 직접 시뮬레이션 한번 해 보시기 바랍니다.


  1. 김도윤 2012.05.02 12:38

    여쭤볼게 있는데요.
    Upsample을 하면 bit간 간격이 넓어지니까 filter도 그에 따라서 time domain에서 늘려줘야 하지 않나요?
    그리고 마지막에 신호의 온타임 지점을 반복한다는 것이 무슨 말인지 잘 이해 못하겠네요.
    혹시 댓글 보시면 답변해주시면 감사하겠습니다.

    • 남성 2012.05.03 00:03 신고

      upsample 이란게 0 을 끼워 넣어주는것이므로 샘플 간의 간격은 좁아지는 겁니다. 즉 sample rate가 그만큼 증가하는 것이죠.

      위 포스팅에서 아래 코드는 필터 계수를 설계하기 위한 것인데~~

      Sqrt_r_filter= rcosine(1e6, 8e6, 'sqrt', 0.2, 3);

      1e6, 8e6 이므로 8배 upsample 에 대해 filter 계수를 생성하는 겁니다.

      예를 들어 4배 upsample 에 대해 filter 계수를 생성하고 싶다면

      rcosine(1, 4, 'sqrt', 0.2, 3); 하시면 됩니다.


      위 포스팅에서 아이패턴은 2심볼 구간 동안그린건데 17 샘플 마다 그려 주게 됩니다.

      즉 온타임 지점은 8 샘플 마다 오니까 예를 들자면

      1 9 17 를 그린 다음에

      17 25 33 를 그리고

      33 41 49 를 그리고

      49 57 65 를 그립니다.

      위 숫자들을 보면 17 33 49 번째 샘플을 한번씩 더 그린거죠. 그렇게 안 그리면 아이패턴을 그릴때마타 한 샘플씩 어긋나게 그리게 되겠죠.



  2. scac 2012.09.20 17:39

    roll-off가 올라가서 eye pattern이 열리는게 아니라, square root raised cosine의 양 옆 진동하는 부분을 너무 짧게 잡아서 저렇게 눈이 닫히는거 아닌가요? 콘볼루션 연산 하는데 저렇게 짧으면 값이 제대로 안나올텐데요;
    포스팅에는 3개 밖에 없는데 이것을 300개로 늘리고 roll-off는 낮은 값으로 해보십시오.

    • 남성 2012.09.20 17:59 신고

      방문해 주셔서 감사합니다.
      위 포스팅의 경우 raised cosine filter 의 group delay 3 을 줬고 oversample rate 가 8 이므로 8*3*2+1 = 49 tap raised cosine filter 를 사용하게 되거든요. 300 이라는 건 group delay 를 말하는 건가요? 보통 group delay 를 그렇게 크게 주진 않는데 어떤걸 300 을 말씀하시는 건지요? group delay 300 이 되면 filter tap 이 4801 tap 이 되는데... raised cosine filter 로 이렇게 긴 필터를 사용하진 않죠.

  3. scac 2012.09.20 18:11

    비트수가 5000개인데 square root raised cosine 의 양 옆 진동하는 부분(그룹 딜레이라고 하는가요?)이 만약 3이면 그 진동하는 부분이 좀 떨어져 있는 다른 비트에 닿지 않습니다. (간섭이 안 일어난다 하더라도 닿게는 해야한다고 생각하는데;)
    그리고 제가 말한 300이라는 숫자는 그냥 막연히 적은 숫자이구요. 저는 rcosfir를 사용해서 좀 헷갈리네요;;

    교과서의 eye pattern은 roll-off를 낮춘다고 해서 간섭이 일어나지 않고 한 점에서 만나거든요.

    • 남성 2012.09.20 21:42 신고

      group delay 3 이란 이전 3 심볼 구간 동안에 영향을 미친다는 겁니다. 따라서 보통 첫번째 온타임 지점은 3 심볼 만큼의 딜레이 후에 잡게 되죠. 각 심볼의 sidelobe 부분이 이전 3개와 심볼과 나중 3개 심볼에 영향을 주는거죠.

      보통 raised cosine filter 는 송수신단 square root raised cosine 로 나누어져서 matched filter 의 형태로 사용되죠.
      말씀하신 rcosfir() 함수는 raised cosine filter 입니다.
      위에 제가 설명한 것은 square root raised cosine 필터이고 raised cosine filter 의 경우에는 roll off 에 상관 없이 점으로 나오는게 맞습니다.
      위에 함수중 rcosine() 함수의 'sqrt' 부분을 'normal' 로 바꾸면 raised cosine filter 가 됩니다.

  4. scac 2012.09.20 21:57

    본문 포스팅에 roll-off 인자가 0.2인데 수신신호 샘플링 타임 지점이 너무 두껍게 나와서 말씀드린 것입니다. 제가 알기로는 roll-off 인자가 0~1에서 어떤 값을 가지든 수신필터, 송신필터가 각각 square root raised cosine일 경우 한 점에서 만나는 것으로 알고 있어서요. 본문 포스팅 이미지는 너무 두껍지 않나 해서..
    그리고 rcosfir도 'sqrt'를 넣어서 사용하는걸 말씀 안드렸네요.. :)


    이해를 돕기 위해 본문 포스팅에 있는 매트랩에서 숫자만 수정해 보았습니다. 말씀하신 그룹딜레이 3을 100으로 늘렸습니다. 본문포스팅 roll-off 가 0.2인 결과물보다 훨씬 나아보입니다.

    Sqrt_r_filter=rcosine(1e6, 8e6, 'sqrt', 0.2, 100);

    bits=randi([0 1], 5000, 1);

    Mapped_signal=bits*2-1;

    Upsampled_signal=upsample(Mapped_signal,8);

    Tx_signal=filter(Sqrt_r_filter,1,Upsampled_signal);
    Rx_signal=filter(Sqrt_r_filter,1,Tx_signal);



    Rx_signal1=Rx_signal((2*8*100)+1);
    Rx_signal_cutted=Rx_signal((2*8*100)+1+1:end);
    N_column=floor(length(Rx_signal_cutted)/16);
    Rx_signalmatrix=reshape(Rx_signal_cutted(1:16*N_column),16,N_column);


    Rx_signalmatrix_all=[[Rx_signal1 Rx_signalmatrix(end,1:end-1)] ; Rx_signalmatrix];

    figure, plot(Rx_signalmatrix_all,'b'), grid on

    • 남성 2012.09.21 09:42 신고

      네 scac 님이 하신 말씀이 맞습니다.
      교과서에서 말하는 raised cosine filter 는 이론적으로 무한대의 탭을 가지는 필터이므로 당연히 딱 점이 되서 나옵니다.

      제가 얘기하는것은 Sqrt_r_filter=rcosine(1e6, 8e6, 'sqrt', 0.2, 100);
      를 하면 Sqrt_r_filter 의 경우 1601 탭이 나오는데. 이것도 물론 이상적인 것은 아니지만, 구현 불가능 할 정도로 긴 필터 탭이라서 실제로 써 먹기는 어렵다는 거구요.

      즉 위에 말씀 하신 필터 같은 경우엔 이론적으로는 가능한데 실제로 써먹기는 거의 불가능 할 정도로 너무 길다는 겁니다.

      실제로 FPGA 등에서 멀티플라이어 수는 제한이 있으니까요.

      따라서 1601 탭 같이 긴 필터 탭으로 만들면.. 실제로 구현이 안되는, 알고리즘을 위한 알고리즘을 하는 꼴이 되는거죠.

  5. 뽀로리 2013.09.15 13:58

    지금 디지털통신공학및설계 수업때문에 15개의 그래프를 위해서 시뮬레이션 돌리기를 시도중입니다.
    처음 만저보는 프로그램이기도 하면서 통신도 잘 몰라 매트랩 하는데 어려움이 많은데요.
    혹시 4진 eye pattern은 어떻게 만들 수 있을까요?
    책은 Introduction to Analog & digital communications 2nd Simon Haykin/Michale Moher 입니다.
    메일 kkb8594@gmail.com입니다.
    제발 조금 도와주십시오. 부탁드립니다.

    • 남성 2013.09.15 17:29 신고

      4진 아이패턴 이라는게 16 QAM 에대한 아이패턴을 말씀하시는건가요?
      위 포스팅은 BPSK 에 대한 예인데.. 16 QAM 이라면 http://iamaman.tistory.com/205 글 보고 mapping 만 바뀌면 됩니다. 일단 변조 방식에 대한 이해가 필요 하실듯 하네요

MATLAB 은 대화형 언어로서 Command 창에 명령어를 치면 결과가 바로 밑에 뚝딱 뚝딱 나옵니다.

그런데 코드가 길어 진다면 Command 창에서 작업을 하는 것은 비 효율적입니다.

그래서 보통은 Editor 창에서 코딩을 하고 M 파일로 저장하고 실행을 시켜서 Command 창에서 확인을 합니다.

Editor 창에서 F5 를 누르거나 초록색 삼각형 모양으로 생긴 실행 단추를 클릭하면 M 파일이 실행이 됩니다.

 

오늘 소개할 내용은 MATLAB Editor 창의 파일 비교 기능 입니다.

M 파일을 버전 별로 작성하다 보면 꼭 변경 내용을 추적해야 될 때가 있습니다.

코드가 짧다면 상호 비교하는데 문제가 안되지만, 코드가 100 줄 넘어가기 시작하면 비교하는 것도 참 눈 아프고 힘든 일입니다.

이럴 때 오늘 설명하는 파일 비교 기능을 이용하면 편리 합니다.

오늘 설명 드리는 Comparison 툴은 파일 뿐만 아니라 폴더 및 하위 폴더 그리고 mat 파일 까지도 비교가 가능합니다.

 

Editor 에서 다음과 같이 두 파일을 열었습니다.

 

 

위 그림의 예에서 왼쪽 파일은 plot_file2.m 파일이고 오른쪽은 plot_file1.m 파일입니다.

양쪽이 다른 점은 빨간색 네모 칸 친 부분이 다릅니다.

plot_file1.m 파일에 커서가 있는 상태에서 다음과 같이 Tools 메뉴를 클릭하고 그 하위에 Compare against 그리고 그 하위에 Choose 를 클릭합니다.

그럼 아래와 같은 창이 나오는데, 그림처럼 파일 이름 등을 지정 하고 text comparison 을 설정하고 Compare 버튼을 누릅니다.

 

 

그럼 또 이런 창이 나옵니다.

 

 

위 창에서 보면 일단 맨 밑의 빨간색 네모 칸을 보면 10줄은 동일 하고 왼쪽에서 2줄이 틀리고 오른쪽에서는 1줄이 틀린다고 나옵니다.

그리고 가운데에 빨간색 네모 칸 친 부분을 보면 < 기호가 있고 x 기호가 있는 게 보입니다.

<, > 기호는 한쪽에만 내용이 있을 때 나타납니다. 그리고 x 기호는 둘 다 있긴 있는데 일부가 다를 때 나타나죠.

 

맨 위쪽에 있는 $ 기호처럼 생긴 부분을 클릭하면 빈 줄은 비교를 하지 않습니다.

틀린 부분들에 대한 이동은 맨 위에 있는 빨간색 네모 부분의 화살표 들을 클릭하면 됩니다.

위 예는 설명을 위해 굉장히 간단하게 만들어서 위 compare 기능이 그리 좋지 않은 것처럼 보이지만 한 100 줄만 넘어가는 코드를 비교할 때는 정말 빠르게 다른 부분을 찾을 수가 있습니다.



MATLAB 에는 다항식 계수부분 분수 사이의 변환을 해주는 함수로서 residue() 함수를 제공해준다.

 

residue() 함수의 원형은 다음과 같다.

 

[r,p,k] = residue(b,a)

 

[b,a] = residue(r,p,k)

 

위 식에서

  • b 는 다항식의 분자 부 계수 이다.
  • a 는 다항식의 분모 부 계수 이다.

 

  • r 은 나머지 성분이다.
  • p 는 pole 성분이다.
  • k 는 몫의 값이다.

 

 

이라는 다항식에 대해 부분 분수 전개를 수행해 보자.

 

수학적으로는 다음과 같다.

 

                                                  (1)

 

이제 residue() 함수를 이용하여 부분 분수 전개를 해보자.

 

  • 분자 부의 계수를 설정한다.

b=[2 1]

b =

2 1

 

  • 분모 부의 계수를 설정한다.

a=[1 3 2 ]

a =

1 3 2

 

  • 부분 분수 전개를 수행한다.

[r,p,k]=residue(b,a)

r =

3

-1

p =

-2

-1

k =

[]

 

결과를 보면 나머지(r) 이 3, -1 임을 알 수 있고 pole(p) 값은 -2, -1 그리고 몫(k)은 없음을 확인 할 수 있다.

 

이를 우리가 수학적으로 계산했던 식 (1) 번의 결과와 동일함을 알 수 있다.

 

  • 다음으로 부분 분수로부터 다항식 값을 얻어 보자.

[b,a] = residue(r,p,k)

b =

2 1

a =

1 3 2

 

위 결과와 같이 우리가 초기에 설정 했던 초기 다항식 계수 들이 그대로 나오는 것을 확인 할 수 있다.



MATLAB 에서 그래프를 그려 보자.

x=1:10

x =

1 2 3 4 5 6 7 8 9 10

 

y=1:10

y =

1 2 3 4 5 6 7 8 9 10

 

figure, plot(x, y, 'r*--'), grid on

 

명령어를 치면 아래와 같이 윈도우 창의 중간 부분 그림이 뜬다.

 

 

 

 

저런 그림을 여러 그릴 경우에는 그림이 겹쳐서 한번에 보기에는 곤란한 경우가 많다.

 

이런 경우에 그림의 위치를 설정해 줄수 있다면 굉장히 편리하다.

 

matlab 에서 화면의 크기를 알려면 다음과 같이 명령어를 입력하면 된다.

 

get(0, 'screensize')

ans =

1 1 1680 1050

 

 

모니터는 현재 1680 바이 1050 픽셀 크기를 가지고 있다고 나온다.

 

그럼 이제 화면상의 하단 상단 두개의 그래프를 넣어 보자.

 

figure('units', 'normalized', 'pos',[0.0 0.0 0.5 0.5])

plot(x,y,'r*--'), grid on

 

 

figure('units', 'normalized', 'pos',[0.5 0.5 0.5 0.5])

plot(x,y,'b>--'), grid on

 

 

명렁과 같이 했을 다음과 같은 화면이 보이는 것을 확인 할 수 있다.

 

 

 

그럼 코드에 대해 설명해 본다.

 

figure('units', 'normalized', 'pos',[0.0 0.0 0.5 0.5]) 에서

 

  • 'units', 'normalized' 부분은 화면의 단위를 1 정규화 하겠다는 뜻이다. 픽셀의 개념으로 쓰는 것이 아니라 가로든 세로든 최대 값은 1 정규화 해서 쓰겠다는 것이다.
  • 'pos',[0.0 0.0 0.5 0.5] 그래프가 놓인 위치를 지정하는 것이다.
    • pos position property 표현한 것이다. position 이라고 줘도 되고 구분만 된다면 예와 같이 pos 정도까지만 줘도 position property 라는 것을 알아 듣는다.
    • [0.0 0.0 0.5 0.5] 부분은 [left, bottom, width, height] 나타낸다. 모니터 화면의 왼쪽 아래 모서리가 0, 0 지점이다.

 

 

이상으로 matlab 에서의 그래프 위치 지정에 대한 설명을 마친다.

 

그래프의 위치를 자유 자재로 지정 함으로써 한번에 여러 그래프를 봐야 매우 편리 하게 이용할 있다.

 


  1. 넬핀 2012.01.04 10:05

    좋은 정보 잘 봤습니다. 단순히 figure는 숫자로 구분하는 기능만 있을 줄 알았는데 위치까지 조절해서 깔끔하게 실행할 수 있겠네요 +_+
    감사합니다~!

    • 남성 2012.01.04 10:48 신고

      도움이 되었다니 다행이네요~ 방문해 주셔서 감사합니다. ^^

  2. 백곰사냥꾼 2013.04.06 10:55

    왕 편한 자료입니다 ㅎㅎ
    그래프 띄울때 항상 구석에 가있어서 불편했는데 좋은자료 ㄳ

오늘은 MATLAB 기본 함수들에 대해 설명하려 한다.

너무 많은 함수들이 있어서 다 설명하기는 힘들고 생각날 때 마다 본 페이지에 추가 하는 방식으로 설명하는 것이 좋을 것 같다.

 

아래와 같은 매트릭스를 예를 들어 보자.

x=magic(4) % 4행 4열의 마 방진 매트릭스 행 또는 열의 모든 합이 같은 매트릭스를 만들어 준다.

x =

16 2 3 13

5 11 10 8

9 7 6 12

4 14 15 1

 

flipud() 함수는 위 아래를 변경 할 때 이용한다. 즉 행의 순서를 바꿀 때 사용한다.

x2=flipud(x)

x2 =

4 14 15 1

9 7 6 12

5 11 10 8

16 2 3 13

 

fliplr() 함수는 좌 우를 변경 할 때 이용한다. 즉 열의 순서를 바꿀 때 사용한다.

x3=fliplr(x)

x3 =

13 3 2 16

8 10 11 5

12 6 7 9

1 15 14 4

 

4행 4열의 매트릭스 x 를 2행 8열로 만들어 보자. 다음과 같이 reshape() 함수를 이용하면 간단히 매트릭스의 행 열을 변환 할 수 있다.

x4= reshape(x, 2, 8)

x4 =

16 9 2 7 3 6 13 12

5 4 11 14 10 15 8 1

 

x4 매트릭스를 하나의 벡터로 표현해 보자.

x4(:)

ans =

16

5

9

4

2

11

7

14

3

10

6

15

13

8

12

1

MATLAB 의 기본 계산은 항상 열을 기준으로 표현된다는 것을 명심해야 한다.

 

반복되는 매트릭스를 만들어 보자.

예를 들면 x 라는 매트릭스를 2행 3열 만큼 만들어 보자.

x5=repmat(x, 2, 3)

x5 =

16 2 3 13 16 2 3 13 16 2 3 13

5 11 10 8 5 11 10 8 5 11 10 8

9 7 6 12 9 7 6 12 9 7 6 12

4 14 15 1 4 14 15 1 4 14 15 1

16 2 3 13 16 2 3 13 16 2 3 13

5 11 10 8 5 11 10 8 5 11 10 8

9 7 6 12 9 7 6 12 9 7 6 12

4 14 15 1 4 14 15 1 4 14 15 1

 

위의 x5 매트릭스의 총 사이즈를 알아 보자.

 

size(x5)

ans =

8 12

 

4행 4열의 매트릭스 x 를 2행 3열로 만들었으므로 위 결과와 같이 8행 12열이 나옴을 알 수 있다.

다음으로 x5의 length()를 확인 해 보자.

length(x5)

ans =

12

 

length() 함수는 size() 의 값 중 큰 값을 리턴 해 준다. size() 값이 8 행 12열 이니까 그 중 큰 값인 12를 반환 주는 것이다.

 

이제 x5 매트릭스의 총 인자의 수를 알아 보자.

numel(x5)

ans =

96

 

size(x5) 값이 8행 12열 이므로 총 인자의 수는 96 이라는 것을 알 수 있다.


Raised cosine filter 는 펄스 쉐이핑 뿐만이 아니라 대역 제한을 위한 필터로서 통신의 가장 기본적인 필터라 할 수 있다. 송 수신단에 각각 square root raised cosine 필터를 적용함으로써 Matched filter 로서 이용하는 게 일반적이다.

오늘은 MATLAB 을 이용하여 raised cosine filter 를 설계하는 방법에 대해 알아 본다.

매트랩의 Communications Toolbox 가 있는 경우 rcosine() 이라는 raised cosine filter 설계 함수를 제공해 준다.

함수의 기본형은 다음과 같다.

[num,den] = rcosine(Fd,Fs,type_flag,r,delay,tol)


위 기본형에서 각 인자에 대해 설명해 본다.

  • num 은 필터의 numerator 값이다. FIR 형태로 필터 설계시에는 num 값만 나오게 된다.
  • den 은 필터의 denomerator 값이다. IIR 형태의 필터 설계시에 값이 나오게 된다.
  • Fd 는 입력 샘플의 샘플링 주파수 이다.
  • Fs 는 필터의 샘플링 주파수 이며, Fs/Fd 값은 반드시 양의 정수가 되어야 한다.
  • type_flag 로는 다음과 같이 4개의 옵션이 있다.
  1. 'default' or 'fir/normal' : 기본적으로 적지 않게 되면 FIR 필터이며 Raised cosine filter 를 설계하게 된다.
  2. 'iir' or 'iir/normal' : IIR 형태의 Raised cosine filter 를 설계하게 된다.
  3. 'sqrt' or 'fir/sqrt' : FIR 형태의 square root raised cosine filter 를 설계한다.
  4. 'iir/sqrt' : IIR 형태의 square root raised cosine filter 를 설계한다.
  • r 은 필터의 roll off 값이다.
  • delay 는 필터 지연 값을 설정 할 수 있다. 이 값은 입력 신호를 기준으로 한다.
  • tol : IIR 필터 설계시의 tolrance 값이며 FIR 형태의 필터 설계시에는 사용하지 않는다.

 

그럼 예제로서 입력 신호의 샘플 주파수가 1 MHz 이고 필터의 샘플링 주파수가 8 MHz, roll off 값이 0.5, 필터의 delay 가 3 심볼인 FIR 형태의 square root raised filter 를 설계 해 보자.

간단하게 다음 과 같이 된다.

num = rcosine(1e6, 8e6, 'sqrt', 0.5, 3);

 

위 필터 계수에 대해 impulse response 를 그려 보면 다음과 같다.

 

주파수 응답은 MATLAB 커멘드 창에 freqz(num,1) 라고 타이핑 하면 다음과 같이 볼 수 있다.




  1. 페인터 2011.12.08 02:51

    좋은글 좋은자료 잘읽고 갑니다. 큰도움이 되었습니다(_ _)

conv() 함수와 filter() 함수의 차이를 알아보고자 한다.

다음과 같이 두 행벡터에 대하여

  

x =[1 2 3]

 

x =

1 2 3

 

y=[3 4 5]

 

y =

3 4 5

 

convolution 해보면 결과는 아래 같다.

C=conv(x,y)                                                 (1)

 

C =

3 10 22 22 15

 

 

그리고 x, y 대하여 filtering 해보면 결과는 아래와 같다.

 

F=filter(x, 1, y)                                             (2)

 

F =

3 10 22

 

봐도 conv() 함수를 적용한게 길게 나온 것을 확인 수가 있다.

conv() 함수는 벡터의 길이를 더한 -1 만큼 길이의 결과가 나온다.

예에서는 x 길이가 3 이고 y 길이가 3 이므로 3+3-1 = 5 만큼의 결과가 나오는 것이다.

 

filter() 함수는 입력값이 길이만큼 출력이 나온다. filter() 함수의 입력 y 길이가 3 이므로 F 결과는 3 나오는 것이다.

디지털 필터에는 FIR IIR 형태의 필터가 있으며 여기서는 numerator 성분 (2) flter() 함수의 두번째 인자가 1 이므로 FIR 필터링을 수행한 것이다.

 

filter() 함수를 이용하여 conv() 결과와 동일하게 나오게 하기 위해서는 입력값 y 에다 0 추가해 주면 된다.

 

F2=filter(x,1,[y 0 0])                                            (3)

 

F2 =

3 10 22 22 15

 

(3) 같이 입력값 y 뒤에 0 두개 추가 주면 (1) 결과와 동일한 것을 확인 있다.

 이제 conv() 함수와 인수분해 전개의 유사점에 대해 알아 보자.

다음과 같은 식이 있을 경우

                                                                                          (4)

 

식은 다음과 같이 전개가 된다.

 

                                                                                                                 (5)

 

MATLAB 에서 이런 전개를 경우에 conv() 이용할 있다.

(4) 괄호 안에 있는 계수 값들을 차수의 순서 대로 묶어서 conv() 취하면 아래와 같은 결과가 나온다.

    y=conv([2 3], [3 2 1])

 

y =

6 13 8 3                                         (6)

 

결과적으로 (6) 결과를 (5) 비교해 보면 계수 값이 같은 것을 있다.

 

  

Matlab Editor 에서 스크립트 코딩시에 cell 모드 이용 방법에 대해 알아 본다.

 

Cell mode 는 코드를 구분하여 실행시킬 때 쓸 수 있고, matlab 코드로부터 report 작성시에도 각 단락 등을 구분 시켜 주는 역할을 한다.

아래 보는 밑줄친 부분(%% + 기호)을 클릭하면 그 부분에 %% 기호가 나타나면서 가로줄이 그어 진다.

또는 단순히 %% 를 써도 되고 control+space 단축키나 cell 메뉴에 들어가서도 셀 설정을 할 수가 있다.

  

셀 모드에서 %% 원본 이라 되어 있는 셀에 커서를 놓은 후에

동그라미 친 부분(Evaluate cell) 을 누르면 해당 셀만 실행이 되고 위 코드의 경우엔 아래와 같은 그래프가 나온다.

 

그 아래 셀로 커서를 옮겨서 14 번째 줄에서 2*x +1 에서 숫자 2를 2 만큼씩 크게 해 가면서 그래프를 그려 보자

물론 이 경우에 2를 4로 바꾼 후 에 현재 셀을 실행 시키고 또 6으로 바꾸고 현재셀을 실행 시키는 형태로도 할 수가 있지만, 여기선 어디까지나 셀 모드의 기능을 익히기 위한 것이므로

위 첫번재 그림의 빨간색 네모 칸에 2를 쓴 후에 그 옆의 + 버튼을 눌러 보자.

+ 버튼을 두번 눌렀을 때 14 번째 줄은

y=6*x+1;

로 바뀌며

위의 그림 역시


이렇게 바뀌는 것을 알 수 있다.

이처럼 cell mode 의 기능들을 이용하면 단순 반복 해야 하는 일들을 보다 효율적으로 수행할 수 있다.



오늘 소개할 MATLAB의 기능은 matlab shortcuts 이다.

보통 MATLAB 유저들이 많이 간과하고 넘어가는 부분인데 이 기능을 잘 이용하면 매우 편하게 current directory를 이동할 수가 있다.

 

보통의 경우 메인 창의 Current folder 창에 작업하고자 하는 폴더를 써 넣거나

Command window 에서


  >> cd c:   


이런 식으로 폴더를 변경하곤 한다.

 

하지만 자주 이용하는 또는 향후에 또 사용하게 되는 폴더라면 인터넷 브라우저에서 북 마크를 하는 것처럼 Matlab Shortcuts 에 등록을 해 놓고 이용을 하는 것이 훨씬 더 편할 것이다.

 

사용방법은 매우 간단하다.

MATLAB 메인 창의 Shortcuts 에서 마우스 오른쪽 클릭을 한다.

 



NEW shortcuts 을 클릭을 하면 아래와 같이 Editor 가 나온다.

이동하고자 하는 폴더 또는 작업에 맞게 적절히 Label 을 붙인다.

나의 경우엔 J 드라이브의 project 폴더에서 자주 작업을 하므로

J:\project 폴더로 이동하도록 하기 위하여 Callback 에

cd J:\project 적었고 Label 은 Project_folter 라고 적었다.



이제 save 를 하면 MATLAB Shortcuts 창에 아래와 같이 Project_folter 가 등록이 된다.

 

 

   

이젠 저 단추만 클릭 해주면 작업하고자 하는 J:\project 폴더로 이동을 하게 된다.


위에서 설정한 shortcut 내용들은 아래 폴더에 저장이 됩니다. shortcut 을 백업해 놓고 싶은 경우 아래 주소에서 shortcuts_2.xml 파일을 백업해 놓으면 됩니다.

C:\Users\<계정>\AppData\Roaming\MathWorks\MATLAB\R2016a

빠른 실행 메뉴에 보이게 하기 위해서는 MATLABQuickAccess.xml 파일 역시도 백업을 해놓고 사용하는게 좋습니다.

움직이는 그래프를 만들어 보자

매트랩에서 그래프를 동영상으로 만드는 방법은 매우 간단하다.

변수에 따라 for 문을 이용하여 그래프를 그리고 각 프레임을 저장을 한다. 그리고 movie() 함수를 이용하여 실행한다.

아래 코드는 사인 함수의 계수값을 증가시키면 어떻게 되는지를 보여주는 코드이다.  


axx = 0:0.01:2*pi;  

for k=1:16
plot(axx, k*sin(axx));

grid on

axis([0 2*pi -16 16])

M(k)=getframe;

end

movie(M,1)  


avifile() 함수를 이용하여 avi object 를 만들어서 avi 파일로 저장도 가능하다.

아래 결과는 좀더 좁은 간격의 k 값에 대하여 실행한 결과이다. 
 



  1. ㅠㅠ 2012.03.20 01:32

    동영상처럼 나오게 하려고 했는데
    for k=0.1:0.1:16으로 하니 Subscript indices must either be real positive integers or logicals. 에러뜨네요 ㅠㅠ 방법좀 알려주시면 감사하겠습니다!

    • 남성 2012.03.20 02:26 신고

      for k=0.1:0.1:16 로 하면 k 값은 0.1, 0.2, 0.3 ... 16 까지 가는거죠.

      그렇다면 아래 부분의 코드를 보면 M(k)=getframe; 부분에서

      M(0.1), M(0.2), .... M(16) 이렇게 들어가게 되겠죠. 벡터 또는 매트릭스의 인덱스는 양의 정수나 logical 값만 들어가야되니깐 에러가 나게 되는거죠.

      좋은 방법은 아니지만 아래와 같이 코딩할 수도 있습니다.


      axx = 0:0.01:2*pi;

      for k=0.1:0.1:16
      plot(axx, k*sin(axx));

      grid on

      axis([0 2*pi -16 16])

      M(round(k*10))=getframe;

      end

      movie(M,1)


      더 좋은 방법은 k를 1:160 으로 하고

      plot(axx, k./10 .*sin(axx)); 이렇게 바꿔주는게 개인적으로는 더 추천 하는 방식입니다.

  2. 2012.03.20 22:24

    비밀댓글입니다

    • 남성 2012.03.20 22:54 신고

      죄송하지만, k가 커졌다 작아졌다 한다는게 정확히 무슨 말씀인지 모르겠네요. 좀더 구체적으로 말씀해 주시는게 좋을듯 합니다.

      k 값이 1에서 부터 16 까지 갔다가 다시 1로 줄어들게 되는걸 말슴하시는건가요? 아님 계속 k 값이 랜덤 하게 들쭉 날쭉 하게 되게 한다는건가요?

    • 2012.03.21 00:52

      비밀댓글입니다

    • 남성 2012.03.21 01:01 신고

      이렇게 하면 됩니다.


      axx = 0:0.01:2*pi;

      temp_k=[1:16 15:-1:1];

      for k=1:length(temp_k)
      plot(axx, temp_k(k)*sin(axx));

      grid on

      axis([0 2*pi -16 16])

      M(k)=getframe;

      end

      movie(M,1)

  3. 2015.07.30 14:52

    비밀댓글입니다

    • 남성 2015.07.30 23:52 신고

      위 코드에서는 M 의 인덱스가 잘못된것 같네요.

      다음과 같이 해 보세요.

      아래 코드가 방문자님이 의도한게 맞나 모르겟네요.

      wt = 0:pi/50:10*pi;

      x = cos(wt + (pi/180));
      y = wt-wt;
      x1 = wt-wt;
      y1 = cos(wt);

      for c = 0:1:10
      figure
      plot3(x-c, y, wt, x1, y1-c, wt)

      axis tight
      grid on
      M(c+1)=getframe;
      end

      movie(M,1)

    • 2015.07.31 11:08

      비밀댓글입니다

    • 남성 2015.07.31 12:15 신고

      다음과 같이 해보세요.

      wt = 0:pi/50:10*pi;

      x = cos(wt + (pi/180));
      y = wt-wt;
      x1 = wt-wt;
      y1 = cos(wt);

      h=figure(1);
      set(h,'visible','off')

      for c = 0:1:10
      plot3(gca, x-c, y, wt)
      axis([-10 10 -10 10 0 32])
      grid on
      M(c+1)=getframe;
      end

      set(h,'visible','on')
      movie(M,1)

    • 2015.08.01 12:45

      비밀댓글입니다

    • 남성 2015.08.01 13:41 신고

      도움되었다니 다행이네요. 방문해 주셔서 감사합니다. ^^

MATLAB 다양한 능력중에 string evaluation 관련한 함수인 eval 이용해 보자. 가령 아래와 같은 예를 수행해 보고 싶은 경우가 있다.  

A1 이라는 변수에는 1행 1열의 uniform sample 을 저장하고
A2 이라는 변수에는 2행 2열의 uniform sample 을 저장하고
A3 이라는 변수에는 3행 3열의 uniform sample 을 저장하고
…. A100 이라는 변수에는 100 행 100 열의 uniform sample 을 저장하고 싶다고 해보자. 

뭐 변수가 한 두개라면 그냥 몇줄 쓰면 되니 그리 문제될게 없지만, 이 경우처럼 100 개나 되는 변수에 저 값들을 일일이 할당한다는건 비효율 적이라는건 누구나 알 수 있다. 이런 경우 eval() 함수를 이용하면 너무나 간단히 위 문제가 해결 된다.  

for n = 1:100
       eval(['A' num2str(n) ' = rand(n)'])
end  

위 코드가 다다 너무나 간단하다 ㅋㅋ
workspace에 보면 A1~A100 의 변수가 생긴걸 확인 할 수 있을 것이다. 위 코드 중에 num2str 함수는 숫자를 문자로 바꿔주는 함수이고 rand 함수는 0~1 사이의 uniform random variable 을 n 행 n 열로 생성하는 함수이다.

  1. 우웅v 2012.02.10 02:05

    for n=1:10
    eval(['A' num2srt(n) ' = rand(n)'])
    end

    Undefined function 'num2srt' for input arguments of type
    'double'.

    이렇게 하면 에러가나는데 뭐가 문제인지 혹시 아시나요...?

  2. 남성 2012.02.10 08:49 신고

    타이핑 실수네요.

    num2srt() 라는 함수는 없습니다.

    num2str() 로 바꿔 주시면 됩니다.

  3. 초보자 2012.12.26 16:16

    한가지 질문 해도될까요..?
    eval 함수를 사용해서 변수에 데이터를 저장하셨잖아요
    그럼 저 데이터가 들어간 변수를 하나의 cell에 저장할 수는 없는건가요?
    저도 비슷하게 eval 함수를 사용하는데
    하나의 cell에 데이터가 들어간 변수를 넣는게 힘들어서요 ㅠ

    • 남성 2012.12.27 06:43 신고

      eval 을 이용해서 cell 에도 저장을 할 수 있습니다.

      질문의 의도를 정확히 이해한건진 모르겠지만 다음과 같이 저장할 수 있을것 같습니다.

      A={};
      for n = 1:100
      eval(['A{' num2str(n) '} = rand(n)'])
      end

      위와 같이 코딩하면

      A{1}, A{2}, .... A{100} 에 값들이 저장됩니다.

MATLAB 을 이용하여 아래와 같은 조건의 시스템에 대한 시스템 응답을 구해 보자.

 

<code>

 

a=[1 -1 1];     % 피드백 루프의 계수

b=[1 0 0];        % 포워드 루프의 계수

n=(0:40)';        % Sample index

x=inline('cos(2*pi*n/6).*(n>=0)','n'); % x(n) 인라인 함수로 생성

z_i=filtic(b,a,[1 2]); % filter initial condition 설정

y=filter(b,a,x(n),z_i); % filter 함수를 이용한 시스템 응답

stem(n,y), grid on, xlabel('n'), ylabel('y[n]') 


FFT 되는 신호의 중간에 upsampling을 하고 그 결과가 어떻게 되는지를 확인해 보자

다음과 같은 sin() 그래프가 있고 이에 대한 FFT 를 수행해보자
pha = linspace(0,2*pi,32)';
x=sin(pha);
y=fft(x);
figure(1), plot(pha,x,'b:.'), grid on

다음과 같이 x 의 값에 3 sample 씩 upsample 을 하고 그에 대해 FFT 를 해보자. 
xZero = upsample(x,3);
yy=fft(xZero); 
x 대한 fft 결과인 y 값과 upsample 신호에 대한 fft 결과인 yy 같이 그리기 위하여 다음과 같이 x 축을 설정하고 FFT 결과에 대한 절대값 그래프 위상값을 확인해 보자
x_axis=0:length(x)-1;
x2_axis= length(x)/length(xZero)*(0:length(xZero)-1); 
h2= figure('units', 'normalized', 'pos',[0.5 0.5 0.4 0.4]);, plot(x_axis,abs(y),'r:*' , x2_axis, abs(yy),'b-->'), grid on, title('| FFT reselt |'),legend('Original','Zero padded')

h3= figure(
'units', 'normalized', 'pos',[0.5 0 0.4 0.4]);, plot(x_axis,unwrap(angle(y)),'r:*' , x2_axis, unwrap(angle(yy)),'b-->'), grid on, title('angle(FFT reselt )'),legend('Original','Zero padded')

 
위의 결과에서 확인할 있는 것과 같이 upsample 만큼 반복 특성이 발생한다는 것을 확인 있다.

  1. 공대생 2015.05.03 23:19

    감사히 잘 봤습니다!
    실례지만 혹시 upsample 함수가 matlab에 존재하는 함수인가요?
    upsample함수 혹시 만드신거면 소스 좀 여쭤봐도 될까요?

    • 남성 2015.05.03 23:38 신고

      upsample 함수는 아래 주소에서 확인 할 수 있듯이 signal processing toolbox 에 있는 함수 입니다.

      http://kr.mathworks.com/help/signal/ref/upsample.html;jsessionid=a315b53bd2a1bf11daf0b2f6bc42

      혹시 Signal Processing Toolbox 가 없어서 만들어 사용하신다면 그냥 간단하게 만들 수도 있겠지만 다른 코드를 참조 하고 싶다면 아래 주소의 octave upsample 을 참조 하시는게 좋을것 같네요.

      http://octave-signal.sourcearchive.com/documentation/1.0.9/upsample_8m-source.html

    • 공대생 2015.05.03 23:47

      이렇게 늦은 시간에 답을 친히 해주시고,
      단순한 질문에 이렇게 상세하게 답변해주셔서 감사합니다.
      좋은 밤 되시길 바라겠습니다^^

    • 남성 2015.05.04 00:31 신고

      네~ 공대생 님도 좋은 밤 되세요~ 방문해 주셔서 감사합니다. ^^

wav 파일을 읽어 들이고 이에 대하여 MATLAB 프로그램에서 실행 시켜 보자.

 같은 폴더 내에 잇는 a.wav 라는 파일을 파일을 읽어들이는 건 아래 함수와 같이 하면 된다.  

[y,Fs,bits]=wavread('a.wav');  

그때 함수로부터 받을 있는 값들은 음성 값인 y 값과 샘플링 주파수 Fs , 그리고 샘플당 비트 수인 bits 등을 받는다.

이렇게 읽어 들인 소리 파일을 재생하는 또한 아래와 같이 매우 간단한다.

sound(y, Fs);


하지만 함수는 경험상 봤을 , 중간에 멈춰 지지가 않는다.
만약 노래 등을 중간에 멈추고 싶을때는 난감한 경우가 많다.
    
이런경우를 대비하여 다음과 같이 audioplayer object 를 이용할 수 있다.
 
ob= audioplayer(Y, Fs) 
 
위와 같이  줄만 써주면 audioplayer  object 만들어 진다. 
 
play(ob), pause(ob), stop(ob) 등과 같은 다양한 method 제공해 준다.
 
이외에도 다양한 함수등을 제공해 주니 Help 참조하기 바란다.


변수의 크기를 알 수 있는 함수는 3가지 정도가 있다 .

다음과 같은 x 라는 변수에 대하여 함수별로 그 특징을 알아보자.  

x=rand(2,7)  

x =

0.8147 0.1270 0.6324 0.2785 0.9575 0.1576 0.9572

0.9058 0.9134 0.0975 0.5469 0.9649 0.9706 0.4854

 size(x)
ans = 2 7

결과와 같이 size() 함수는 Dim 별로의 크기를 나타내 준다 

length(x)
ans =  7      
Dim. 크기중 수를 나타내 준다

numel(x)
ans = 14

전체 element 개수를 나타낸다.

MATLAB 으로 PDF 를 그려보자.

PDF(Probability Density Function)라는게 어케 보면 Histogram이랑 개형은 같지만

엄밀하게 얘기하면 좀 다르다고 할 수 있다.

그 차이는 바로 Normalization에 있다.

PDF 의 특징은 그 적분값이 1 이어야 하므로 당연히 Histogram 을 그린후 그 넓이를 Normalization 을 해야한다.

간단하게 [0 100] 구간의 Uniform pdf 를 그려 보자.

MATLAB editor 에 다음과 같이 타이핑 하고 실행해 보자

<MATLAB CODE>


N = 3000;

Uniform_sample=rand(N,N)*100; % 0~ 100 사이의 Uniform sample 을 발생시킨다.

[pdf,X]= hist(Uniform_sample(:), 500); % hist 함수를 이용하여 구간별로 쌓아준다

resol=X(2)-X(1);             % resolution을 계산한다.

pdf= pdf/(N^2.*resol); % Normalization 해 준다.

figure(1), bar(X,pdf), grid on, title('Uniform PDF') % 그래프를 그린다.

그래프를 보면 딱 봐도 넓이가 100 * 0.01 =1 이라는 것을 알 수 있다. 이게 제대로 된 PDF 이다.

 

그럼 좀더 나아가서 rand() 함수를 이용하여 Gaussian sample 을 만들고 PDF 를 그려 보자.

 

<MATLAB CODE>

 

N = 3000;

std = 2; % 표준편차

std_normal=sqrt(-2*log(rand(N))).*cos(2*pi*rand(N)); % 평균이 0 이고 분산이 1인 Gaussian Sample

NormalSample=std*std_normal; % 표준편차가 2인 Gaussian Sample

[pdf,X]= hist(NormalSample(:), 500); % hist 함수를 이용하여 구간별로 쌓아준다

resol=X(2)-X(1); % resolution을 계산한다.

pdf= pdf/(N^2.*resol); % Normalization 해 준다.

figure(2), bar(X,pdf), grid on, title('Gaussian PDF') % 그래프를 그린다.


FFT 가 되는 신호에 대한 zero padding 효과가 어떻게 되는지를 알아보자

다음과 같은 sin() 그래프가 있고 이에 대한 FFT 를 수행해보자

pha = linspace(0,2*pi,32)';

x=sin(pha);

y=fft(x);

figure(1), plot(pha,x,'b:.'), grid on

아래 코드와 같이 x 값의 뒤에 0 을 넣어 보자. 이를 zero padding 이라 한다.

x 값의 크기 만큼, 이 예에서는 32 만큼의 0을 뒤에 넣어 보자. 그리고 이에 대해 fft 를 해보자


 xZero = [x ; zeros(32,1)];

yy=fft(xZero);


 x 대한 fft 결과인 y 값과 zero padding 신호에 대한 fft 결과인 yy 같이 그리기 위하여 다음과 같이 x 축을 설정하고 FFT 결과에 대한 절대값 그래프를 확인해 보자.


 x_axis=0:length(x)-1;

x2_axis= length(x)/length(xZero)*(0:length(xZero)-1);

 figure(2), plot(x_axis,abs(y),'r:*' , x2_axis, abs(yy),'b-->'), grid on, title('| FFT reselt |'),legend('Original','Zero padded')


 

 

다음과 같이 Original 신호의 FFT 결과의 중간 중간에 점이 하나씩 생긴 것을 확인 있을 것이다.

, fft 이전 신호의 부분에 zero padding 하게 되면 FFT 이후에 그만큼 oversample 된다는 것을 있다.

나아가서 예에서 Original 신호의 중간에 두개씩이 생기게 하고 싶다면 개의 zero 넣어줘야 할까에 대해서도 생각해 있다.

x 뒷부분에 64개의 0 추가 준다면 점이 두개씩 생겨나게 된다.


+ Recent posts