Filters change the relative amplitudes of frequency components in a
signal. The problem of finding coefficients to achieve a particular
frequency response is solved; in the literature, by available programs,
and by a little known method invented by (I believe) Greg Berchin, who
frequents this group. To use them, you need to determine the frequency
response or alternatively the impulse response that will solve your
problem. Presumably, that will mimic the impulse response of the system
you want to model.

The correlation-based method that was mentioned probably
uses the following property.  Suppose that the input
signal is x[n], the system h[n], and the output y[n].

Let rxx[n] be the autocorrelation of x[n], and
rxy[n] be the cross-correlation of x[n] and y[n].
Then it can be shown that rxy[n] = rxx[n] * h[n].

So if you can estimate rxx and ryx, then you can
do a de-convolution to get h[n].  Further, it can
be shown that ryy[n] = rxx[n] * (h[n] * h[-n]).

Note that in this last one, the phase information
is gone and you can only solve up to h[n]*h[-n].
Since you are interested in FIR estimates, most
of the time people do what is called minimum-phase
decomposition.

Anyways, these are just a bunch of ideas,
extrapolated from your question.

- AR models (What does this stand for? Is it autoregressive?)
- The autocorrelation method I'd heard about (thanks for the
explanation)
- LMS Filters (I've found some good stuff at wikipedia to get me
started)

Thanks all!

A couple good books on these topics:

"Statistical Digital Signal Processing and Modeling" by Hayes
"Adaptive Filter Theory" by Haykin

You are given x and y, assuming that y = h * x + e  (* = convolution),
e is some zero-mean noise process and h the impulse response of an LTI
system. Determine h.

One approach not mentioned yet is via inverse filtering. If you can
choose x (this is not always the case, or sometimes the possible
inputs x are limited), there are much better choices than the obvious
impulse, because you have to somehow get rid of the noise e.

Let's say x was chosen such that there existed an inverse filter x^-1,
with

x * x^-1 = delta_n,

where delta_n is the impulse delayed by n samples. Then h can be found
by convolving y with x^-1. The longer x^-1, the more the noise gets
averaged, and the better your approximation of h.

Also, before you mix everything up: there are two ways to approach
this task. One is to stipulate a model for h. A model is essentially
concerned with data compression. You want to explain the
characteristics of h with as few numbers (called parameters) as
possible. AR modeling is just one of many alternatives.

System ID via Wiener filter (of which LMS is just one of many possible
adaptive variants), cross-correlation methods (in time or frequency
domain) or inverse filtering gives a direct approximation of h. No
data compression. This requires more effort (for the same accuracy),
because you have to estimate many more numbers than for modeling.

Obviously, both methods have their pros and cons. The best choice
mainly depends on what you intend to do use h for.