오늘은 MATLAB 에서 sample rate conversion 방법에 대해 설명하려 합니다.

 

저는 보통 음원 파일에 대해 sample rate conversion 을 수행할때는 아래 포스팅에서 소개한 ffmpeg 이나 sox 등을 주로 사용하지만 MATLAB에서도 sample rate conversion 을 수행해야 할 경우가 종종 있습니다.


2013/10/21 - [유틸] - Windows 에서 sox 로 음악 파일 변환하기


2014/03/07 - [유틸] - ffmpeg 을 사용하여 rmvb 파일을 avi 파일로 변환하기

 

MATLAB 에서 sample rate conversion 을 할 때는 보통 resample 함수를 사용합니다. Resample 함수는 MATLAB Signal Processing Toolbox 내에 포함되어 있으므로 resample 함수를 사용할 수 없다면 MATLAB command window 에 ver 명령어를 쳐서 Signal Processing Toolbox 가 설치되어 있는지 확인 해 보시기 바랍니다.

 

일반적으로 Resample 함수는 다음과 같이 사용합니다.

 

y = resample(x,p,q)

 

x 라는 입력 신호에 대해 p/q 배로 sample rate 를 변경 하는 겁니다. p, q 값은 당연히 양의 정수 이어야 합니다.

 

보통 sample rate conversion 을 수행할때는 anti-aliasing 필터링을 하게 되는데~ 이러한 과정도 내부에서 자동으로 수행되기 때문에 사용자가 굳이 신경 쓸 필요가 없습니다.

 

아래 주소에서 resample 함수의 help를 보면 기본적으로 firls을 사용하여 필터를 디자인 하고 Kaiser window 를 사용하고 윈도우의 beta 값은 5를 사용 한다고 하는군요.

 

http://www.mathworks.co.kr/kr/help/signal/ref/resample.html

 

anti-aliasing 필터의 계수를 변경하고 싶다면 아래 원형의 b 값에 원하는 필터의 coefficients 값을 넣어주면 됩니다.

 

y = resample(x,p,q,b)

 

resample 함수를 사용하여 8192 Hz 의 음원을 CD 음질인 44100 Hz 로 변경해 보죠~

 

다음과 같이 MATLAB command window 에 명령어를 칩니다. chirp 이라는 음원은 MATLAB 에서 제공하는 음원으로soundsc 함수를 통해 재생해보면 새 소리가 납니다.

 

 

>> load chirp % chirp 신호 load

>> Fs % sample rate 확인

 

Fs =

 

8192

 

>> y2=resample(y,44100,Fs); % 8192 에서 44100 Hz 로 sample rate conversion

>> soundsc(y,Fs) % 8192 Hz 음원 소리 확인

>> soundsc(y2,44100) % 44100 Hz 음원 소리 확인

 

위 예제를 실행해보면 8192 Hz 음원과 44100 Hz 음원의 소리가 같은 것을 확인 할 수 있습니다.

 

하지만 각 데어터의 길이를 확인해 보면 8192 Hz 음원인 y 는 13129 인데 비해 44100 Hz 인 y2 음원의 길이는 70678 라는 것을 확인 할 수 있습니다.

 

>> length(y)

ans =

 

13129

 

>> length(y2)

ans =

 

70678

 

이는 8192/44100 = 13129 /70678 = 0.1858 이기 때문입니다.


전자과에서 수학을 접하다 보면 sinc 함수를 접하게 됩니다. 특히 신호처리 과목을 듣는 사람들이라면 거의 백퍼 접하게 되는게 sinc 함수 입니다.

 

Sinc 함수는 아래 주소에 설명이 잘 나와있습니다.

 

http://ko.wikipedia.org/wiki/%EC%8B%B1%ED%81%AC%ED%95%A8%EC%88%98

 

sinc 함수를 푸리에 변환하면 직사각형 함수(Rectangular Function)가 되게 되고~ 처음 이 사실을 알고 참~ 신기하다는 생각을 한적이 있습니다.

 

오늘은 위에 설명한 사실을 MATLAB 을 사용하여 확인 해 보려 합니다.

 


MATLAB 에서 sinc 함수는 Signal Processing Toolbox에 들어있습니다. Signal Processing Toolbox가 없는 분들은 간단한 함수이니 위 수식대로 만들어서 사용해도 될 겁니다.

 

다만 x=0 인 경우 분모도 0 이 되므로 이 부분만 주의해서 사용하시면 됩니다.

 

따라서 다음과 같이 작성 할 수 있습니다. 함수명은 기존 sinc 함수와의 충돌을 방지하기 위해 sinc_ 로 했습니다.


function y=sinc_(x)

i=find(x==0);

x(i)= nan;

y = sin(pi*x)./(pi*x);

y(i) = 1;

 

