We consider the problem of decomposing an even order symmetric tensor with positive eigenvalues, into a sum of rank-1 components (also knows as canonical polyadic or CP decomposition). These tensors naturally arise in many signal processing applications (such as ICA, blind source separation, source localization) when we compute higher order cumulants of the measurements. We show that when these components possess certain harmonic structures, it is possible to design clever sampling techniques to identify $O(N^{2q})$ rank-1 factors using a tensor of order $2q$ (where $q$ is an integer) and dimension $N$. This is made possible by exploiting the idea of a higher order difference set that can be associated with the cumulant tensors. For unstructured even order tensors, we show that under some mild conditions, the problem of CP decomposition is equivalent to solving a system of quadratic equations as long as the rank of the tensor is $O(N^q)$. We finally propose two algorithms, one based on convex relaxation, and the other utilizes non-convex, Jacobi iteration to solve the resulting quadratic system.