We propose a high order numerical method for computing time dependent 4-D Wigner equation with unbounded potential and study a canonical quantum double-slit interference problem. To address the difficulties of 4-D phase space computations and higher derivatives from the Moyal expansion of nonlocal pseudo-differential operator for unbounded potentials, an operator splitting technique is adopted to decompose the 4-D Wigner equation into two sub-equations, which can be computed analytically or numerically with high efficiency. The first sub-equation contains only linear convection term in $(\bm x, t)$-space and can be solved with an advective method, while the second involves the pseudo-differential term and can be approximated by a plane wave expansion in $\bm k$-space. By exploiting properties of Fourier transformation, the expansion coefficients for the second sub-equation have explicit forms and the resulting scheme is shown to be unconditionally stable for any higher derivatives of the Moyal expansion, ensuring the feasibility of the 4-D Wigner numerical simulations for quantum double-slit interferences. Numerical experiments demonstrate the spectral convergence in $(\bm x, \bm k)$-space and provide highly accurate information on the number, position, and intensity of the interference fringes for different types of slits, quantum particle masses, and initial states (pure and mixed).