fft

고속 푸리에 변환(Fast Fourier Transform)

구문

Y = fft(X)

Y = fft(X,n)

Y = fft(X,n,dim)

설명

Y = fft(X)는 고속 푸리에 변환(FFT) 알고리즘을 사용하여 X의 이산 푸리에 변환(DFT)을 계산합니다. Y는 X와 크기가 동일합니다.

• X가 벡터인 경우 fft(X)는 벡터의 푸리에 변환을 반환합니다.

• X가 행렬인 경우 fft(X)는 X의 열을 벡터로 취급하고 각 열의 푸리에 변환을 반환합니다.

• X가 다차원 배열인 경우 fft(X)는 크기가 1이 아닌 첫 번째 배열 차원의 값을 벡터로 취급하고 각 벡터의 푸리에 변환을 반환합니다.

Y = fft(X,n)은 n포인트 DFT를 반환합니다.

• X가 벡터이고 X의 길이가 n보다 작은 경우 X는 n 길이까지 후행 0으로 채워집니다.

• X가 벡터이고 X의 길이가 n보다 큰 경우 X는 n 길이로 잘립니다.

• X가 행렬인 경우 각 열은 벡터에서 취급되는 것처럼 취급됩니다.

• X가 다차원 배열인 경우 크기가 1이 아닌 첫 번째 배열 차원은 벡터에서 취급되는 것처럼 취급됩니다.

Y = fft(X,n,dim)은 차원 dim을 따라 푸리에 변환을 반환합니다. 예를 들어, X가 행렬인 경우 fft(X,n,2)는 각 행의 n포인트 푸리에 변환을 반환합니다.

예제

잡음이 있는 신호

푸리에 변환을 사용하여 잡음에 묻힌 신호의 주파수 성분을 찾고 피크 주파수의 진폭을 구합니다.

샘플링 주파수 1kHz와 신호 지속 시간 1.5초로 신호의 파라미터를 지정합니다.

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector

진폭 0.8의 DC 오프셋, 진폭 0.7의 50Hz 정현파, 진폭 1의 120Hz 정현파를 포함하는 신호를 생성합니다.

S = 0.8 + 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

분산 4, 평균 0의 랜덤 잡음으로 손상시킵니다.

X = S + 2*randn(size(t));

잡음이 있는 신호를 시간 영역에 플로팅합니다. 주파수 성분이 플롯에서 시각적으로 뚜렷하게 보이지 않습니다.

plot(1000*t,X)
title("Signal Corrupted with Zero-Mean Random Noise")
xlabel("t (milliseconds)")
ylabel("X(t)")

신호의 푸리에 변환을 계산합니다.

Y = fft(X);

푸리에 변환에는 복소수가 포함되므로 fft 스펙트럼의 복소수 크기를 플로팅합니다.

plot(Fs/L*(0:L-1),abs(Y),"LineWidth",3)
title("Complex Magnitude of fft Spectrum")
xlabel("f (Hz)")
ylabel("|fft(X)|")

이 플롯은 DC 오프셋에 대한 0Hz에서의 피크를 포함하여 5개의 주파수 피크를 보여줍니다. 이 예제에서 신호는 0Hz, 50Hz, 120Hz에서 세 개의 주파수 피크를 가질 것으로 예상됩니다. 여기서 플롯의 두 번째 절반은 0Hz에서의 피크를 제외하고 첫 번째 절반이 거울 반사된 형태입니다. 그 이유는 시간 영역 신호의 이산 푸리에 변환이 주기적 특성을 갖고 있기 때문입니다. 여기서 스펙트럼의 첫 번째 절반은 양의 주파수이고 두 번째 절반은 음의 주파수이며, 이 때 첫 번째 요소는 영주파수에 사용하도록 예약됩니다.

