How to fit built in filter functions

I have a signal (e.g. square wave) that I am passing through a "black box" resulting in a modified output signal (e.g. saw tooth). I want to deduce the frequency response of the black box using Matlab's filter functions (e.g. butter, chebyl, ellip). How can I fit the filter inputs to the output when the input is variable?

In otherwords, how can I build a function that can fit the filter parameters in order to transform a known input signal into a known output signal in order to model the "black box".

I want to fit, [xout, yout] = func(beta, xin, yin), where beta are the filter parameters and xin == xout. I only need xin so that I know the sampling frequency of yin and therefore know the frequency components of the input signal.

Proposed steps

1. Read in input signal
2. Read in output signal
3. Apply filter to input signal
4. Compare filtered signal to output signal
5. Adjust filter parameters
6. Apply filter to input signal
7. Compare filtered signal to output signal
8. repeat until filtered and output signals match