줄리아에서 고속 푸리에 변환(FFT) 사용하는 방법
개요 1 2
The Fastest Fourier Transform in the West(FFTW)는 매사추세츠 공과대학(MIT)에서 Matteo Frigo와 Steven G. Johnson에 의해 개발된 이산 푸리에 변환을 계산하기 위한 소프트웨어 라이브러리이다. FFT 구현을 위한 줄리아 프레임워크로 AbstractFFTs.jl라는 패키지가 존재하지만, 이는 이 자체로 쓰기 위해 존재하는 것이 아니라 FFTW.jl과 같이 고속 푸리에 변환을 구현할 때 도움을 주기 위해 만들어진 것이다.
This package is mainly not intended to be used directly. Instead, developers of packages that implement FFTs (such as FFTW.jl or FastTransforms.jl) extend the types/functions defined in AbstractFFTs. This allows multiple FFT packages to co-exist with the same underlying fft(x) and plan_fft(x) interface.3
요약
- 푸리에 변환:
fft()- $2$차원 배열에서 열 단위로 변환:
fft( ,[1]) - $2$차원 배열에서 행 단위로 변환:
fft( ,[2]) - 배열의 특정한 차원에 대해서만 변환:
fft( ,[n₁, n₂,...])
- $2$차원 배열에서 열 단위로 변환:
- 푸리에 역변환:
ifft() - $0$ 주파수 중앙 정렬:
fftshift()- 역변환:
ifftshift()
- 역변환:
- 주파수 샘플링:
fftfreq(n, fs=1)
코드
푸리에 변환
푸리에 변환의 표기로는 $\mathcal{F}[f]$, $\hat{f}$가 주요하게 쓰이는데, 줄리아에서는 이 표기를 코드에서 그대로 쓸 수 있다. 진동수가 $100$, $200$, $350$인 사인파를 $1/1000$ 간격으로 샘플링하고 이를 더하자.
using FFTW
using Plots
using LaTeXStrings
Fs = 1000 #진동수
T = 1/1000 #샘플링 간격
L = 1000 #신호의 길이
x = [i for i in 0:L-1].*T #신호의 도메인
f₁ = sin.(2π*100*x) #진동수가 100인 사인파
f₂ = 0.5sin.(2π*200*x) #진동수가 100인 사인파
f₃ = 2sin.(2π*350*x) #진동수가 100인 사인파
f = f₁ + f₂ + f₃
정의에 따라 $f$의 푸리에 변환 $\mathcal{F}f$는 $50$, $100$, $200$에서만 값을 갖는다. 또한 이산 푸리에변환의 정의에 의해 $y$축 대칭인 값을 얻으나, 기본적으로는 주파수가 $0$일 때의 값이 첫번째 값으로 되어있다. 따라서 신호의 주파수와 진폭을 확인하기 위함이라면 앞의 절반만 plot해도 무관하다.
Fs = 1000 # 샘플링 주파수
ℱf = fft(f) # 푸리에 변환
ξ = Fs*[i for i in 0:L/2-1]/L #주파수 도메인(절반)
plot(ξ, abs.(ℱf[1:Int(L/2)])*2/L, title=L"Fourier transform of $f$", label="")
xlabel!("frequency")
ylabel!("amplitude")
savefig("fft.png")
$0$ 주파수 중앙 정렬
푸리에 변환의 출력은 기본적으로 주파수가 $0$일 때의 값이 첫번째 값으로 오게끔 되어있다. 주파수가 $0$일 때의 값이 중앙에 오도록 하려면 fftshift()를 쓰면 된다. 이를 다시 되돌릴 때는 ifftshift()를 쓰면 된다. 참고로 ifftshift는 ifft + shift가 아니라, i(nverse) + fftshift이다. 즉 fftshift()의 역변환이므로 헷갈리지 말자.
p1 = plot(ξ, abs.(ℱf), title=L"$\mathcal{F}f$", label="", yticks=[], xticks=[0, 100, 200, 350, 500, 1000])
p2 = plot(ξ.-500, abs.(fftshift(ℱf)), title=L"$\mathrm{fftshift}(\mathcal{F}f)$", label="", yticks=[], xticks=[-500,-350,-200,-100,0,100,200,350,500])
plot(p1, p2, size=(800,400))
savefig("fftshift.png")
고차원 푸리에 변환
2차원 푸리에 변환의 값과 비교하기 위해 $x = [1\ 2\ 3\ 4]^{T}$의 푸리에 변환 값을 먼저 계산해두자.
julia> x = [1.0; 2; 3; 4]
4-element Vector{Float64}:
1.0
2.0
3.0
4.0
julia> fft(x)
4-element Vector{ComplexF64}:
10.0 + 0.0im
-2.0 + 2.0im
-2.0 + 0.0im
-2.0 - 2.0im
fft()는 2차원 배열을 입력으로 받으면 알아서 2차원 푸리에 변환을 반환한다. 혹은 첫번째와 두번째 차원에서 변환을 계산한다는 의미에서 fft(, [1,2])도 같은 결과를 반환한다.
julia> y = [x x x x]
4×4 Matrix{Float64}:
1.0 1.0 1.0 1.0
2.0 2.0 2.0 2.0
3.0 3.0 3.0 3.0
4.0 4.0 4.0 4.0
julia> fft(y)
4×4 Matrix{ComplexF64}:
40.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
-8.0+8.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
-8.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
-8.0-8.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
julia> fft(y, [1,2])
4×4 Matrix{ComplexF64}:
40.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
-8.0+8.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
-8.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
-8.0-8.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
따라서 각 열마다 변환을 하고 싶으면 fft(, [1]), 각 행마다 변환을 하고 싶으면 fft(, [2])를 쓰면 된다.
julia> fft(y, [1])
4×4 Matrix{ComplexF64}:
10.0+0.0im 10.0+0.0im 10.0+0.0im 10.0+0.0im
-2.0+2.0im -2.0+2.0im -2.0+2.0im -2.0+2.0im
-2.0+0.0im -2.0+0.0im -2.0+0.0im -2.0+0.0im
-2.0-2.0im -2.0-2.0im -2.0-2.0im -2.0-2.0im
julia> fft(y, [2])
4×4 Matrix{ComplexF64}:
4.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
8.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
12.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
16.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
주파수 샘플링
fftfreq(n, fs=1)
길이가 $n$이고 간격이 $fs/n$인 주파수 도메인을 반환한다. 앞서 설명했듯이 푸리에 변환은 $0$ 주파수 값이 가장 처오기 때문에 fftfreq()로 샘플링한 주파수의 첫번째 값은 $0$이다. 또한 앞의 절반이 양수, 뒤의 절반이 음수이다. 따라서 fftshift()를 쓰면 인덱스에 따른 오른차순으로 정렬된다.
julia> fftfreq(4, 1)
4-element Frequencies{Float64}:
0.0
0.25
-0.5
-0.25
julia> fftfreq(5, 1)
5-element Frequencies{Float64}:
0.0
0.2
0.4
-0.4
-0.2
julia> fftshift(fftfreq(4, 1))
-0.5:0.25:0.25
환경
- OS: Windows11
- Version: Julia 1.8.2, FFTW 1.5.0