📂푸리에해석

줄리아에서 이산 푸리에 변환 행렬 구현하기

설명

$N$차원 벡터 $\mathbf{x}$의 이산 푸리에 변환은 아래와 같은 행렬 곱으로 나타난다.

$$ \widehat{\mathbf{x}} = F \mathbf{x} $$

이때 $F$는, $\omega = e^{-i2\pi/N}$이라 나타내면,

$$ F = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \cdots & \omega^{N-1} \\ 1 & \omega^2 & \omega^4 & \cdots & \omega^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{N-1} & \omega^{2(N-1)} & \cdots & \omega^{(N-1)^{2}} \end{bmatrix} $$

역변환은 $F^{-1} = \dfrac{1}{N} \overline{F}$로 나타난다. 이때 $\overline{F}$는 $F$의 켤레 복소 행렬이다.

구현

1D 이산푸리에변환

줄리아에서는 다음과 같이 구현할 수 있다.

function dft_matrix(n)
    ω = exp(-2π*im/n)
    F = [ω^(i*j) for i in 0:n-1, j in 0:n-1]
    
    return F
end

물론 다음과 같이 한 줄로 작성해도 된다.

function dft_matrix(n)
    return [exp(2π *im*i*j/n) for i in 0:n-1, j in 0:n-1]
end

고속푸리에변환 패키지인 FFTW을 불러와 확인해보면 다음과 같다.

using FFTW
using Plots

F = dft_matrix(128)
x = rand(128)

ξ = fft(x)
η = F*x

plot(
    plot(abs.(ξ), label="FFTW"),
    plot(abs.(η), label="DFT matrix"),
    plot(abs.(ξ .- η), label="Difference"),
    layout=(3,1)
)

2D 이산푸리에변환

2차원 변환을 구현하기 위해 이미지를 하나 불러오자.

using Tomography

p = phantom(128, 2)
heatmap(p, yflip=true, c=:viridis)

FFTW와 비교해보면 다음과 같다.

F = dft_matrix(256)
ξ = fft(p)
η = F*p*transpose(F)

plot(
    heatmap(abs.(ξ), yflip=true, title="FFTW"),
    heatmap(abs.(η), yflip=true, aspect_ratio=1, title="DFT matrix"),
    heatmap(abs.(ξ .- η), yflip=true, aspect_ratio=1, title="Difference"),
    layout=(1,3)
)

역변환

역변환은 변환 행렬의 켤레복소행렬이므로 다음과 같이 구현할 수 있다.

function dft_inverse(F)
    n = size(F, 1)
    return conj(F) / n
end

F = dft_matrix(128)
x = rand(128)
F⁻¹ = dft_inverse(F)
plot(
    plot(x, label="Original signal x"),
    plot(real(F⁻¹*F*x), label="Reconstructed signal F⁻¹*F*x"),
    plot(abs.(x .- F⁻¹*F*x), label="Difference"),
    layout = (3,1)
)

2차원의 경우도 마찬가지로,

p = phantom(256, 2)
F = dft_matrix(256)
F⁻¹ = dft_inverse(F)
p_reconstructed = F⁻¹*F*p*transpose(F)*transpose(F⁻¹)

plot(
    heatmap(p, yflip=true, title="Original image"),
    heatmap(real(p_reconstructed), yflip=true, title="Reconstructed image"),
    heatmap(abs.(p .- real(p_reconstructed)), yflip=true, title="Difference", clims=(0, 1e-10)),
    layout=(1,3)
)

fftshift

이산 푸리에 변환은 기본적으로 주파수가 $0$인 성분이 처음에 나온다. 푸리에변환과 관련된 패키지에서는 $0$ 주파수 성분이 중앙에 오도록 바꿔주는 함수가 존재하는데 (대개 fftshift라는 이름을 가진다) 이를 구현해보자. 앞의 절반 성분과 뒤의 절반 성분을 바꿔주면 되므로 다음과 같이 쉽게 구현할 수 있다.

Fs = 1000                    # 진동수
T = 1/1000                   # 샘플링 간격
L = 1000                     # 신호의 길이
x = [i for i in 0:L-1].*T    # 신호의 도메인

F = dft_matrix(L)        # 0 주파수 성분이 첫번째 성분에 오도록 하는 푸리에 변환 행렬
S = shifted_matrix(F)    # 0 주파수 성분이 중앙에 오도록 하는 푸리에 변환 행렬

