Well, I did it a few times...
If individual pile stiffness matrix coefficients are:
k11 - axial (along X-axis) [k/ft]
k22 - translation along Y-axis [kip/ft]
k33 - translation along Z-Axis [kip/ft]
k44 - torsional [kip-ft/rad]
k55 - bending about Y-axis [klp-ft/rad]
k66 - bending about Z-axis [klp-ft/rad]
k26, k62, k35, k53 - off diagonal terms [kip] (or neglect)
Center of pilecap coordinates Y=0, Z=0
and y, z - individual pile coordinates,
Then equivalent matrix stiffness of Pile foundation:
K11 = Sum(k11)
K22 = Sum(k22)
K33 = Sum(k33)
K44 = Sum(k44)+Sum(k11*(y^2+z^2))
K55 = Sum(k55)+Sum(k11*z^2)
K66 = Sum(k66)+Sum(k11*y^2)
K26=K62=Sum(k26) (NEGATIVE!)
K35=K53=Sum(k35) (POSITIVE!)
The trouble comes when you will find out that each of your piles has different stiffness due to pile group effect and stiffer outside piles are overloaded. But it's a whole different issue...
Yakov