It is easily possible to design FIR filters with a complex frequency response specification, i.e., with a desired frequency response
$$D(\omega)=|D(\omega)|e^{j\phi(\omega)}\tag{1}$$
where both the desired magnitude $|D(\omega)|$ as well as the desired phase $\phi(\omega)$ can be prescribed. A least squares approximation of $(1)$ by an FIR filter just involves the solution of a system of linear equations. I've written a Matlab/Octave function solving that problem: lslevin.m.
Note that it is important to choose the specification $(1)$ in such a way that it can be realistically achieved by an FIR filter of a given length $N$. E.g., the average delay implied by the desired phase $\phi(\omega)$ should be close to $(N-1)/2$, which is the delay of a linear phase FIR filter of length $N$.