(Matlab) DFT 함수 구현

지난주에 함수 사용법을 배우던 날 결석을 해서 

이 것을 만들때 함수 만드는거에서부터 막혀서 눈치보면서 코딩을 했다.


일단 함수를 만드는 법 부터 설명

function X=mydft(x,N) 라고 선언을 해주고 내용을 작성하면 된다

X는 return value

x와 N은 input value이다.(parameter라고 보면 되지)

중요한것은 m-file로 저장을 따로 해야 하며

저장시에는 함수 이름 그대로 mydft라고 저장을 해야한다는거!!!

아!!또한 같은 디렉토리안에 파일들이 존재해야한다!!



자아...이제 DFT에 대한 간단한 설명을 해보자.

일단 MATLAB에서는 DTFT를(DTFS도 커버) 사용한다.  

그렇지만 무한대는 MATLAB에서 사용을 할 수 없다.

즉, 이 식을 컴퓨터로 표현하는 것은 불가능 하다는 것

그래서 식을 이렇게 바꾼 것을 DFT(Discrete Fourier Transform)이라고 한단다.

일단 DFT 함수의 기능을 테스트하기 위한 코드
(이 코드에 대해서도 연구할게 있다. 그건 일단 나중에 설명)

100MHz의  cos 그래프를 테스트 모델로 만들어보았다.
샘플링 주파수를 2로 하니까 그래프가 저렇게 나오는구나. 암튼 고고!!

% My DFT Function Test
f0=100e3;     % 100Mhz
fs=2*f0;      % sampling frequency
Ts=1/fs;
t=0:Ts:0.01;

x=cos(2*pi*f0*t);
N=500;

y=mydft(x,N);
% plot(y); % complex값은 그래프로 못찍기때문에 제대로 결과값이 나오지 않는다

x_mag=abs(y);
plot(0:N-1,x_mag);    % 이렇게 해주어야 한다

자 이제 이 식을 Matlab에 mydft라는 함수를 만들어서 그 안에 집어 넣어보자!

function X=mydft(x,N)
X=zeros(1,N);
for k=0:N-1
    for n=0:N-1
        X(k)= X(k) + x(n)*exp(-j*(2*pi/N)*k*n);
    end
end

위의 식은 DFT식에 그대로 적용시켜서 넣어 만들어본 코드다.

외관상은 완벽하다. 그럼 한번 돌려볼까?

에러가 난다. 이유는...?????

제일 윗 줄을 보면 X(0)은 인정이 되지 않는다고 한다.

그렇지. 매틀랩은 철저한 행렬을 바탕으로 만들어진 프로그램이다.

행렬을 첫번째 값은 1부터 시작하겠지!!

즉 X라는 행의 0번째 값은 존재하지 않으니까 에러가 난다는 말씀!!

그럼 처음 입력 되는 곳을 0이 아니라 1부터 시작해야겠네???

간단하네~!!

식을 이렇게 한번 바꾸어 보자.

function X=mydft(x,N)
X=zeros(1,N);
for k=0:N-1
    for n=0:N-1
        X(k+1)= X(k+1) + x(n+1)*exp(-j*(2*pi/N)*k*n);
    end
end

호~위와 같은 그래프가 나오네!!!

그렇지만 교수님께서는 이 코드가 완벽하지 않다고 하신다.

어디가 부족하지??? 

N-1(노란부분)을 이렇게 바꾸라고요?

function X=mydft(x,N)
X=zeros(1,N);
for k=0:N-1
    for n=0:length(x)-1
        X(k+1)= X(k+1) + x(n+1)*exp(-j*(2*pi/N)*k*n);
    end
end

이것 역시 컴퓨터와 수식의 차이점이라고 한다.

N의 값 500개의 샘플만으로는 부족하다는 것??

솔직히 교수님이 이것에 대해 설명을 하셨던 순간 

잡생각을 해서 지금도 왜인지 잘 와닿지 않는 돌머리...

아무튼 이 것과 바로 위의 그래프를 비교해 보면 값이 크게 다르다.

샘플의 양이 현저히 많아 졌기 때문에 값도 약 4배가량 차이가 난다.

여기서 하나의 의문점이 생긴다.

