Most computationally efficient way to find convolution of a matrix kernel with impulse response? Let say if we wish to filter an input sequence x[n1, n2, n3] of NxNxN points using an Linear Shift Invariance system with impulse response h[n1, n2, n3], where the filter is a separable sequence, h[n1, n2, n3] = a[n1]b[n2]c[n3].  
How can we develop a computationally efficient way to obtain the    output y[n1, n2, n3]?
My thought:  since we can represent the Fourier Transform of the Filter as: 
$$ H(w_1, w_2, w_3) = \sum \sum \sum h[n_1,n_2,n_3]*e^{-jw_1n_1} *e^{-jw_2n_2}*e^{-jw_3n_3} $$
$$ H(w_1, w_2, w_3) = A(w_1)B(w_2)C(w_3)$$, where
$$ A(w_1) =   \sum h[n_1,n_2,n_3]*e^{-jw_1n_1} $$
$$ B(w_2) =   \sum h[n_1,n_2,n_3]*e^{-jw_2n_2} $$
$$ C(w_3) =   \sum h[n_1,n_2,n_3]*e^{-jw_3n_3}$$
then,
$$ y(n_1, n_2, n_3) = x(n_1, n_2, n_3)*h(n_1, n_2, n_3)$$
$$ Y(w_1, w_2, w_3) = X(w_1, w_2, w_3)H(w_1, w_2, w_3)$$
$$ Y(w_1, w_2, w_3) = X(w_1, w_2, w_3)A(w_1)B(w_2)C(w_3)$$
I got stuck in this point and I am not too sure how to proceed through the process.  I don't quite see how to separate the sequence and how to multiply through the vectors to get the convolution sum of $Y(w_1, w_2, w_3)$
can anyone please guide me to the right direction?  Thanks!

x * h = IFT[FT[x * h] and FT[x * h]  is Ft[x] FT[h]. Since h is separable, it might be computationally efficient to expand its transform as a convolution Ft[x]*FT[b]*FT[c].
