You can equalize magnitude and phase simultaneously by defining a desired complex frequency response
with magnitude $M(\omega)$ and phase $\phi(\omega)$ chosen such that they compensate for the given magnitude and phase distortions.
An FIR filter approximating $(1)$ can be designed by using the following error function:
where $\omega_k$ are frequencies on a dense grid, $W_k$ are some non-negative weights allowing to emphasize certain frequency regions compared to others, and $H(\omega)$ is the frequency response of the FIR filter to be designed:
Minimizing $(2)$ with respect to the filter coefficients $h[n]$ can be achieved by solving a system of linear equations. I've written a Matlab/Octave function that does just that: lslevin.m
Note that it may be necessary to add a linear phase (i.e., a delay) to the required $\phi(\omega)$ in order for a causal filter to better approximate the desired complex frequency response. As a guideline, the average delay implied by the desired phase should be about half the chosen filter length $N$. It is worthwhile experimenting with different additional delays because they may have a significant influence on the resulting approximation error.