간단하게 x=0 인 경우만 아래 포스팅에서 배웠던 NaN(Not-a-Number)을 넣어줬다가 그 값만 1 로 바꿔주는 거죠. Nan 이 아니라 아무 값이나 넣고 싶은 값을 넣어도 결과는 같고~ 그냥 x(i)= nan; 을 없애도 상관 없습니다. 


2011/04/20 - [programming language/MATLAB] - MATLAB NaN

 

그럼 sinc 함수를 그려보죠~ 


x=-6:0.1:6;

y=sinc_(x);

 

figure,

plot(x,y), grid on

 

 

 

다음으로 위에서 얘기했던 sinc() 함수의 푸리에 변환이 구형함수가 되는지 확인해 보죠~ 구형 함수를 나타내기 위해 좀 긴~~ sinc 함수의 x 축을 좀 길게 잡습니다.  


x=-80000*pi:0.1:8000*pi;

y=sinc_(x);

ffty=fft(y);

 

fftShifted=fftshift(ffty);

plot(abs(fftShifted)), grid on

 

다음과 같이 sinc 함수의 푸리에 트랜스폼 결과가 직사각파가 되는 것을 확인 할 수 있습니다. Sinc 함수 … 참~~ 신기하죠?? 나만 신기한가…



MATLAB 을 이용한 통신 시뮬레이션 코드를 찾다 보니 아래 주소의 블로그를 발견했다.

 

http://www.dsplog.com/

 

대단히 자세한 이론 설명과 훌륭한 MATLAB 예제 코드들을 보면 디지털 통신을 공부하는데 정말 부족함이 없어 보인다.

 

학부 시절에 FSK(Frequency Shift Keying)를 공부하면서~ BER(Bit Error Rate) 시뮬레이션 까지는 못해봤던 것 같은데~

 

아래 주소에 Frequency Shift Keying BER 시뮬레이션 코드가 있었다.

 

http://www.dsplog.com/2007/08/30/bit-error-rate-for-frequency-shift-keying-with-coherent-demodulation/

 

이론적 설명도 자세하고 세부 코드까지 포함되어 있어서 굉장히 도움이 많이 되는 포스팅이라는 생각이 든다.

 

아래 코드는 저자의 시뮬레이션 코드에 스펙트럼 확인을 하기 위해 몇 가지 코드를 추가한 예이다. 내가 추가한 부분은 빨간색으로 표시했다.

 

아래 코드에서 T 변수의 주석에 symbol duration 이라고 되어 있어서 약간 헷갈렸지만~ 변수 T 가 Sampling Rate(Fs)이다.

 

암튼 저자의 FSK simulation 코드를 보면 Coherent FSK 복조까지 명확하게 알 수 있어서 FSK 를 공부하는 분들이 참조하기에 대단히 좋은 코드라는 생각이 든다.

 

아래 코드에서 내가 스펙트럼을 확인하기 위해 사용한 pwelch() 함수를 사용하기 위해서는 MATLAB Signal Processing Toolbox 가 있어야 한다.

 

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

% 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        : 05 June 2008

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

 

% Simple Matlab example of binary Frequency Shift Keying using

% coherent demodulation

 

clear

N = 10^5 % number of bits or symbols

T = 8; % symbol duration

Fs=T; % Sampling Rate

t = [0:1/T:0.99]; % sampling instants

tR = kron(ones(1,N),t); % repeating the sampling instants

 

Eb_N0_dB = [0:11]; % multiple Eb/N0 values

 

for ii = 1:length(Eb_N0_dB)

% generating the bits

ip = rand(1,N)>0.5; % generating 0,1 with equal probability

freqM = ip+1; % converting the bits into frequency, bit0 -> frequency of 1, bit1 -> frequency of 2

freqR = kron(freqM,ones(1,T)); % repeating

x = (sqrt(2)/sqrt(T))*cos(2*pi*freqR.*tR); %generating the FSK modulated signal

 

% noise

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

 

% coherent receiver

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

op1 = conv(y, sqrt(2/T)*cos(2*pi*1*t)); % correlating with frequency 1

op2 = conv(y, sqrt(2/T)*cos(2*pi*2*t)); % correlating with frequency 2

 

% demodulation

ipHat = [real(op1(T+1:T:end)) < real(op2(T+1:T:end))]; %

nErr(ii) = size(find([ip - ipHat]),2); % counting the number of errors

 

end

simBer = nErr/N;

theoryBer = 0.5*erfc(sqrt((10.^(Eb_N0_dB/10))/2)); %theoretical BER

 

close all

figure

semilogy(Eb_N0_dB,theoryBer,'b-');

hold on

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

axis([0 11 10^-4 0.5])

grid on

legend('theory:fsk-coh', 'sim:fsk-coh');

xlabel('Eb/No, dB')

ylabel('Bit Error Rate')

