[MATLAB] FFT 돌리기
MATLAB

[MATLAB] FFT 돌리기

 

 

 

0. 관련 글 리스트

 

 

 

 

 

 

MATLAB에서 FFT를 계산하는 방법 중 인터넷에 가장 널리 알려진 code는 PositiveFFT이다.

그냥 MATLAB에 내장된 기본적인 fft 함수만 이용하면

우리가 흔히 알고있는 spectrum이 (회로쟁이 기준) 정상적으로 뽑혀나오질 않기 때문에 추가적인 작업이 필요한데,

PositiveFFT 라는 것은 이를 위한 것이다.

​

먼저, sinewave를 이용해 기본적인 fft를 돌려보면 아래와 같다. 4Hz sinewave를 예시로 잡았다.

 

<<코드가 보이지 않는다면 오른쪽 아래의 달모양 버튼으로

테마를 다크모드에서 화이트모드로 바꿔주세요>>

 

fin = 4                 % Frequency of the sine wave
Fs = 100               	% Sampling rate
Ts = 1/Fs              	% Sampling time interval
t = 0:Ts:100-Ts        	% Data duration
                        % 100s of 4Hz sinewave
y = 2*sin(2*pi*fin*t)	% The sine curve

 

이렇게 만들어낸 y vector를 fft function을 이용해 값을 뽑아보면 이런 식이 된다.

 

>>stem(abs(fft(y)))

 

Y-Axis 값도 안 맞고 X-Axis도 frequency가 아니라 그저 sample 개수만큼의 정보를 보여줄 뿐이다.

이 FFT 결과는 negative frequency부터 positive frequency 까지의 결과를 plot 한 것에 불과하기 때문이다.

FFT의 특성 상, 전체 범위 10,000 중 절반부터 (5,000) 그 이상은 negative frequency 부분이 반복되어 나타난 결과이다.

Sinewave를 fft 했는데 peak이 두 개가 뜬 것으로 이를 판단할 수 있다.

​

이런 경우, PositiveFFT를 사용할 수 있다.

아무데나 폴더 하나 파고 새로만들기 -> 라이브함수 (.mlx 파일)를 선택해서 아래 코드를 기입한다.

 

 

 

 

 

 

1. Positive FFT 함수

% Make a New Function File (.mlx)
function [freq, X] = positiveFFT(x,Fs)
    N = length(x);              % Get the Number of Points
    k = 0:N-1;                  % Create a Vector from 0 to N-1
    T = N/Fs;                   % Get the Minimum Frequency
    freq = k/T;                 % Create a Vector fot the Frequency Range

    X = fft(x)/N;				% Normalize the FFT Data
    cutOff = ceil(N/2);         % Only want the First Half of the FFT
                                % Otherwise it is Redundant
    X = X(1:cutOff)*2;          % Take the First Half of Spectrum
    freq = freq(1:cutOff);
end

 

PositiveFFT는 input vector(​x)와 sampling frequency(Fs)를 argument로 받는다.

주어진 Fs로 표현할 수 있는 max. frequency의 sinewave는 Fs이다.

Successive한 두 sample 사이에서 sinewave가 한 주기를 겪으면 그게 Fs이기 때문이다.

그 이상으로 빠른 sinewave는 주어진 Fs 가지고는 표현할 방법이 없다.

​

위 코드를 보면, 우선, x라는 vector를 받아와서 sample 수가 몇 개인지 세어본다.

이 짓을 하는 이유는 FFT의 bin width (Min. Freq.)를 알아내기 위함이다.

(+bin size라고 표현하면 bin의 개수를 말하는건지 min. frequency를 말하는건지 모호하므로 width라 칭하겠음)

Sample 수가 N개일 때, 표현 가능한 min. frequency는 Fs/N이다. 더 상세한 자료는 아래 논문 참고.

 

9086회 인용 ㄷㄷ

 

On the use of windows for harmonic analysis with the discrete Fourier transform

This paper makes available a concise review of data windows and their affect on the detection of harmonic signals in the presence of broad-band noise, and in the presence of nearby strong harmonic interference. We also call attention to a number of common

ieeexplore.ieee.org

 

Sample의 시작부터 sample의 끝까지 sinewave가 한 개의 period를 겪으면 그게 min. frequency이며,

이것보다 더 느린 frequency의 sinewave는 표현이 불가능하다.

k는 N개의 vector를 뜻한다.

T는 min. frequency (bin size)의 period를 뜻한다.

freq는 간격이 Fs/N인 bin이 N개 존재하는 frequency range를 뜻한다.

즉, 아래와 같은 range에서 N개의 array가 Fs/N의 간격으로 존재하는 것이다.

 

(모바일에서는 수식들이 LaTeX code로 표기됩니다. 데스크탑 버전으로 봐주세요.)

 

$$ \left[0, \left( {N-1 \over N} \right)F_{S}\right] $$

 

요까지가 FFT의 X-Axis(Frequency)를 만들어주는 과정이다.

​

