Characterizing an unknown LTI system
One option from the realms of adaptive filtering is the Least Mean Squares (LMS) filter depicted below:



The idea is you take the output of the unknown system, compare it with the output of your adaptive filter and minimize the difference by tweaking the filter coefficients, using a LMS algorithm. When the error $e(n)$ is zero (or more often, lower than a threshold you define) you found your unknown system's response.

There are countless resources where you can learn about LMS. In this case, given

the input vector: $$X(i) = [x(i) \space x(i-1) \ldots x(i-N+1)]^T$$ and the coefficients vector: $$W = [w_0 \space w_1 \ldots w_{N-1}]^T$$

The output of the adaptive filter is: $$\hat{y}(i)=X(i)^TW$$

Now define the cost function to be: $$e(i)^2 = \big(d(n)-\hat{y}(i)\big)^2$$
It can be shown that the gradient of the cost function with respect to the coefficients is:$$\hat{\nabla}(i) = \frac{\partial e(i)^2}{\partial W} = -2e(i)X(i)$$
Since the gradient points towards the maximum of the cost function you update the coefficients in the opposite direction: $$\hat{W}(i+1) = \hat{W}(i)-\mu\hat{\nabla}(i) = \hat{W}(i)+2\mu e(i)X(i)$$
where $\mu$ controls the adaptation speed. I've ignored a lot of details but this is the idea.

Note that the filter will produce a frequency response that is adapted to you input signal $x(n)$. It's not guaranteed for any other frequency. For example, in the response depicted below, my $x(n)$ is composed of 3 tones. You can see only at those particular frequencies does the filter (N=3,6,10) match the true response (H1(z)). So consider feeding it a wideband signal like white noise.



If you want to learn more you can check for example:

Wikipedia article

B. Widrow, Adaptive Signal Processing, Prentice-Hall, 1985