Как мне умножить выходные векторы scipy.fftpack вместе?

Вам не нужно возвращаться к np.float64 и hstack. Вы можете создать пустой целевой массив той же формы, что и sfft.rfft(Y) и sfft.rfft(X), затем создать его представление np.complex128 и заполнить это представление результатом умножения. Это автоматически заполнит целевой массив по желанию.
Если я повторю ваш пример:

import numpy as np
import scipy.fftpack as sfft

X = np.random.normal(size = 2000)
Y = np.random.normal(size = 2000)
Xf = np.fft.rfft(X)
Xf_cpx = Xf[1:-1].view(np.complex128)
Yf = np.fft.rfft(Y)
Yf_cpx = Yf[1:-1].view(np.complex128)

Zf = np.empty(X.shape)
Zf_cpx = Zf[1:-1].view(np.complex128)

Zf[0] = Xf[0]*Yf[0]

# the [...] is important to use the view as a reference to Zf and not overwrite it
Zf_cpx[...] = Xf_cpx * Yf_cpx 

Zf[-1] = Xf[-1]*Yf[-1]

Z = sfft.irfft.irfft(Zf)

вот и все! Вы можете использовать простой оператор if, если хотите, чтобы ваш код был более общим и обрабатывал нечетные длины, как описано в ответе Хайме. Вот функция, которая делает то, что вы хотите:

def rfft_mult(a,b):
    """Multiplies two outputs of scipy.fftpack.rfft"""
    assert a.shape == b.shape
    c = np.empty( a.shape )
    c[...,0] = a[...,0]*b[...,0]
    # To comply with the rfft support of multi dimensional arrays
    ar = a.reshape(-1,a.shape[-1])
    br = b.reshape(-1,b.shape[-1])
    cr = c.reshape(-1,c.shape[-1])
    # Note that we cannot use ellipses to achieve that because of 
    # the way `view` work. If there are many dimensions, one should 
    # consider to manually perform the complex multiplication with slices.
    if c.shape[-1] & 0x1: # if odd
        for i in range(len(ar)):
            ac = ar[i,1:].view(np.complex128)
            bc = br[i,1:].view(np.complex128)
            cc = cr[i,1:].view(np.complex128)
            cc[...] = ac*bc
    else:
        for i in range(len(ar)):
            ac = ar[i,1:-1].view(np.complex128)
            bc = br[i,1:-1].view(np.complex128)
            cc = cr[i,1:-1].view(np.complex128)
            cc[...] = ac*bc
        c[...,-1] = a[...,-1]*b[...,-1]
    return c

Вы можете просмотреть фрагмент возвращаемого массива, например:

>>> scipy.fftpack.fft(np.arange(8))
array([ 28.+0.j        ,  -4.+9.65685425j,  -4.+4.j        ,
        -4.+1.65685425j,  -4.+0.j        ,  -4.-1.65685425j,
        -4.-4.j        ,  -4.-9.65685425j])
>>> a = scipy.fftpack.rfft(np.arange(8))
>>> a
array([ 28.        ,  -4.        ,   9.65685425,  -4.        ,
         4.        ,  -4.        ,   1.65685425,  -4.        ])
>>> a.dtype
dtype('float64')
>>> a[1:-1].view(np.complex128) # First and last entries are real
array([-4.+9.65685425j, -4.+4.j        , -4.+1.65685425j])

Вам нужно будет обрабатывать четные или нечетные БПФ по-разному:

>>> scipy.fftpack.fft(np.arange(7))
array([ 21.0+0.j        ,  -3.5+7.26782489j,  -3.5+2.79115686j,
        -3.5+0.79885216j,  -3.5-0.79885216j,  -3.5-2.79115686j,
        -3.5-7.26782489j])
>>> a = scipy.fftpack.rfft(np.arange(7))
>>> a
array([ 21.        ,  -3.5       ,   7.26782489,  -3.5       ,
         2.79115686,  -3.5       ,   0.79885216])
>>> a.dtype
dtype('float64')
>>> a[1:].view(np.complex128)
array([-3.5+7.26782489j, -3.5+2.79115686j, -3.5+0.79885216j])