Modelling an existing DSP filter by analysing the differences between its inputs and outputs
If I understood correctly, you have a signal $x[n]$ and an unknown discrete-time LTI filter, so that you can look at the output of the filter $y[n]$. Now you are looking for the impulse response of the filter $h[n]$.

The output follows a convolution rule
$$y[n]=x[n]*h[n]=\sum_{m=-\infty}^{\infty}x[m]h[n-m]$$
The process to find $h[n]$, given $x[n]$ and $y[n]$ is called deconvolution. There are some approaches to do that.

For example, if the Fourier transform of the signals exist, using convolution theorem, the convolution can be converted to product in the frequency domain: $$Y(\omega)=X(\omega)H(\omega)$$
Therefore, $$H(\omega)=\frac{Y(\omega)}{X(\omega)}$$
You should however be careful about the zeros of $X(\omega)$. If $x[n]$ can be an option, then you should choose it such that it can excite the system at all frequencies such that you can probe the desired frequency range that you are interested in $H(\omega)$. Alternatively, you can excite the filter in repeated number of times at different frequencies, which can be interpreted as sampling the frequency response. This will give you a sampled spectrum, so if you know the filter does not have singularities between the sample frequencies, then you can model it by an LTI filter with your desired order (depending on your acceptable approximation error).  

In response to your second question, after you acquired an approximation such as $\hat{h}[n]$ for the impulse response, then you can apply it to any desired input by using convolution.