CODEWITHSUNDEEP
pythonnumpyfftnormalizationnorm
sten
sten

Reputation: 7476

fft normalization vs the norm-option

i have to often normalize the result of circular-convolution, so I 'borrowed'-and-modfied the following routine :

def fft_normalize(x):# Normalize a vector x in complex domain.
    c = np.fft.rfft(x)
    # Look at real and image as if they were real
    ri = np.vstack([c.real, c.imag])
    # Normalize magnitude of each complex/real pair
    norm = np.linalg.norm(ri, axis=0)
    if np.any(norm==0): norm[norm == 0] = np.float64(1e-308) #!fixme
    ri= np.divide(ri,norm)
    c_proj = ri[0,:] + 1j * ri[1,:]
    rv = np.fft.irfft(c_proj, n=x.shape[-1])
    return rv

def fft_convolution(a, b):
    return np.fft.irfft(np.fft.rfft(a) * np.fft.rfft(b))

so that I do this :

fft_normalize(fft_convolution(a,b))

i see in the numpy docs there is a 'norm' option. Is this equivalent to what i'm doing ? And if yes, which option should I use ?

def fft_convolution2(a, b):
    return np.fft.irfft(np.fft.rfft(a) * np.fft.rfft(b), norm='ortho')

When I test it it behaves better when I do fft_normalize()

Second, i had to add this line, but it does not seem, right. any ideas ?

    if np.any(norm==0): norm[norm == 0] = np.float64(1e-308) #!fixme

As a side note, if you know !! numpy docs mentions that they promote float32 to float64 and that scipy.fftpack does not !! Would fftpack be faster ! because scipy says fftpack is obsolete and there is no info on the new scipy ?

@cris are you sugesting i do it this way :

def fft_normalize(x):# Normalize a vector x in complex domain.
    c = np.fft.rfft(x)
    ri = np.vstack([c.real, c.imag])
    norm = np.abs(c)
    if np.any(norm==0): norm[norm == 0] = MIN #!fixme
    ri= np.divide(ri,norm)
    c_proj = ri[0,:] + 1j * ri[1,:]
    rv = np.fft.irfft(c_proj, n=x.shape[-1])
    return rv

Upvotes: 0

Views: 908

Answers (1)

Cris Luengo
Cris Luengo

Reputation: 60504

The norm argument to the FFT functions in NumPy determine whether the transform result is multiplied by 1, 1/N or 1/sqrt(N), with N the number of samples in the array. Normally, the inverse transform is normalized by dividing by N, and the forward transform is not. Specifying “ortho” here causes both transforms to be normalized by 1/sqrt(2). Specifying “forward” causes only the forward transform to be normalized by 1/N.

These normalizations are very different from the one you apply, where you normalize each element in the frequency domain.

Upvotes: 2

Related Questions