y₁ = sin.(2π*100*x)            # 진동수가 100인 사인파
y₂ = 0.5sin.(2π*200*x)         # 진동수가 200인 사인파
y₃ = 2sin.(2π*250*x)           # 진동수가 250인 사인파
y = y₁ + y₂ + y₃ + randn(L)    # 잡음이 섞인 신호

ξ = Fs*[i for i in 0:L-1]/L    # 주파수 도메인

plot(
    plot(ξ       , abs.(F*y), label="DFT of the signal"        ),
    plot(ξ .- 500, abs.(S*y), label="Shifted DFT of the signal"),
    layout=(2,1)
)

2차원에서의 결과는 다음과 같다.

p = phantom(256, 2)
F = dft_matrix(256)
S = shifted_matrix(F)

ξ = F*p*transpose(F)
η = S*p*transpose(S)

plot(
    heatmap(abs.(ξ), yflip=true, clims=(0, 200)),
    heatmap(abs.(η), yflip=true, clims=(0, 200)),
    layout=(1,2),
)

코드 전문

function dft_matrix(n)
    ω = exp(-2π*im/n)
    F = [ω^(i*j) for i in 0:n-1, j in 0:n-1]
    
    return F
end


function dft_matrix(n)
    return [exp(2π *im*i*j/n) for i in 0:n-1, j in 0:n-1]
end


using FFTW
using Plots

F = dft_matrix(128)
x = rand(128)

ξ = fft(x)
η = F*x

plot(
    plot(abs.(ξ), label="FFTW"),
    plot(abs.(η), label="DFT matrix"),
    plot(abs.(ξ .- η), label="Difference"),
    layout=(3,1)
)


using Tomography

p = phantom(128, 2)
heatmap(p, yflip=true, c=:viridis)


F = dft_matrix(256)
ξ = fft(p)
η = F*p*transpose(F)

plot(
    heatmap(abs.(ξ), yflip=true, title="FFTW"),
    heatmap(abs.(η), yflip=true, aspect_ratio=1, title="DFT matrix"),
    heatmap(abs.(ξ .- η), yflip=true, aspect_ratio=1, title="Difference"),
    layout=(1,3)
)


function dft_inverse(F)
    n = size(F, 1)
    return conj(F) / n
end

F = dft_matrix(128)
x = rand(128)
F⁻¹ = dft_inverse(F)
plot(
    plot(x, label="Original signal x"),
    plot(real(F⁻¹*F*x), label="Reconstructed signal F⁻¹*F*x"),
    plot(abs.(x .- F⁻¹*F*x), label="Difference"),
    layout = (3,1)
)


p = phantom(256, 2)
F = dft_matrix(256)
F⁻¹ = dft_inverse(F)
p_reconstructed = F⁻¹*F*p*transpose(F)*transpose(F⁻¹)

plot(
    heatmap(p, yflip=true, title="Original image"),
    heatmap(real(p_reconstructed), yflip=true, title="Reconstructed image"),
    heatmap(abs.(p .- real(p_reconstructed)), yflip=true, title="Difference", clims=(0, 1e-10)),
    layout=(1,3)
)


Fs = 1000                    # 진동수
T = 1/1000                   # 샘플링 간격
L = 1000                     # 신호의 길이
x = [i for i in 0:L-1].*T    # 신호의 도메인

F = dft_matrix(L)        # 0 주파수 성분이 첫번째 성분에 오도록 하는 푸리에 변환 행렬
S = shifted_matrix(F)    # 0 주파수 성분이 중앙에 오도록 하는 푸리에 변환 행렬

y₁ = sin.(2π*100*x)            # 진동수가 100인 사인파
y₂ = 0.5sin.(2π*200*x)         # 진동수가 200인 사인파
y₃ = 2sin.(2π*250*x)           # 진동수가 250인 사인파
y = y₁ + y₂ + y₃ + randn(L)    # 잡음이 섞인 신호

ξ = Fs*[i for i in 0:L-1]/L    # 주파수 도메인

plot(
    plot(ξ       , abs.(F*y), label="DFT of the signal"        ),
    plot(ξ .- 500, abs.(S*y), label="Shifted DFT of the signal"),
    layout=(2,1)
)


p = phantom(256, 2)
F = dft_matrix(256)
S = shifted_matrix(F)

ξ = F*p*transpose(F)
η = S*p*transpose(S)

plot(
    heatmap(abs.(ξ), yflip=true, clims=(0, 200)),
    heatmap(abs.(η), yflip=true, clims=(0, 200)),
    layout=(1,2),
)

환경

• OS: Windows11
• Version: Julia 1.10.0, FFTW v1.8.0, Plots v1.40.3, Tomography v0.1.5