title('Bit error probability curve')

 

 

%%

 

figure

pwelch(y,1024,512,1024,Fs)

 

위 MATLAB 코드를 실행시켜보면 다음과 같이 Coherent Frequency Shift Keying BER 과 FSK Power Spectral Density 그래프를 확인 할 수 있다.


< Coherent Frequency Shift Keying BER >

< Coherent FSK Power Spectral Density >

 

정확히 주파수 1, 2 Hz 에 톤이 뜨는 것을 확인 할 수 있다.


  1. Genius 2014.02.03 13:53

    좋은 자료 감사합니다.
    담아갈께요. ^^

    • 남성 2014.02.04 10:38 신고

      방문 감사합니다. 공지에 밝혔다 시피 링크만 허용합니다.

  2. ryu 2014.06.19 17:29

    conv에서 t 설정은 어떻게 해주는건지 잘 모르겠어요ㅠㅠ

    • 남성 2014.06.20 19:58 신고

      질문이 무슨 뜻인지 잘 몰겠네여... 위에 보시다시피 t = [0:1/T:0.99];

      t 는 1/T 간격으로 해 준건데요....

  3. chobo 2016.06.10 21:35

    안녕하세요 좋은자료 잘봣습니다.
    혹지 저기서 fs 는 뭘 말하는 건가요?
    FSK에서의 sampling rate가 의미하는걸 잘 모르겟어서..

    • 남성 2016.06.12 21:26 신고

      FSK 에서의 의미라기 보다는 디지털에서의 의미이라고 볼 수 있을것 같네요. 초당 Fs 샘플이 있다는 거죠.

요즘 영상이나 음성과 같음 미디어 컨텐츠들이 많이 사용되고 있습니다.

 

이런 미디어 파일들은 데이터 량이 굉장히 크기 때문에 손실 압축 방식으로 그 데이터를 줄이는 압축 기술들이 많이 사용되는데~

 

이런 기술에 많이 이용되는 것이 바로 DCT(Discrete cosine transform) 라고 합니다.

 

mp3, jpg 같은 파일 형식들이 다~~ DCT 를 사용한다고 하니깐 정말 우리 생활과 너무나 밀접한 기술이라 할 수 있을 것 같네요.

 

DCT 위키 피디아 : http://en.wikipedia.org/wiki/Discrete_cosine_transform

 

위 주소의 내용을 보면 DFT 는 periodic 신호의 비연속 특성 때문에 고주파가 많이 올라오는 반면에 DCT 는 연속적이어서 고주파 성분이 적고 ~ 그에 따라 저주수에 신호 파워가 몰리는 특징이 있더군요.

 

따라서 특정 Threshold 보다 작은 고주파 신호를 없앤다 해도 DFT 에 비해 신호의 왜곡이 덜 할 것이라는 것을 예상해 볼 수 있습니다.

 

DCT 도 DFT 와 유사하기 때문에 FAST DCT 알고리즘이 존재 합니다.

 

FAST DCT 에 대한 내용들이 정리되어 있는 페이지 : http://fourier.eng.hmc.edu/e161/lectures/dct/node2.html

 

오늘은 위 주소의 내용들을 참조하여 FAST DCT 에 대해 간단히 소개해 보려 합니다.

 

일반적으로 흔히 말하는 DCT 수식은 다음과 같고 DCT-II 라고 합니다.

 

 

Inverse DCT 는 DCT-III를 말하며 다음 수식과 같죠.

 

 

FAST DCT 의 과정은 다음과 같이 세 단계로 구성되더군요.

 

FAST DCT 과정

 

 

 

 

FAST IDCT 과정

 

 

 

 

이론을 알아 봤으니깐 간단하게 MATLAB으로 확인해 볼까요~

 

MATLAB Signal Processing Toolbox 가 있으신 분들은 dct(), idct() 라는 함수를 사용하실 수 있습니다. 두 함수들은 DCT,IDCT matrix 를 orthogonal 하게 만들어 주기 위해서 특정 scale factor 들을 곱해 줬는데요.

 

위에서 서명했던 DCT-II, DCT-III 수식을 MATLAB 의 orthogonal 한 DCT 수식들과 같게 만들기 위해서는 다음과 같이 scale factor 를 곱해 주어야 합니다.

 

DCT-II 의 경우는 X(0) 에 을 곱해 주고 X 전체에 를 곱해 주면 되고~

 

DCT-III 의 경우는 X(0) 에 을 곱해 주고 전체에 곱해 주면 MATLAB 내장 함수들과 동일한 결과를 얻을 수 있습니다.

 

코드는 다음과 같습니다.


 

 

실행 결과에 대해 command 창에서 확인해 보면 에러가 10-15 정도니까~~~ 0 이나 마찬가지라는 것을 확인 할 수 있죠~

 


+ Recent posts