실수 신호의 경우 fft 스펙트럼은 양의 주파수 스펙트럼이 음의 주파수 스펙트럼의 켤레 복소수인 양방향 스펙트럼입니다. 양의 주파수 및 음의 주파수 fft 스펙트럼을 표시하려면 fftshift를 사용하면 됩니다. L의 길이가 짝수인 경우 주파수 영역은 나이퀴스트 주파수 -Fs/2의 음수부터 Fs/2-Fs/L까지이며 Fs/L의 간격 또는 주파수 분해능을 갖습니다.

plot(Fs/L*(-L/2:L/2-1),abs(fftshift(Y)),"LineWidth",3)
title("fft Spectrum in the Positive and Negative Frequencies")
xlabel("f (Hz)")
ylabel("|fft(X)|")

세 개의 주파수 피크의 진폭을 구하려면 Y의 fft 스펙트럼을 단방향 진폭 스펙트럼으로 변환합니다. fft 함수는 원래 신호와 변환된 신호 간의 스케일링 인자 L을 포함하므로 Y를 L로 나누어 다시 스케일링합니다. fft 스펙트럼의 복소수 크기를 구합니다. 양의 주파수 스펙트럼이 음의 주파수 스펙트럼의 켤레 복소수인 양방향 진폭 스펙트럼 P2의 피크 진폭은 시간 영역 신호의 피크 진폭의 절반입니다. 단방향 스펙트럼으로 변환하려면 양방향 스펙트럼 P2의 첫 번째 절반을 취합니다. 양의 주파수 스펙트럼에 2를 곱합니다. P1(1)과 P1(end)에 2를 곱할 필요는 없습니다. 이들 진폭은 각각 영주파수와 나이퀴스트 주파수에 해당하고 음의 주파수에 켤레 복소수 쌍을 갖지 않기 때문입니다.

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

단방향 스펙트럼의 주파수 영역 f를 정의합니다. 단방향 진폭 스펙트럼 P1을 플로팅합니다. 예상대로, 진폭이 0.8, 0.7, 1에 가깝지만 잡음이 추가되었기 때문에 정확하지는 않습니다. 대부분의 경우, 신호가 길수록 더 나은 주파수 근삿값을 얻을 수 있습니다.

f = Fs/L*(0:(L/2));
plot(f,P1,"LineWidth",3) 
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

이제, 손상되지 않은 원래 신호의 푸리에 변환을 사용하여 정확한 진폭 0.8, 0.7, 1.0을 가져옵니다.

Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

plot(f,P1,"LineWidth",3) 
title("Single-Sided Amplitude Spectrum of S(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")

가우스 펄스(Gaussian Pulse)

가우스 펄스를 시간 영역에서 주파수 영역으로 변환합니다.

샘플링 주파수 44.1kHz와 신호 지속 시간 1ms로 신호의 파라미터를 지정합니다. 표준편차가 0.1ms인 가우스 펄스(Gaussian Pulse)를 생성합니다.

Fs = 44100;         % Sampling frequency
T = 1/Fs;           % Sampling period
t = -0.5:T:0.5;     % Time vector
L = length(t);      % Signal length

X = 1/(0.4*sqrt(2*pi))*(exp(-t.^2/(2*(0.1*1e-3)^2)));

펄스를 시간 영역에 플로팅합니다.

plot(t,X)
title("Gaussian Pulse in Time Domain")
xlabel("Time (t)")
ylabel("X(t)")
axis([-1e-3 1e-3 0 1.1]) 

fft의 실행 시간은 변환의 길이에 따라 다릅니다. 작은 소인수만으로 이루어진 변환 길이는 큰 소인수로 이루어진 변환의 길이보다 훨씬 빠르게 실행됩니다.

이 예제에서 신호 길이 L은 아주 큰 소인수인 44,101입니다. fft 성능을 개선하려면 원래 신호 길이에서 2의 다음 거듭제곱이 되는 입력 길이를 식별하십시오. 이 입력 길이의 fft를 호출하면 펄스 X를 지정된 변환 길이까지 후행 0으로 채웁니다.

n = 2^nextpow2(L);

가우스 펄스를 주파수 영역으로 변환합니다.

Y = fft(X,n);

주파수 영역을 정의하고 고유한 주파수를 플로팅합니다.