그 다음 코드를 보면, X는 FFT의 output을 뜻하는데, sample의 개수인 N으로 나누어주어야 normalize가 되면서 제대로 된 값이 나온다.

cutOff는 negative frequency part를 날려주기 위해 절반영역, 즉, Nyquist band 까지만 선택을 하겠다는 뜻으로 사용한다.

또한, ceil function을 이용해 first half를 잘라내면서 amplitude 정보가 절반 날아가버리므로 x2를 해준다.

마지막으로 X (FFT의 output)와 freq (FFT의 X-Axis)를 cutOff에 맞춰 first half만 잘라내는 것으로 function은 끝이 난다.

 

 

 

 

 

 

2. FFT 돌리기

그 후, PositiveFFT.m 파일과 같은 폴더에 위치한 곳에 라이브 스크립트 (.mlx) 파일을 아무렇게나 하나 만들고

아래 코드와 같이 FFT를 돌릴 수 있다.

Sinewave 만드는 code는 이 글 제일 위에 있는거 그대로 갖다 붙인 것이다.

 

clc                 		% Clear the Logs in Command Window
clear               		% Clear the Data in Cache
close all           		% Close the Windows (e,g, figures)

%% Generate Sinewave
fin = 4                 	% Frequency of the Sinewave
Fs = 100               		% Sampling Rate
Ts = 1/Fs              		% Sampling Time Interval
t = 0:Ts:100-Ts        		% Data Duration
                        	% 100s of 4Hz Sinewave
y = 2*sin(2*pi*fin*t)		% The Sinewave

%% Calculate FFT
y = y'						% Convert to Row-Wise Vector
[x_freq, y_FFTout] = PositiveFFT(y, Fs)
x_freq = x_freq'							% Convert to Row-Wise Vector
y_FFTout_dB = 20.*log10(abs(y_FFTout))		% Convert Y-Axis in dB Scale

%% Plot FFT
figure(1)
subplot(2,1,1)
stem(x_freq, abs(y_FFTout))			% Plot the Magnitude in Linear X-Axis
subplot(2,1,2)
semilogx(x_freq, y_FFTout_dB)		% Plot the dB Scale Magnitude in Log X-Axis

나름대로 깔끔하게 잘 나오며, amplitude도 내가 설정했던 2가 그대로 잘 나온다.

2는 dB scale로 6.02니까 정확하다. 

 

 

 

 

 

 

3. Window 씌워서 FFT 돌리기

사실, 실제 data를 이용해 FFT를 돌릴 때는 sample 개수가 2N을 만족해야 올바른 값을 얻을 수가 있다.

또한, 우리가 실제 data를 다룰 때는 time duration을 위와 같이 100초면 100초 딱 정해놓는 접근이 아니라

data의 sampling rate가 정해졌을 때 sample 개수를 정하는 방식으로 FFT를 돌리기 때문에,

time-axis를 sample 개수를 기준으로 다시 설정하고 FFT를 돌릴 필요가 있다.

 

이에 맞춰, 아래와 같이 code를 수정해 보았다.

Number of sample (Nsamp)는 앞서 진행한 10,000과 비슷한 값인 8,192로 정했다.

 

clc                 						% Clear the Logs in Command Window
clear               						% Clear the Data in Cache
close all           						% Close the Windows (e,g, figures)

%% Generate Sinewave
fin = 4                 					% Frequency of the Sinewave
Fs = 100               						% Sampling Rate
Ts = 1/Fs              						% Sampling Time Interval
Nsamp = 2^13								% 2^13 = 8,192
t = linspace(0, Nsamp*Ts, Nsamp)			% Time-Axis
y = 2*sin(2*pi*fin*t)						% The Sinewave

%% Calculate FFT
y = y'										% Convert to Row-Wise Vector
[x_freq, y_FFTout] = PositiveFFT(y, Fs)
x_freq = x_freq'							% Convert to Row-Wise Vector
y_FFTout_dB = 20.*log10(abs(y_FFTout))		% Convert Y-Axis in dB Scale

%% Plot FFT
figure(1)
subplot(2,1,1)
stem(x_freq, abs(y_FFTout))					% Plot the Magnitude in Linear X-Axis
subplot(2,1,2)
semilogx(x_freq, y_FFTout_dB)				% Plot the dB Scale Magnitude in Log X-Axis

 

​분명 frequency는 4Hz로 맞긴 한데, amplitude가 1.75로 이상하게 잡히면서 뭔가 잘못되었음을 알 수가 있다.

또한, dB-scale에서 보이는 noise floor도 -50dB 근처인데, 앞선 그래프는 -300dB 근처였으므로 이 또한 문제가 있다는 것을 의미한다.

마치 skirt 처럼 생긴 저 현상을 spectral leakage라고 부른다.

Oscillator의 phase noise를 표현할 때도 skirt란 표현을 쓴다.

저걸 해결하기 위한 방법으로 흔히 알려진 것이 window이다.

 

아래와 같이, PositiveFFT 함수에 Hann Window를 위한 code를 추가하고

fft 파트의 code를 수정하면 된다.

