I am trying to program a DSP filter to study prototypes of loudspeaker filters, which will be implemented as passive analog circuits in the final speaker system. This process involves simulation of the electro-acoustics of the entire loudspeaker system, taking into account the real-world impedance curves of the loudspeaker drivers, their acoustic frequency response, and the acoustics of the loudspeaker enclosure. The resulting transfer curves of the (simulated) analog filters therefore do not follow any "textbook filter functions".
I tried fitting IIR/biquad filter coefficients to the target function defined by the analog filter. To this end I used the invfreqz function in GNU Octave. This method can achieve a very good fit of the IIR/biquad filter function to the target function of the analog filter. However, the IIR/biquad filters are usually not stable, so they are useless in my situation.
I think the Matlab implementation of the invfreqz function is a bit better at estimating stable filters, but I guess that might not be the definitive answer to the problem (I don't have access to Matlab, so I can't check). I would also like to implement the solution using free software, possibly by taking a somewhat different approach than "brute-fore application of invfreqz".
I am not very experienced with DSP filter methods, so I am a bit clueless as how to tackle this problem. What could be a practical approach to estimate IIR/biquad filter coefficients that give a transfer function that is "sufficiently similar" to an arbitrary target function defined by an analog filter?