f = Fs*(0:(n/2))/n;
P = abs(Y/sqrt(n)).^2;

plot(f,P(1:n/2+1)) 
title("Gaussian Pulse in Frequency Domain")
xlabel("f (Hz)")
ylabel("|P(f)|")

코사인파

시간 영역과 주파수 영역의 코사인파를 비교합니다.

샘플링 주파수 1kHz와 신호 지속 시간 1초로 신호의 파라미터를 지정합니다.

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

각 행이 스케일링된 주파수의 코사인파를 나타내는 행렬을 만듭니다. 결과값 X는 3×1000 행렬이 됩니다. 첫 번째 행의 파동 주파수는 50이고, 두 번째 행의 파동 주파수는 150이고, 세 번째 행의 파동 주파수는 300입니다.

x1 = cos(2*pi*50*t);          % First row wave
x2 = cos(2*pi*150*t);         % Second row wave
x3 = cos(2*pi*300*t);         % Third row wave

X = [x1; x2; x3];

단일 Figure에 있는 X의 각 행에서 처음 100개 요소를 순서대로 플로팅하고 주파수를 비교합니다.

for i = 1:3
    subplot(3,1,i)
    plot(t(1:100),X(i,1:100))
    title("Row " + num2str(i) + " in the Time Domain")
end

X의 행을 따라, 즉 각 신호에 대해 fft를 사용하도록 dim 인수를 지정합니다.

dim = 2;

신호의 푸리에 변환을 계산합니다.

Y = fft(X,L,dim);

각 신호의 양방향 스펙트럼과 단방향 스펙트럼을 계산합니다.

P2 = abs(Y/L);
P1 = P2(:,1:L/2+1);
P1(:,2:end-1) = 2*P1(:,2:end-1);

주파수 영역에서 단일 Figure의 각 행에 대해 단방향 진폭 스펙트럼을 플로팅합니다.

for i=1:3
    subplot(3,1,i)
    plot(0:(Fs/L):(Fs/2-Fs/L),P1(i,1:L/2))
    title("Row " + num2str(i) + " in the Frequency Domain")
end

정현파(Sinusoidal)의 위상

주파수 15Hz와 40Hz의 정현파 2개로 구성된 신호를 만듭니다. 첫 번째 정현파는 위상이 $-\pi /4$인 코사인파이고, 두 번째 정현파는 위상이 $\pi /2$ 인 코사인파입니다. 신호를 100Hz에서 1초 동안 샘플링합니다.

Fs = 100;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*15*t - pi/4) + cos(2*pi*40*t + pi/2);

신호의 푸리에 변환을 계산합니다. 변환의 크기를 주파수의 함수로 플로팅합니다.

y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*Fs;
stem(f,abs(z))
title("Double-Sided Amplitude Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("|y|")
grid

크기가 작은 변환 값을 제거하여 변환 위상을 계산합니다. 위상을 주파수의 함수로 플로팅합니다.

tol = 1e-6;
z(abs(z) < tol) = 0;

theta = angle(z);

stem(f,theta/pi)
title("Phase Spectrum of x(t)")
xlabel("Frequency (Hz)")
ylabel("Phase/\pi")
grid

FFT의 보간

신호의 푸리에 변환을 0으로 채워 보간합니다.

샘플링 주파수 80kHz와 신호 지속 시간 0.8초로 신호의 파라미터를 지정합니다.

Fs = 80;
T = 1/Fs;
L = 65;
t = (0:L-1)*T;

2Hz 정현파 신호와 더 높은 고조파의 중첩을 만듭니다. 신호에는 2Hz 코사인파, 4Hz 코사인파, 6Hz 사인파가 포함됩니다.

X = 3*cos(2*pi*2*t) + 2*cos(2*pi*4*t) + sin(2*pi*6*t);

신호를 시간 영역에 플로팅합니다.

plot(t,X)
title("Signal superposition in time domain")
xlabel("t (ms)")
ylabel("X(t)")

신호의 푸리에 변환을 계산합니다.

Y = fft(X);

신호의 단방향 진폭 스펙트럼을 계산합니다.

f = Fs*(0:(L-1)/2)/L;
P2 = abs(Y/L);
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);

