CODEWITHSUNDEEP
pythonnumpy
Creator
Creator

Reputation: 33

Why I get different result for inbuilt and defined FFT in python?

I have the code below for fft2 performed by numpy and a 2d fft performed by direct code. an anyone point out why they are different? My inputmatreix is rA.

def DFT_matrix(N):
    i, j = np.meshgrid(np.arange(N), np.arange(N))
    omega = np.exp( - 2 * math.pi * 1J / N )
    W = np.power( omega, i * j ) / np.sqrt(N)
    return W

sizeM=40
rA=np.random.rand(sizeM,sizeM)
rAfft=np.fft.fft2(rA)

rAfftabs=np.abs(rAfft)+1e-9


dftMtx=DFT_matrix(sizeM)
dftR=dftMtx.conj().T
mA=dftMtx*rA*dftR

mAabs=np.abs(mA)+1e-9

print(np.allclose(mAabs, rAfftabs))

Upvotes: 1

Views: 344

Answers (1)

mandulaj
mandulaj

Reputation: 763

There are a few problems with your implementation.

1. DFT Matrix formula

First of all, as explained here, the formula for computing the DFT X of a MxN signal x is defined as:

Since you are computing the DFT for a MxM input, you just need to compute the DFT_Matrix once. Also note that due to the way W is defined, there is no need for conjugation and since W is symmetric and unitary there is no need for any transpose.

2. Matrix multiplication

When it comes to actually multiplying the matrixes together, you have to make sure to use the matrix multiplication operator @ instead of the element wise multiplier *

3. DFT_matrix normalization

By default the fft functions don't normalize your output. This means that before comparing the two outputs, you either have to divide the np.fft.fft2 result by sqrt(M*M) = M or you drop the np.sqrt(N) in your DFT_matrix function.

Summary:

Her is your example with the appropriate fixes for a MxN input. At the end, the magnitudes and angles are compared.

import numpy as np

def DFT_matrix(N):
    i, j = np.meshgrid(np.arange(N), np.arange(N))
    omega = np.exp( - 2 * np.pi * 1j / N )
    W = np.power( omega, i * j ) # Normalization by sqrt(N) Not included
    return W

sizeM=40
sizeN=20
np.random.seed(0)
rA=np.random.rand(sizeM,sizeN)

rAfft=np.fft.fft2(rA)


dftMtxM=DFT_matrix(sizeM)
dftMtxN=DFT_matrix(sizeN)

# Matrix multiply the 3 matrices together 
mA = dftMtxM @ rA @ dftMtxN


print(np.allclose(np.abs(mA), np.abs(rAfft)))
print(np.allclose(np.angle(mA), np.angle(rAfft)))

Both checks should evaluate to True. However note that the complexity of this algorithm, assuming M=N is N³ while the library's fft2 brings that down to N²log(N)!

Upvotes: 2

Related Questions