Independent component analysis (ICA) is a popular technique for demixing multi-channel data. The performance of a typical ICA algorithm strongly depends on the presence of additive noise, the actual distribution of source signals, and the estimated number of non-Gaussian components. Often a linear mixing model is assumed and source signals are extracted by data whitening followed by a sequence of plane (Jacobi) rotations. In this article, we develop a novel algorithm, based on the quaternionic factorization of rotation matrices and the Newton-Raphson iterative scheme. Unlike conventional rotational techniques such as the JADE algorithm, our method exploits $4 \times 4$ rotation matrices and uses approximate negentropy as a contrast function. Consequently, the proposed method can be adjusted to a given data distribution (e.g. super-Gaussians) by selecting a suitable non-linear function that approximates the negentropy. Compared to the widely-used, the symmetric FastICA algorithm, the proposed method does not require an orthogonalization step and is more accurate in the presence of multiple Gaussian sources.