2013-03-03 2 views
0

Я пытаюсь расширить fft-код, который отлично работает для 1D-массивов в python для изображений. На самом деле я знаю, что проблема связана с логикой в ​​расширении. Я мало знаю о FFT, и я должен отправлять задания для обработки изображений. Я буду благодарен за любые подсказки или решенияКак продлить 1D FFT-код для вычисления БПФ изображений (2D) в python?

Вот код, На самом деле, я пытаюсь создать модуль для FFT в python, и он уже отлично работает для 1D с помощью сайта rosetta Code.

from cmath import exp, pi 
from math import log, ceil 

def fft(f): 
    N = len(f) 
    if N <= 1: return f 
    even = fft(f[0::2]) 
    odd = fft(f[1::2]) 
    return [even[k] + exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)] + \ 
      [even[k] - exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)] 

def pad(f): 
    n = len(f) 
    N = 2 ** int(ceil(log(n, 2))) 
    F = f + [0] * (N - n) 
    return F, n 

def unpad(F, n): 
    return F[0 : n] 

def pad2(f): 
    m, n = len(f), len(f[0]) 
    M, N = 2 ** int(ceil(log(m, 2))), 2 ** int(ceil(log(n, 2))) 
    F = [ [0]*N for _ in xrange(M) ] 
    for i in range(0, m): 
     for j in range(0, n): 
      F[i][j] = f[i][j] 
    return F, m, n 

def fft1D(f): 
    Fu, n = pad(f) 
    return fft(Fu), n 

def fft2D(f): 
    F, m, n = pad2(f) 
    M, N = len(F), len(F[0]) 
    Fuv = [ [0]*N for _ in xrange(M) ] 
    for i in range(0, M): 
     Fxv = fft(F[i]) 
     for j in range(0, N): 
      Fuv[i][j] = (fft(Fxv))[j] 
    return Fuv, [m, n] 

Я назвал этот модуль тис код:

from FFT import * 
f= [0, 2, 3, 4] 
F = fft1D(f) 
print f, F 
X, s = fft2D([[1,2,1,1],[2,1,2,2],[0,1,1,0], [0,1,1,1]]) 
for i in range(0, len(X)): 
    print X[i] 

Это результат:

[0, 2, 3, 4] ([(9+0j), (-3+2j), (-3+0j), (-3-2j)], 4) 
[(4+0j), (4-2.4492935982947064e-16j), (4+0j), (8+2.4492935982947064e-16j)] 
[(8+0j), (8+2.4492935982947064e-16j), (8+0j), (4-2.4492935982947064e-16j)] 
[0j, -2.33486982377251e-16j, (4+0j), (4+2.33486982377251e-16j)] 
[0j, (4+0j), (4+0j), (4+0j)] 

Первый 1д прекрасно, как я проверить результат с выходом Matlab, но для 2-го один выход Matlab:

>> fft([1,2,1,1;2,1,2,2;0,1,1,0;0,1,1,1]) 

ans = 

    3.0000    5.0000    5.0000    4.0000   
    1.0000 - 2.0000i 1.0000     0 - 1.0000i 1.0000 - 1.0000i 
    -1.0000    1.0000   -1.0000   -2.0000   
    1.0000 + 2.0000i 1.0000     0 + 1.0000i 1.0000 + 1.0000i 

Результат отличается, что означает, что я делаю что-то неправильно в логике кода. Пожалуйста, помогите, не беспокоясь, поскольку я до сих пор не изучал БПФ, так что я не могу полностью понять математику, возможно, после того, как я ее изучил, я могу понять проблему.

+1

есть причина вы не хотите использовать 'numpy'? он уже имеет п-мерный код FFT – isedev

+0

Да, так как это назначение для обработки изображений, использование любой библиотеки не будет принято, если это абсолютно необходимо. поэтому я не использую numpy, scipy или PIL для этого ................ – bistaumanga

ответ

2

Ваш код немного трудно поддается наблюдению, но похоже, что вы делаете БПФ в одном направлении одновременно. Посмотрите на интеграл от FT, вы увидите, что объединения x и y независимы. То есть (извините, это обозначение ужасно, ' указывает на функцию в пространстве Фурье)

FT(f(x, y), x) -> f'(k, y) 
FT(f'(k, y), y) -> f''(k, w) 

Так что вы хотите сделать, это принять FFT каждой строки (чем сделать N 1D FFTs) и засунуть результаты в новый массив (который берет вас от f(x, y) -> f'(k, y)). Затем возьмите FFT каждого столбца этого результата массива (делает М 1D БПФ) и засунуть эти результаты в другой новый массив (который берет вас от f'(k, y) -> f''(k, w).

Смежные вопросы