[SciPy-User] fft convolutions

[Alex Flint]

I am trying to perform 2d convolutions between a large 2d array A and a bunch of small 2d arrays B1...Bn. My approach is roughly:

a = fft(A,size)

for b in bs:

   ans = ifft(fft(b,size)*A)

   slow = convolve2d(A, b, 'same')

However, as implemented above, ans is offset an inconsistent amount from the answer produced by convolve2d, presumably because convolve2d is treating b as if the origin is in the center whereas fft treats b as if the origin is at the top left (but it doesn't seem to be quite as simple as this). What am I missing?

Also, I'm not using fftconvolve at the moment because I want to compute the fft of A just once, then use it repeatedly for each b.

Cheers,

Alex

[Sturla Molden]

Den 28.06.2011 20:17, skrev Alex Flint:
> I am trying to perform 2d convolutions between a large 2d array A and
> a bunch of small 2d arrays B1...Bn. My approach is roughly:
>
> a = fft(A,size)
> for b in bs:
> ans = ifft(fft(b,size)*A)
> slow = convolve2d(A, b, 'same')
>
> However, as implemented above, ans is offset an inconsistent amount
> from the answer produced by convolve2d, presumably because convolve2d
> is treating b as if the origin is in the center whereas fft treats b
> as if the origin is at the top left (but it doesn't seem to be quite
> as simple as this). What am I missing?

You are not doing 2D convolution with the FFT. You want fft2 (or rfft2).

Sturla

_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

[Alex Flint]

oops, I am actually using fft2/ifft2, I just forgot to write it in my pseudo code

[Sebastian Berg]

Hi,

I assume its *a inside the loop, but don't you have to conjugate one of
the two?
a = np.conj(np.fft2(A))
I did this once, and the scipy fftconvolve seemed to forget the
conjugation?
Also use np.fft.fftshift to undo the shifting introduced by the fft.

Regards,

Sebastian

[Sturla Molden]

Den 30.06.2011 14:26, skrev Sebastian Berg:
> Hi,
>
> I assume its *a inside the loop, but don't you have to conjugate one of
> the two?

No, he just has to multiply in rectangular form.

He might be confused by circular connvolution though.

>>> import numpy as np
>>> a = np.zeros(100)
>>> b = np.zeros(100)
>>> a[1] = 1
>>> b[:10] = 2

>>> np.convolve(a,b)
array([ 0., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 0., 0.,

>>> np.round(np.fft.ifft(np.fft.fft(a)*np.fft.fft(b)).real)
array([ 0., 2., 2., 2., 2., 2., 2., 2., 2., 2., 2., 0., 0.,

[Edward Montague]

 Took awhile to locate this group.

 My question is about the 3d instance of fftconvolve.

I may have a correct  result, however there's no way for me to check this.

The resources  to do so just don't appear to be readily available from the internet.

[Edward Montague]

 To resolve some of these issues I'm examining the 1d instance, using the following code.

  One difference  between the numerical integration and the fftconvolve, is the magnitude of the

convolution ;  this is usually larger for fftconvolve. I have yet to find a way to bring them into alignment.

 ```

# Some code from :
 
# https://dspillustrations.com/pages/posts/misc/convolution-examples-and-the-convolution-integral.html
 
 
import matplotlib.pyplot as plt

import plotly.graph_objects as go
from plotly.subplots import make_subplots

import numpy as np
from scipy.fftpack import fftn, ifftn, fftshift
from scipy.integrate import simps
import scipy.signal as sg

# ----------------------------------------------------------------------

 
def showConvolution(f1, f2, t0):
    # Calculate the overall convolution result using Simpson integration
    convolution = np.zeros(len(t))
    for n, t_ in enumerate(t):
        prod = lambda tau: f1(tau) * f2(t_-tau)
        convolution[n] = simps(prod(t), t)
   
    # Create the shifted and flipped function
    f_shift = lambda t: f2(t0-t)
    prod = lambda tau: f1(tau) * f2(t0-tau)

    # Plot the curves
   
    #plt.subplot(211)
    #plt.plot(t, f1(t), label=r'$f_1(\tau)$')
   # plt.plot(t, f_shift(t), label=r'$f_2(t_0-\tau)$')
   # plt.plot(t, prod(t), 'r-', label=r'$f_1(\tau)f_2(t_0-\tau)$')
   
    # plot the convolution curve
    # plt.subplot(212)
    plt.plot(t, convolution, label='$(f_1*f_2)(t)$')
   
    # recalculate the value of the convolution integral at the current time-shift t0
    current_value = simps(prod(t), t)
   
   
    plt.plot(t0, current_value, 'ro')  # plot the point

# ----------------------------------------------------------------------
    
Fs = 50  # our sampling frequency for the plotting
Fs = 128
T = 5    # the time range we are interested in
t = np.arange(-T, T, 1/Fs)  # the time samples
 
   
#f1 = lambda t: 0.7*(abs(t)<0.5).astype(float)
#f2 = lambda t: 0.8*(abs(t-0.62)<0.75).astype(float)

f1 = lambda t: np.maximum(0, 1-abs(t))
f2 = lambda t: (t>0) * np.exp(-t)


plt.figure(figsize=(8,3))
plt.subplot(211)
plt.plot(t, f1(t), label='$f_1(t)$')
plt.plot(t, f2(t), label='$f_2(t)$')


for i in range(20):
    t0=-5+10*i/20
    showConvolution(f1, f2, t0)
   
# plt.show()

#  Now attempting discretization .

# t = np.arange(-T, T, 1/Fs)  # the time samples

# f1 = lambda t: np.maximum(0, 1-abs(t))
# f2 = lambda t: (t>0) * np.exp(-2*t)



def g1(t0):
    return np.maximum(0, 1-abs(t0))
   
   
def g2(t0):
    if(t0>0):
        return np.exp(-t0)
    else:
        return 0


# f1 = lambda t: (abs(t)<0.5).astype(float)
# f2 = lambda t: 0.8*(abs(t-0.2)<0.75).astype(float)

pwr=7
n=int(2**pwr)
#n=128 # 2^7
a1=np.zeros((n))
b1=np.zeros((n))

for i in range(n):
    t0=-5+10*float(i)/(n-1)
    a1[i]=g1(t0)
    b1[i]=g2(t0)
   
c1 = sg.fftconvolve(a1,b1, mode = 'same')  # use mode same
mc = c1.max()
print(' c1  max ')
print(mc)


#c1 /= pwr
#c1 /= 2
c1 /= mc
ma  = a1.max()
mb  = b1.max()

print(' a1  max ')
print(ma)

print(' b1  max ')
print(mb)


mb  = ma*mb
#c1 =  c1*mb/2#/((pwr+1)/2
#c1 =  c1*mb
#c1 /= mc
c1=c1/2
plt.subplot(212)
plt.plot(a1,'k',color='steelblue')
# plt.plot(t, f1(t), label='$f_1(t)$')


plt.plot(b1,'k',color='orange')
plt.plot(c1,'k',color='red')

plt.show()

c1=np.zeros((n))

c1=conv1(a1,b1,c1)

plt.plot(c1,'k',color='red')

plt.show()

```

[Edward Montague]

 Probably best, in the first instance, to divide fftconvolve by n; where n is the number of samples.

That's if it hasn't already been divided by this, within  the routine.

[Edward Montague]

Also dismiss conv1 as this routine wasn't included with the code.

[Edward Montague]

  So there's some correspondence between the continuous  and discrete convolutions.

 Now  as fftconvolve is a rather involved routine, I decided to use  something that I'm more familiar

with,, this just happens to showcase some of the capabilities  of plotly's 3d graphing.

```

#   Three volume plots ; using subplots.
#   Some from a 3d array .
#   Now using 3d array and fftconvolve .

#   A n x n  convolution allways returns a 2n result.

import plotly.graph_objects as go
from plotly.subplots import make_subplots

import numpy as np
from scipy.fftpack import fftn, ifftn, fftshift

import scipy.signal as sg

# ------------------------------------------------------------------

def convn(a1,b1):
    # Restrict data to the first half, then similar to mode = 'same' ?.
    # a1 and b1 real, easily updated to complex.
   
    A1=fftn(a1)
    B1=fftn(b1)
    C1=np.multiply(A1,B1)
    c1=ifftn(C1)
    c1=c1.real
    return c1
   
# ------------------------------------------------------------------
   

n=32
m=32
p=32


# Generate data, within first half only.
p1=p
m1=m
n1=n
p1=int(p1/4)
m1=int(m1/3)
n1=int(n1/2)

X, Y, Z = np.mgrid[:n, :m, :p]
vol = np.zeros((n, m, p))

# need to center the generated objects, somewhat ;
# add to the index a value.
ns=0
for k in range(p1):
    for j in range(m1):
        for i in range(n1):
            vol[i+ns,j+ns,k+ns]=1

volf = np.zeros((n, m, p))

for k in range(m1):
    for j in range(p1):
        for i in range(n1):
            volf[i+ns,j+ns,k+ns]=1


# ------------------------------

vol /= vol.max()
volf /= volf.max()

# Compare methods here .

#volg = sg.fftconvolve(vol,volf, mode = 'same')

volg = convn(vol,volf)


volg = np.real(volg)
volg /= volg.max()


ma = vol.max()
mb = volf.max()

mb = ma*mb
volg = volg*mb

# ------------------------- Graphing results ---------------------------

# Initialize figure with 3 3D subplots

fig = make_subplots(
    rows=3, cols=1,
    specs=[[{'type': 'Volume'}], [{'type': 'Volume'}],
           [{'type': 'Volume'}]],
          
    subplot_titles=("data1 3d", "data2 3d", "3d Real(fftconvolve)")      
           )
# ------------------------------------------------------------------

# adding volumes to subplots.
   
fig.add_trace(
    go.Volume(
    x=X.flatten(),
    y=Y.flatten(),
    z=Z.flatten(),
    value=vol.flatten(),
    isomin=0.1,
    isomax=1.0,
    opacity=0.1, # needs to be small to see through all surfaces
    surface_count=17 # needs to be a large number for good volume rendering
    ),
    row=1, col=1
    )   
   
   
fig.add_trace(
    go.Volume(
    x=X.flatten(),
    y=Y.flatten(),
    z=Z.flatten(),
    value=volf.flatten(),
    isomin=0.1,
    isomax=1.0,
    opacity=0.1, # needs to be small to see through all surfaces
    surface_count=17 # needs to be a large number for good volume rendering
    ),
    row=2, col=1
    )

fig.add_trace(
    go.Volume(
    x=X.flatten(),
    y=Y.flatten(),
    z=Z.flatten(),
    value=volg.flatten(),
    isomin=0.1,
    isomax=1.0,
    opacity=0.1, # needs to be small to see through all surfaces
    surface_count=17 # needs to be a large number for good volume rendering
    ),
    row=3, col=1
    )
   

fig.update_layout(
    title_text='3D subplots with different colorscales',
    height=1800,
    width=1024
)

fig.show()

# ------------------------------------------------------------------

'''
modestr {‘full’, ‘valid’, ‘same’}, optional

    A string indicating the size of the output:

    full

        The output is the full discrete linear convolution of the inputs.
        (Default)
       
    valid

        The output consists only of those elements that do not rely on the
         zero-padding. In ‘valid’ mode, either in1 or in2 must be at least
          as large as the other in every dimension.
         
    same

        The output is the same size as in1, centered with respect to the
         ‘full’ output.

'''

```

[Edward Montague]

 I'm now using python vedo 4.1  to plot volumes, may also try mayavi.

I wrote my own function to perform convolution, for me this is less confusing than fftconvolve.

```

#   Three volume plots from 3d arrays.
#
#   Now using 3d array, convn and fftconvolve .

#   A n x n  convolution allways returns a 2n result.

import numpy as np
from scipy.fftpack import fftn, ifftn, fftshift
import scipy.signal as sg

# ------------------------------------------------------------------

def convn(a1,b1):
    # Restrict data to the first half, then similar to mode = 'same' ?.
    # a1 and b1 real, easily updated to complex.
   
    A1=fftn(a1)
    B1=fftn(b1)
    C1=np.multiply(A1,B1)
    c1=ifftn(C1)
    c1=c1.real
    return c1
   
# ------------------------------------------------------------------
   

n=128
m=n
p=n

# Generate data, within first half only.
p1=p
m1=m
n1=n
p1=int(p1/4)
m1=int(m1/3)
n1=int(n1/2)

# X, Y, Z = np.mgrid[:n, :m, :p]

vol = np.zeros((n, m, p))

for k in range(p1):
    for j in range(m1):
        for i in range(n1):

            vol[i,j,k]=1

volf = np.zeros((n, m, p))

for k in range(m1):
    for j in range(p1):
        for i in range(n1):

            volf[i,j,k]=1


# ------------------------------

ma = vol.max()
mb = volf.max()

vol /= vol.max()
volf /= volf.max()

# Compare methods here .

#volg = sg.fftconvolve(vol,volf, mode = 'same')

volg = convn(vol,volf)

#volg = np.real(volg)
volg /= volg.max()

mb = ma*mb
volg = volg*mb
#
# ------------------------- Graphing results ---------------------------
#

from vedo import Volume, show

va = Volume(vol).mode(1).c('Dark2').alpha([0,0.6,0.8,1])
vb = Volume(volf).mode(1).c('Dark2').alpha([0,0.6,0.8,1])
vc = Volume(volg).mode(1).c('Dark2').alpha([0,0.6,0.8,1])

show([va,vb,vc], shape=[1,3], axes=1)

quit()

 

```