I have the following transfer function
$$ H(s) = \frac{R}{sRC + 1} $$
where $V$ is the output and $I$ is the input. We can write
$$ V(s) = I(s) H(s) $$
or
$$ V( j \omega ) = I(j \omega) H(j \omega) $$
Now I have the output $V(t)$ as array of samples. As I have written above, I also know the transfer function. In this case, how can I obtain the input $I(t)$ as array of samples?
how can I obtain the input I(t) as array of samples?
Not easily.
Let's look at the transfer function from voltage to current
$$H(\omega) = \frac{sRC+1}{R} = 1/R + sC$$
This is basically the sum of a constant and an ideal differentiator. If you have $v(t)$ as an closed form expression, you can simply solve it as
$$i(t) = v(t)/R + C\frac{\partial v(t)}{\partial t}$$
I have the output V(t) as array of samples.
Assuming you have sampled this at a uniform interval we can express this as
$$v[n] = v(nT)$$
where T is the sample period. That's where the sampling theorem comes in: The sample rate needs to higher then twice the highest frequency in the signal or you get aliasing.
Here is where it get's tricky: the transfer function is decidedly NOT bandlimited. Quite the contrary it approaches infinity for $\omega \rightarrow \infty$ so it can't be sampled without significant amount of aliasing.
Your best shot is to somehow approximate the differentiator in discrete time. There are multiple ways of doing this, but the easiest would be something like
$$i[n] = v[n]/R + C/T \cdot (v[n]-v[n-1])$$
This will give good results for low frequencies but will be very wrong at high frequencies. The trick here is to choose a sample rate that's high enough so that the "wrong frequencies" are outside your band of interest.
The matlab and python code below seems to do what you need.
The plot above has the output voltage in black, the input current in magenta, and the estimated current in green circles. As you can see the input and estimated currents overlay pretty well.
The main problem is that the transfer function:
$$ H(s) = \frac{R}{sRC + 1} $$
does not have a proper inverse.
So I used the suggestion here and replaced $s$ with
$$\frac{s}{1 + 0.0001s}$$
which approximates the derivative $s$ with something that's proper.
syms s;
R = 10;
C = 1* 10^-6;
H = R/(R*C*s + 1);
[Hnum,Hden] = numden(H);
Htf = tf(sym2poly(Hnum), sym2poly(Hden));
Hinv = inv(H);
% https://stackoverflow.com/a/20539050/12570
proper_derivative = s/(1 + 0.0001*s);
Hinv2 = subs(Hinv,{s},{proper_derivative});
[Hinv2_num,Hinv2_den] = numden(Hinv2);
Hinv2tf = tf(sym2poly(Hinv2_num),sym2poly(Hinv2_den));
t = [0:999]*0.01;
i_in = cumsum(randn(1,1000));
v_out = lsim(Htf, i_in, t);
i_est = lsim(Hinv2tf, v_out, t);
clf
plot(v_out, 'k');
hold on;
plot(i_est, 'go');
plot(i_in, 'm');
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
R = 10
C = 0.000001
H = signal.lti(R,[R*C,1])
t = np.linspace(0, 5, 1000)
i_in = np.cumsum(np.random.normal(0,1,1000))
tout, v_out, x = signal.lsim(H, i_in, t)
Hinv = signal.lti([R*C, 1], R)
s/(1000000*(s/10000 + 1)) + 1/10
Hinv2 = signal.lti([11, 100000], [100, 1000000])
tout2, i_est, x2 = signal.lsim(Hinv2, v_out, t)
plt.plot(t, v_out, 'k')
plt.plot(t, i_est, 'go')
plt.plot(t, i_in, 'm')
