Странная фаза для гауссовского пучка 2D - PullRequest
Новая
       1

Странная фаза для гауссовского пучка 2D

0 голосов
/ 07 апреля 2019

У меня есть гауссов луч в 2D:

2D Gaussian

После выполнения fft2 и угла я получаю странные результаты:

magnitude of FFT

phase of FFT

def finite2D(x,y, N, M, a, hx):
    f = np.array([[0.0]*N]*N)
    for i in range(len(x)):
        for k in range(len(y)):
            f[i][k] = np.exp(-(x[i]*x[i] + y[k]*y[k]))

    D1 = fftpack.fft2(f)
    D2 = fftpack.fftshift(D1)

    b = N*N/(4*a*M)
    x = np.linspace(-b, b, N)
    y = np.linspace(-b, b, N)
    xx, yy = np.meshgrid(x, y)
    plt.imshow(np.abs(D2))
    plt.show()

    plt.imshow(np.angle(D2))
    plt.show(True)
    return D2, phas

a = 5
N = 128
M = 256
b = N*N/(4*a*M)
hx = 2*a/N
x = np.linspace(-a, a, N)
y = np.linspace(-a, a, N)
finite2D(x,y, N, M, a, hx)

Это должна быть фаза 0 или близка к 0. Почему это не так, и как мне это исправить?

/// Обновлено: 

def finite2D(x,y, N, M, a, hx):
    f = np.array([[0.0]*N]*N)
    for i in range(len(x)):
        for k in range(len(y)):
            f[i][k] = np.exp(-(x[i]*x[i] + y[k]*y[k]))

    f = fftpack.ifftshift(f)
    D1 = fftpack.fft2(f)
    D2 = fftpack.fftshift(D1)

    b = N*N/(4*a*M)
    x = np.linspace(-b, b, N)
    y = np.linspace(-b, b, N)
    xx, yy = np.meshgrid(x, y)
    plt.imshow(np.abs(D2))
    plt.show()

    plt.imshow(np.angle(D2))
    plt.show(True)
    return D2

a = 5
N = 128
M = 256
b = N*N/(4*a*M)
hx = 2*a/N
x = np.linspace(-a, a, N, endpoint=False)
y = np.linspace(-a, a, N, endpoint=False)
finite2D(x,y, N, M, a, hx)

Фаза: Phase

1 Ответ

1 голос
/ 07 апреля 2019

БПФ предполагает, что источник находится в верхнем левом углу изображения.Таким образом, вы вычисляете БПФ гауссова со смещением на половину размера изображения.Этот сдвиг приводит к высокочастотному фазовому сдвигу в частотной области.

Чтобы решить эту проблему, вам необходимо сместить источник вашего гауссовского сигнала в верхний левый угол изображения.ifftshift делает это: 

f = fftpack.ifftshift(f)
D1 = fftpack.fft2(f)
D2 = fftpack.fftshift(D1)

Обратите внимание, что там, где величина очень мала, фаза определяется ошибками округления, не ожидайте нулевую фазу там.

Обновленный результат выглядит хорошо, но в центральной области все еще остается очень маленький градиент.Это вызвано полупиксельным сдвигом гауссиана.Этот сдвиг задается определением координат x и y: 

N = 128
x = np.linspace(-a, a, N)
y = np.linspace(-a, a, N)

Для N четного размера сделайте 

x = np.linspace(-a, a, N, endpoint=False)
y = np.linspace(-a, a, N, endpoint=False)

таким, чтобыобразец где x==0.

...