Hann Window는 sinewave와 같은 파형을 잘 뽑아내는 성질을 가지는 것으로 알려져있다.

(이거보다 좀 더 성능 좋은 window는 Blackman-Harris Window)

 

...
w = hann(N);			% Hann Window
X = fft(x.*w)/N;        % Normalize the Windowed FFT Data
...

window를 적용하지 않은 파란색 그래프에 비해 window를 적용한 주황색 그래프의 noise floor가 훨씬 좋아진 것을 알 수 있다.

여기서 의문점이 들 수 있는데, amplitude가 2에서 오히려 더 멀어졌다는 점이다.

하지만 Hann Window의 성분은 fundamental 주위의 2개의 점을 추가로 더 포함해야 정확하게 계산이 가능하다.

위키피디아에도 잘 설명이 되어 있으니 참고. (3개의 non-zero value라고 설명되어 있음)

 

Hann function - Wikipedia

From Wikipedia, the free encyclopedia Jump to navigation Jump to search Hann function (left), and its frequency response (right) The Hann function, used to perform Hann smoothing,[1] is named after the Austrian meteorologist Julius von Hann, is a window fu

en.wikipedia.org

ΔΣ ADC의 교과서인 아래 책의 Appendix (Spectral Estimation)에도 관련 내용이 잘 나와있으니 참고.

S. Pavan, R. Schreier, G. C. Temes, "Understanding Delta-Sigma Data Converters"

 

>>figure(2)
>>stem(x_freq, abs(y_FFTout))
>>xlim([3.5 4.5])

 

4Hz를 중심으로 3개의 인접한 point가 있는 것을 알 수 있고,

0.950349 + 0.707285 + 0.300212 = 1.957846

이므로 우리가 설정했던 amplitude인 2에 상당히 근접해졌다는 것도 알 수 있다.

 

 

 

 

 

 

4. Coherent Sampling을 적용해 FFT 돌리기

뭔가 꽤 많이 나아진 것 같기는 한데, 위 결과로도 만족이 되지 않는다.

skirt가 생기는 것이 영 거슬려서 좀 더 깔끔하게 만들어보고 싶다면 coherent sampling을 이용하면 된다.

Maxim Integrated 사의 자료에 자세한 설명이 나와있으니 참고.

 

 

Coherent Sampling vs. Window Sampling | Maxim Integrated

Abstract: One of the most useful techniques for evaluating the dynamic performance of fast and ultra-fast data converters is coherent sampling. This technique increases the spectral resolution of a Fast Fourier Transform (FFT) and eliminates the need for w

www.maximintegrated.com

 

$$M_{cycle} = N_{samp} \times {f_{in} \over F_{s}}$$

 

이므로 $M_{cycle}$은 327.68이 되어야 한다. 여기에 가장 가까운 소수 (prime number)는 331이므로,

이 값을 넣어 위 식에다 다시 $f_{in}$을 역산해보면 4.0405가 나온다.

(혹은 $F_{s}$를 역산해서 MATLAB code에 넣고 돌려도 된다.)

 

 

위와 같이 조금 더 나은 모양의 fft가 나온다. 점점 더 좋아지고 있다.

 

 

Linear X-Axis의 데이터를 뽑아보면 위와 같다. Coherent sampling을 하지 않은 경우에 비해

딱 3개의 point만 도드라지게 나타나고 나머지 녀석들은 좀 더 죽어버린 것을 알 수 있다.

3 point의 amplitude를 합쳐보면

0.998939 + 0.528711 + 0.471502 = 1.999152

로, 이 정도면 그냥 값이 2라고 봐도 무방하다.

 

 

 

 

 

 

5. Spectral Leakage를 조금 더 최소화 해서 FFT 돌리기

난 이 정도의 skirt도 보기 싫어 죽겠어 ! 라고 한다면 남은 방법이 하나 더 있다. Sample 개수를 늘리면 된다.

$N_{samp}$를 216인 65,536으로 잡고 hann window도 적용하고

coherent sampling도 적용해서 $f_{in}$을 3.9993286 ($M_{cycle}$이 2621일 때)로 잡으면 아래와 같이 fft 결과가 나온다.

 

 

Sample 개수가 213인 파란색 그래프에 비해 216인 주황색 그래프의 noise floor가 훨씬 더 좋아졌고,

skirt도 훨씬 더 적게 발생하는 것을 볼 수 있다.

 

 

 

 

 

 

6. Common Error

요청된 131072x131072(128.0GB) 배열이 최대 배열 크기 기본 설정(7.9GB)을 초과합니다.
이로 인해 MATLAB이 응답하지 않을 수 있습니다.

 

위와 같은 error가 뜬다면, fft 함수의 input vector가 column vector인지 체크해보아야 한다.

예를 들어 input vector A가 1x131072 double과 같은 row vector data라면,

아래 code를 이용해 column vector로 바꾼 뒤 사용을 하면 위 error가 뜨지 않는다.

 

% Transpose
A = A';

 

 

 

저작자표시 비영리 변경금지

    티스토리툴바