Finding the Transfer Function of a Multiplicative Distortion in MATLAB
Basically the transfer function is given by (Continuous time case):
$$ H \left( j \omega \right) = K \frac{ \left( j \omega - {q}_{1} \right) \left( j \omega - {q}_{2} \right) }{ \left( j \omega - {p}_{1} \right) \left( j \omega - {p}_{2} \right) } $$
Now, given your frequency data sampled an $ N $ points $ {\left\{ {\omega}_{i} \right\}}_{i = 1}^{N} $ we need to use the discrete-time form of the transfer function over the Z domain (Z Transform).
Hence the model becomes:
$$ H \left( z \right) {\mid}_{z = {e}^{j \omega}} = K \frac{ \left( z - {q}_{1} \right) \left( z - {q}_{2} \right) }{ \left( z - {p}_{1} \right) \left( z - {p}_{2} \right) }{\Bigg|}_{z = {e}^{j \omega}} = K \frac{ \left( {e}^{j \omega} - {q}_{1} \right) \left( {e}^{j \omega} - {q}_{2} \right) }{ \left( {e}^{j \omega} - {p}_{1} \right) \left( {e}^{j \omega} - {p}_{2} \right) } $$
you may set a system of $ N $ complex equations and use non linear solver to find the 5 parameters ($ K $, $ {q}_{1} $, $ {q}_{2} $, $ {p}_{1} $, $ {p}_{2} $)  you need.
In case the measurements are noisy you may use Least Squares.
Remark
The above solves the problem from the data of @jmcd in frequency domain. In practice, using the data domain one could use simpler System Identification (See System Identification: An Overview) methods to get the equivalent ARMA model.
Regarding Frequency Domain based methods for System Identification have a look at Rik Pintelon - System Identification: A Frequency Domain Approach.