주파수 영역에서 단방향 스펙트럼을 플로팅합니다. 신호의 시간 추출이 매우 짧기 때문에 푸리에 변환의 주파수 분해능은 4Hz 근처의 피크 주파수를 나타낼 만큼 정밀하지 않습니다.

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Original Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

피크 주파수를 더 정확히 평가하려면 원래 신호를 0으로 채워서 분석 윈도우의 길이를 늘릴 수 있습니다. 이 메서드는 보다 정확한 주파수 분해능으로 신호의 푸리에 변환을 자동 보간합니다.

원래 신호 길이 다음의 2의 거듭제곱이 되는 새 입력 길이를 식별하십시오. 신호 X를 후행 0으로 채워 길이를 연장합니다. 0으로 채워진 신호의 푸리에 변환을 계산합니다.

n = 2^nextpow2(L);
Y = fft(X,n);

채워진 신호의 단방향 진폭 스펙트럼을 계산합니다. 신호 길이 n이 65에서 128로 늘어났기 때문에 주파수 분해능은 Fs/n(0.625Hz)가 됩니다.

f = Fs*(0:(n/2))/n;
P2 = abs(Y/L);
P1 = P2(1:n/2+1);
P1(2:end-1) = 2*P1(2:end-1);

채워진 신호의 단방향 스펙트럼을 플로팅합니다. 이 새로운 스펙트럼은 0.625Hz의 주파수 분해능 내에서 2Hz, 4Hz, 6Hz 근처의 피크 주파수를 보여줍니다.

plot(f,P1,"-o") 
title("Single-Sided Spectrum of Padded Signal")
xlabel("f (Hz)")
ylabel("|P1(f)|")

입력 인수

X — 입력 배열

벡터 | 행렬 | 다차원 배열

입력 배열로, 벡터, 행렬, 다차원 배열 중 하나로 지정됩니다.

X가 빈 0×0 행렬인 경우 fft(X)는 빈 0×0 행렬을 반환합니다.

데이터형: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical
복소수 지원 여부: 예

n — 변환 길이

[] (디폴트 값) | 음이 아닌 정수 스칼라

변환 길이로, [] 또는 음이 아닌 정수 스칼라로 지정됩니다. 양의 정수 스칼라로 변환 길이를 지정하면 fft의 성능을 향상시킬 수 있습니다. 이 길이는 일반적으로 2의 거듭제곱으로 지정되거나 작은 소인수의 곱(7보다 크지 않은 소인수)으로 지정됩니다. n이 신호의 길이보다 작은 경우 fft는 n 번째 요소 이후의 나머지 신호 값을 무시하고 잘린 결과를 반환합니다. n이 0인 경우 fft는 빈 행렬을 반환합니다.

예: n = 2^nextpow2(size(X,1))

데이터형: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

dim — 연산을 수행할 차원

양의 정수 스칼라

연산을 수행할 차원으로, 양의 정수 스칼라로 지정됩니다. 차원을 지정하지 않을 경우 디폴트 값은 크기가 1보다 큰 첫 번째 배열 차원이 됩니다.

• fft(X,[],1)은 X의 열을 따라 동작하고 각 열의 푸리에 변환을 반환합니다.

  

• fft(X,[],2)는 X의 행을 따라 동작하고 각 행의 푸리에 변환을 반환합니다.

  

dim이 ndims(X)보다 큰 경우 fft(X,[],dim)은 X를 반환합니다. n이 지정된 경우 fft(X,n,dim)은 차원 dim을 따라 X를 길이 n까지 채우거나 자릅니다.

데이터형: double | single | int8 | int16 | int32 | uint8 | uint16 | uint32 | logical

출력 인수

Y — 주파수 영역 표현

벡터 | 행렬 | 다차원 배열

주파수 영역 표현으로, 벡터, 행렬, 다차원 배열 중 하나로 반환됩니다.