신호와 시스템에서 배우기로는 분명 cos 함수를 FT하게 되면

0을 기준으로 양 옆에 간격을 두고 두개가 뽀죡하게 튀어 나와야 하는데.

하나뿐이 나오지 않는다??

샘플링 주파수를 높이면 되네?왜지?

샘플링 주파수를 정확하게 두배를 하게 된다면

다른 주파수 하나는 옆에 아슬아슬 걸려서 반만 나온다고 한다!!

아래는 매틀랩 카페에서 어떤 분의 설명!!

사실 지금 주파수분석을 해보는 대상신호가 한 가지 주파수성분만 가지고 있으므로 원래 하나만 나와야 하는데 folding frequency (sampling frequency의 절반)를 중심으로 대칭되게 2개로 보이는 것이죠. 마치 거울의 허상처럼 .. 그런데 우리가 물체를 거울에 갖다 붙여버리면 물체와 허상이 완전히 겹쳐서 하나가 되어 버리지 않겠습니까? 바로 우리가 관찰하려고 하는 신호의 최고주차수의 2배보다 Sampling주파수가 더 높아야 한다는 Sampling Theorem을 아슬아슬하게 만족시키지 못하는 경우중 하나입니다. 마치 2차방정식의 근이 원래는 2개이지만 경우에 따라 그 두 근이 겹쳐저 중근이 되는 것처럼 말이죠.

그래서 보통 2.2배를 한다고!!


어찌되었든 샘플링 주파수를 2배로 하지 않고 2.2배로 조금 더 높여보면

아래의 그래프가 출력 된다.

이제 mydft 함수 자체에는 문제가 없다.

아까 가장 위에 있던 본 테스트 소스를 고쳐서

0을 중심으로 대칭적인 그래프를 만들어보자.

일단 이 그래프의 값이 180이나 나오게 된다.

원래대로라면 1이고 이 값을 둘이서 나눠먹어 각각 0.5 씩 출력되어야 하는데...

이렇게 높은 숫자가 나오는 이유는??

mydft 함수에서 x의 개수만큼 계속해서 더해 버렸기 때문이다.

즉 본코드에서 x의 개수만큼 다시 나누어 주어야 한다는 말이다.

그래서 코드를 다시 수정해서 출력을 해보면


% My DFT Function Test-first
f0=100e3;     % 100Mhz
fs=2.2*f0;      % sampling frequency
Ts=1/fs;
t=0:Ts:0.01;

x=cos(2*pi*f0*t);
N=500;

y=mydft(x,N);
% plot(y); % complex값은 그래프로 못찍기때문에 제대로 결과값이 나오지 않는다

x_mag=abs(y)/length(x);
plot(0:N-1,x_mag);  

0.5가 안나오네...




그래서 샘플링 주파수를 4배로 해서 다시 그래프를 보니

% My DFT Function Test-first
f0=100e3;     % 100Mhz
fs=4*f0;      % sampling frequency
Ts=1/fs;
t=0:Ts:0.01;

x=cos(2*pi*f0*t);
N=500;

y=mydft(x,N);
% plot(y); % complex값은 그래프로 못찍기때문에 제대로 결과값이 나오지 않는다

x_mag=abs(y)/length(x);
plot(0:N-1,x_mag);  

0.5가 나온다!!

샘플링 주파수에 따라서 진폭이 변한다?????

이건 말이 안되는데???

왜지???????????



http://cafe.naver.com/matlablove

이 카페에 물어본 후의 어떤분의 답변이다.

mydft내에서 length(x)대신 N이 들어가면 적어도 두번째 경우처럼 샘플링주파수를
신호주파수의 정수배로 잡았을 때는 MATLAB 내장함수인 fft()와 동일하게 님이
예상했던 것처럼 spectrum의 진폭이 0.5로 나오죠.
그런데 첫번째 경우에는 샘플링주파수가 신호주파수의 정수배가 아니므로
소위 spectral leakage현상이 발생하여 spectrum이 선명하게 나오지 않음으로 해서 정확히 
0.5가 나오지 못하는 겁니다. 그런데 대부분의 경우에는 신호주파수를 몰라서
그걸 알기 위해서 DFT를 취하는 것이므로 사실 첫 번째 경우가 일반적인 spectrum입니다.