Recover Filter Coefficients from Filtered Noise
It's an interesting problem.  What you have there is what's known as a blind deconvolution problem.  These are well known "hard" problems, but not necessarily impossible.  Finding an algorithm to solve it relies on using some prior knowledge you have about the filter or the noise source driving it.  

It's an ill-posed problem mathematically, so if there is a solution, there's no guarantee it's unique.  In other words you might get a set of FIR coefficients that seem to fit the data, but it's not guaranteed that they are the "right" coefficients that model the filter well in the real world.

Having said that, I've tried using the EM algorithm and got some promising results.  It's not perfect, but I think it's possible to get a solution.  

I've changed the model slightly so that $$S = N*F + v$$where $v$ is Gaussian noise with some unknown variance.  If we find a good solution, this variance should end up being small compared to the signal.

In the EM algorithm you need some hidden variables: in this case, $N$, and some unknown parameters.  Those are $F$, $\sigma^2$ the unknown noise variance, and $p$ the probability of $N_i\in{0,1}$ being equal to one.  This is the prior distribution on $N$.

The trick to it I found was to use the alternative description of the EM algorithm where the hidden variables are described by an approximating distribution $q$.  We want $q$ to approximate the posterior distribution given the observations and the parameters: $p(N|S,F,p,\sigma^2)$.  This almost certainly doesn't factorize as a product of probabilities $q(N_0)q(N_1)...q(N_n)$ so this is an approximation.  

In this case $q$ consists of $q_i$ for $i=0..n$ where $q_i$ is the probability that $N_i=1$ under this approximating distribution.  

The EM algorithm is an iterative one, where each iteration involves alternately updating $q$ based on the current value of the parameters, then updating the parameters based on the current value of $q$.  

I haven't said anything about $m$ yet: it's unknown, but if we make it too big, the extra FIR coefficients should be set to near zero.  These might be at the start or the end, it just depends where the algorithm decides to converge to.

I simulated some data of my own with the same pulse shape as in your synthetic example, with these results.  The true value of $m$ was 200, so I ran the EM algorithm with $m=400$ to see how it did.  I made $S$ 10000 samples long, to be sure I had enough data.  The result came out at about half the amplitude of the original, and as explained above there are some extra coefficients at the start and the end, but shifting and scaling to compensate you can see that it got the shape of the impulse response fairly accurately:


Zooming in on the middle section:


From your file synthetic_data.txt I got something similar:





As it stands it seems to be overestimating $p$ i.e. putting in too many impulses and underestimating the magnitude of the FIR coefficients to compensate.  

Now on test_data.txt, the algorithm gave this impulse response:



It's harder to tell how well the algorithm has done this time because we don't know the true impulse response.  One thing we can do it try to reconstruct $S$ from our estimate $\hat F$ and some estimate of $N$.  

I've used $q_i$ to estimate $N_i$ by thresholding: if $q_i&gt;0.5$ I put $\hat N_i=1$ as my estimate of $N_i$, otherwise I use $N_i = 0$.

Convolving that sequence $\hat N$ with $\hat F$ gave a decent match to $S$.  I've zoomed in on a typical short subsequence so you can see the detail:


One caveat I have is that the rate of impulses seem a bit high: about 10% of the reconstructed $\hat N_i$ samples were equal to $1$, which is more than your 1 in 20 to 40 estimate.  

But overall I think that's doing well for a blind deconvolution problem, so I'd say it's doable with an EM algorithm or something similar.

There are quite a few other algorithms for blind deconvolution, see the references, and one of the others might be better for you.  Variational Bayesian EM would be worth trying, and there are others which might have some advantages in other ways.  It always depends on how much prior knowledge you have, and on other constraints like runtime.

The code below works in Octave: I haven't tried it in Matlab (because I don't have a licence) but it should work in Matlab too.  Very much prototype code, don't rely on it for anything critical without some more work (sanity checking, catching potential divide-by-zeros etc.)  Runs in about 30 minutes on my two-year-old laptop.

Thanks for uploading the example datasets, which was very useful.

function [Fhat,q] = blindDeconvImpulseNoise(filename)

% Load data from file
S = load(filename);
% Renormalize
scaleFac = sqrt(mean(S.^2));
S = S/scaleFac;

% Estimate (assume filter has less than 400 taps)
tic();
[Fhat,q] = estimate(S,400);
toc();

Fhat = Fhat * scaleFac;

end


% Run the estimation: the EM algorithm
% (see ftp://ftp.cs.toronto.edu/pub/radford/emk.pdf)
function [Fhat,q] = estimate(S,m)

L = length(S);

% Initial values
q                      = ones(L+m-1,1)/2;

% Fhat is the estimated filter.  
Fhat = zeros(m,1);
Fhat(max(1,floor(m/2))) = 1;

% Initial values of the scalar parameters: noise variance and 
sigmasq_hat = 1;
Rhat = 0.01;

last_sig = inf;

starttime = now;
timeout   = 1800;

while 1

  % Estimate hidden variables
  % Start by converting hidden states into 
  Fmtx  = getConvmtx(Fhat,length(q));
  Fmtx  = Fmtx(m:end-m+1,:);

  for i=1:length(q)

    ind = max(1,i - m):min(i+m-1,size(Fmtx,1));
    ind2 = max(1,i - 2*m):min(i+2*m-1,size(Fmtx,2)); 
    ind2 = ind2(ind2~=i);
    Fi = Fmtx(ind,i);
    resid_less_i = S(ind) - Fmtx(ind,ind2) * q(ind2);
    l = (Fi'*Fi - 2*Fi'*resid_less_i)/(2*sigmasq_hat) ...
        - log(Rhat) + log(1-Rhat);

    q(i) = 1./(1+exp(l));

  end

  qmtx                  = getConvmtx(q,m);
  qmtx                  = qmtx(m:end-m+1,:);

  % Estimate parameters
  Sigma                 = (qmtx'*qmtx) ;
  % Matrix diagonal
  Sigma(1:m+1:end)      = sum(qmtx);
  % Store common term Q'*S
  NS                    = (qmtx'*S);
  % Pseudoinverse estimate of Fhat
  Fhat                  = Sigma\NS;
  % Noise estimate, based on error residual
  sigmasq_hat           = (S'*S - NS'*Fhat)/length(S);
  % Pulse rate estimate
  Rhat                  = mean(q);

  if abs(sigmasq_hat - last_sig) &lt; 1e-10
     % No change on this iteration: stop here
     break
  end
  last_sig = sigmasq_hat;

  if 24*60*60*(now - starttime) &gt; timeout
     warning('Timeout!');
     break
  end
end

end


function Nmtx = getConvmtx(N,m)

Nmtx = zeros(length(N)+m-1,m);

for idx=1:m
    Nmtx(idx:end-m+idx,idx) = N;
end

end