X가 single형인 경우 fft는 기본적으로 단정밀도로 계산되며 Y 또한 single형이 됩니다. 그렇지 않으면, Y가 double형으로 반환됩니다.

Y의 크기는 다음과 같습니다.

• Y = fft(X) 또는 Y = fft(X,[],dim)의 경우 Y와 X의 크기는 같습니다.

• Y = fft(X,n,dim)의 경우 size(Y,dim)의 값은 n과 같고 다른 모든 차원의 크기는 X 차원의 크기 그대로 유지됩니다.

X가 실수인 경우 Y는 켤레 대칭이고 Y에 있는 고유한 점 개수는 ceil((n+1)/2)입니다.

데이터형: double | single

세부 정보

벡터의 이산 푸리에 변환

Y = fft(X)는 푸리에 변환을 구현하고 X = ifft(Y)는 푸리에 역변환을 구현합니다. 이 변환은 길이 n의 X와 Y에 대해 다음과 같이 정의됩니다.

$$\begin{array}{l}Y(k)={\displaystyle \sum _{j=1}^{n}X}(j)W_n^{(j-1)\text{}(k-1)}\\ X(j)=\frac{1}{n}{\displaystyle \sum _{k=1}^{n}Y}(k)W_n{}^{-(j-1)\text{}(k-1)},\end{array}$$

여기서

$$W_n=e^{(-2\pi i)/n}$$

은 n개의 단위근 중 하나입니다.

팁

• fft의 실행 시간은 변환의 길이에 따라 다릅니다. 작은 소인수만(7보다 크지 않음)으로 이루어진 변환의 길이는 소수이거나 큰 소인수로 이루어진 변환의 길이보다 훨씬 빠르게 실행됩니다.

• 대부분의 n 값의 경우 실수 입력 DFT는 복소수 입력 DFT의 계산 시간보다 약 2배 더 빠릅니다. 그러나 n이 큰 소인수인 경우 속도 차이는 약간 있거나 아예 없을 수 있습니다.

• 유틸리티 함수 fftw를 사용하여 fft의 속도를 높일 수 있습니다. 이 함수는 특정 크기와 차원의 FFT를 계산하는 데 사용되는 알고리즘의 최적화를 제어합니다.

알고리즘

FFT 함수(fft, fft2, fftn, ifft, ifft2, ifftn)는 FFTW ,[2] 라이브러리를 기반으로 합니다.

참고 문헌

