How to use the spectral.T1FFT.T1FFT function in spectral

def synthesize(self, f_hat):

        f2f = T1FFT.synthesize(f_hat.conj(), axis=1).conj()

        # Multiply f_2 by a phase factor:
        m_min = -np.floor(f2f.shape[2] / 2)
        m_max = np.ceil(f2f.shape[2] / 2) - 1
        psi = np.linspace(0, 2 * np.pi, f2f.shape[1], endpoint=False)  # may not need this many points on every circle
        factor = np.exp(1j * psi[:, None] * np.arange(m_min, m_max + 1)[None, :])

        f2 = f2f * factor[None, ...]

        f1p = T1FFT.synthesize(f2.conj(), axis=2).conj()

        f1c = self.resample_p2c_3d(f1p)


        # delta = -0.5  # we're shifting from [0, 1) to [-0.5, 0.5)
        # xi1 = np.arange(-np.floor(f1c.shape[0] / 2), np.ceil(f1c.shape[0] / 2))
        # xi2 = np.arange(-np.floor(f1c.shape[1] / 2), np.ceil(f1c.shape[1] / 2))
        # XI1, XI2 = np.meshgrid(xi1, xi2, indexing='ij')
        # phase = np.exp(-2 * np.pi * 1j * delta * (XI1 + XI2))
        # f1c_shift = f1c / phase[:, :, None]


        #f = T2FFT.synthesize(f1c, axes=(0, 1))
        #f = T2FFT.synthesize(f1c_shift, axes=(0, 1))
        f = shift_ifft(f1c)

#delta = -0.5  # we're shifting from [0, 1) to [-0.5, 0.5)
        #xi1 = np.arange(-np.floor(f1c_shift.shape[0] / 2.), np.ceil(f1c_shift.shape[0] / 2.))
        #xi2 = np.arange(-np.floor(f1c_shift.shape[1] / 2.), np.ceil(f1c_shift.shape[1] / 2.))
        #XI1, XI2 = np.meshgrid(xi1, xi2, indexing='ij')
        #phase = np.exp(-2 * np.pi * 1j * delta * (XI1 + XI2))
        #f1c = f1c_shift * phase[:, :, None]

        f1c = shift_fft(f)

        # Change from Cartesian (c) to a polar (p) grid:
        f1p = self.resample_c2p_3d(f1c)
        # This gives f1p[r, varphi, theta]

        # FFT along rotation angle theta
        # We conjugate the argument and the ouput so that the complex exponential has positive instead of negative sign
        f2 = T1FFT.analyze(f1p.conj(), axis=2).conj()
        # This gives f2[r, varphi, q]
        # where q ranges from q = -floor(f1p.shape[2] / 2) to q = ceil(f1p.shape[2] / 2) - 1  (inclusive)

        # Multiply f2 by a (varphi, q)-dependent phase factor:
        m_min = -np.floor(f2.shape[2] / 2.)
        m_max = np.ceil(f1p.shape[2] / 2.) - 1
        varphi = np.linspace(0, 2 * np.pi, f2.shape[1], endpoint=False)  # may not need this many points on every circle
        factor = np.exp(-1j * varphi[None, :, None] * np.arange(m_min, m_max + 1)[None, None, :])
        f2f = f2 * factor

        # FFT along polar coordinate of frequency domain
        f_hat = T1FFT.analyze(f2f.conj(), axis=1).conj()
        # This gives f_hat[r, p, q]

        return f, f1c, f1p, f2, f2f, f_hat