ジュリアで高速フーリエ変換（FFT）を使用する方法 📂ジュリア

ジュリアで高速フーリエ変換（FFT）を使用する方法

概要 ¹ ²

The Fastest Fourier Transform in the West(FFTW)は、マサチューセッツ工科大学(MIT)のMatteo FrigoとSteven G. Johnsonによって開発された、離散フーリエ変換を計算するためのソフトウェアライブラリです。FFTの実装用にAbstractFFTs.jlというパッケージがありますが、これは自体を直接使用するためではなく、FFTW.jlなどの高速フーリエ変換を実装するときに支援するために作られたものです。

このパッケージは主に直接使われることを意図していません。代わりに、FFTs（例えばFFTW.jlやFastTransforms.jl）を実装するパッケージの開発者がAbstractFFTsで定義された型/関数を拡張します。これにより、同じ基本的なfft(x)とplan_fft(x)インターフェイスを持つ複数のFFTパッケージが共存できます。³

概要

フーリエ変換: fft()
- $2$次元配列で列ごとに変換: fft( ,[1])
- $2$次元配列で行ごとに変換: fft( ,[2])
- 配列の特定の次元だけを変換: fft( ,[n₁, n₂, ...])
フーリエ逆変換: ifft()
$0$ 中央の周波数: fftshift()
- 逆変換: ifftshift()
周波数サンプリング: fftfreq(n, fs=1)

コード

フーリエ変換

Juliaでは、フーリエ変換の表記として$\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₃

フーリエ変換: fft()
フーリエ逆変換: ifft()

定義によると、$f$のフーリエ変換$\mathcal{F}f$は$50$, $100$, $200$でのみ非ゼロ値を持ちます。また、離散フーリエ変換の定義により、$y$軸の周りに対称な値を得ますが、基本的には周波数が$0$の値が最初の値になっています。したがって、信号の周波数と振幅を確認することが目的なら、前半だけをプロットしてもよいでしょう。

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 ▷eq10◁", label="") 
xlabel!("frequency")
ylabel!("amplitude")
savefig("fft.png")

$0$ 中央の周波数

フーリエ変換の出力は、基本的に周波数が$0$の値を最初に置きます。周波数が$0$の値を中央にしたい場合は、fftshift()を使います。これを元に戻す場合はifftshift()を使いますが、これはifft+shiftではなく、逆変換 + fftshift、つまりfftshift()の逆操作であるため、混乱しないようにしましょう。

p1 = plot(ξ, abs.(ℱf), title=L"▷eq11◁", label="", yticks=[], xticks=[0, 100, 200, 350, 500, 1000]) 
p2 = plot(ξ.-500, abs.(fftshift(ℱf)), title=L"▷eq22◁", 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