In this paper, using Miller's approach and Dougall's identity, we derive new infinite series representations for the quadrivariate Nakagami-m joint density function, cumulative distribution function (cdf) and characteristic functions (chf). The classical joint density function of exponentially correlated Nakagami-m variables can be identified as a special case of the joint density function obtained here. Our results are based on the most general arbitrary correlation matrix possible. Moreover, the trivariate density function, cdf and chf for an arbitrary correlation matrix are also derived from our main result. Bounds on the series truncation error are also presented. Finally, we develop several representative applications: the outage probability of triple branch selection combining (SC), the moments of the equal gain combining (EGC) output signal to noise ratio (SNR) and the moment generation function of the generalized SC(2,3) output SNR in an arbitrarily correlated Nakagami-m environment. Simulation results are also presented to verify the accuracy of our theoretical results.