[1] FFTW (https://www.fftw.org)

[2] Frigo, M., and S. G. Johnson. “FFTW: An Adaptive Software Architecture for the FFT.” Proceedings of the International Conference on Acoustics, Speech, and Signal Processing. Vol. 3, 1998, pp. 1381-1384.

확장 기능

C/C++ 코드 생성

MATLAB® Coder™를 사용하여 C 코드나 C++ 코드를 생성할 수 있습니다.

사용법 관련 참고 및 제한 사항:

• 가변 크기 데이터와 관련된 제한 사항을 보려면 Variable-Sizing Restrictions for Code Generation of Toolbox Functions (MATLAB Coder) 항목을 참조하십시오.

• MEX 출력의 경우, MATLAB^® Coder™는 MATLAB이 FFT 알고리즘에 사용하는 라이브러리를 사용합니다. 독립 실행형 C/C++ 코드의 경우, 기본적으로 코드 생성기는 FFT 라이브러리 호출을 생성하는 대신 FFT 알고리즘에 대한 코드를 생성합니다. 설치된 특정 FFTW 라이브러리에 대한 호출을 생성하려면 FFT 라이브러리 콜백 클래스를 제공하십시오. FFT 라이브러리 콜백 클래스에 대한 자세한 내용은 coder.fftw.StandaloneFFTW3Interface (MATLAB Coder)를 참조하십시오.

• 시뮬레이션 소프트웨어는 MATLAB이 FFT 알고리즘에 사용하는 라이브러리를 사용하여 MATLAB Function 블록을 시뮬레이션합니다. C/C++ 코드 생성 시, 기본적으로 코드 생성기는 FFT 라이브러리 호출을 생성하는 대신 FFT 알고리즘에 대한 코드를 생성합니다. 설치된 특정 FFTW 라이브러리에 대한 호출을 생성하려면 FFT 라이브러리 콜백 클래스를 제공하십시오. FFT 라이브러리 콜백 클래스에 대한 자세한 내용은 coder.fftw.StandaloneFFTW3Interface (MATLAB Coder)를 참조하십시오.

• CRL(Code Replacement Library)을 사용하여, Neon 확장이 적용된 ARM^® Cortex^®-A 프로세서에서 실행되는 최적화된 코드를 생성할 수 있습니다. 이 최적화된 코드를 생성하려면 Embedded Coder^® Support Package for ARM Cortex-A Processors (Embedded Coder)를 설치해야 합니다. ARM Cortex-A에 대해 생성된 코드는 Ne10 라이브러리를 사용합니다. 자세한 내용은 Ne10 Conditions for MATLAB Functions to Support ARM Cortex-A Processors (Embedded Coder)를 참조하십시오.

• CRL(Code Replacement Library)을 사용하여, ARM Cortex-M 프로세서에서 실행되는 최적화된 코드를 생성할 수 있습니다. 이 최적화된 코드를 생성하려면 Embedded Coder Support Package for ARM Cortex-M Processors (Embedded Coder)를 설치해야 합니다. ARM Cortex-M에 대해 생성된 코드는 CMSIS 라이브러리를 사용합니다. 자세한 내용은 CMSIS Conditions for MATLAB Functions to Support ARM Cortex-M Processors (Embedded Coder)를 참조하십시오.

GPU 코드 생성
GPU Coder™를 사용하여 NVIDIA® GPU용 CUDA® 코드를 생성할 수 있습니다.

스레드 기반 환경

MATLAB®의 backgroundPool을 사용해 백그라운드에서 코드를 실행하거나 Parallel Computing Toolbox™의 ThreadPool을 사용해 코드 실행 속도를 높일 수 있습니다.

이 함수는 스레드 기반 환경을 완전히 지원합니다. 자세한 내용은 스레드 기반 환경에서 MATLAB 함수 실행하기 항목을 참조하십시오.

GPU 배열

Parallel Computing Toolbox™를 사용해 GPU(그래픽스 처리 장치)에서 실행하여 코드 실행 속도를 높일 수 있습니다.

fft 함수는 GPU 배열을 부분적으로 지원합니다. 입력 데이터를 gpuArray (Parallel Computing Toolbox)로 지정하면 이 함수의 일부 구문이 GPU에서 실행됩니다. 사용법 관련 참고 및 제한 사항:

• 출력값 Y는 허수부가 모두 0인 경우에도 항상 복소수입니다.

자세한 내용은 GPU에서 MATLAB 함수 실행하기 (Parallel Computing Toolbox) 항목을 참조하십시오.

분산 배열

Parallel Computing Toolbox™를 사용하여 대규모 배열을 클러스터의 결합된 메모리에 걸쳐 분할할 수 있습니다.

사용법 관련 참고 및 제한 사항:

• 분산 배열에 대해 fft는 병렬 FFT 알고리즘을 사용하는 대신 소수 길이 FFT를 수행하기 위해 단일 워커에서 벡터를 수집합니다. 큰 소수 길이 벡터에 대한 FFT에서는 메모리 부족 오류가 발생할 수 있습니다.

자세한 내용은 분산 배열을 사용하여 MATLAB 함수 실행 (Parallel Computing Toolbox) 항목을 참조하십시오.

버전 내역

R2006a 이전에 개발됨

참고 항목

함수

• fft2 | fftn | fftw | fftshift | ifft | interpft

도움말 항목

• 푸리에 변환(Fourier Transform)

외부 웹사이트

• Fourier Analysis (MathWorks Teaching Resources)
• Convolution in Digital Signal Processing (MathWorks Teaching Resources)