get poles and zeros of frequency response
This is a very complicated problems and I don't think there exists a one-size-fits-all solution. You can try Matlab's $invfreqz()$ and see if it works for your purposes https://www.mathworks.com/help/signal/ref/invfreqz.html

In general this is a error minimization problem but the actual data and the way you set the your error function and the search strategy makes a huge difference in the outcome. Things to watch out for


How good is your data? What's the frequency coverage range, density and shape of frequency grid, signal to noise ratio, etc.
What do you care about in the result. Dips vs peaks ? Phase vs amplitude, ?
Stability and causality constraints. How do you deal with non-causalities in the data (real or noise induced)
Search space representation, search strategy &amp; algorithms. 
Error surface shaping and optimization.
Convergence tracking, stop criteria, metrics for success (or lack thereof)


I developed one of these that is optimized around the specific constraints and requirements of my application and found invfreqz() to be totally useless. Your mileage may vary.