Recover Filter Coefficients from Filtered Noise I have a digital signal which may be represented as a pulse noise source filtered with an FIR (finite impulse response) filter. Suppose that the noise consists of discrete pulses (nonzero samples where most samples are zero), and there is an equal probability of a pulse at any sample (so for example the number of pulses in an interval follows the Poisson distribution). Further suppose the strength of each pulse is equal.

Can the coefficients of the FIR filter be recovered from the filtered signal, and how?

Now assume the noise is impulsive, but may be somewhat correlated with itself (it doesn't have perfect Poisson statistics). Also assume that the strength of each pulse may vary. Can the coefficients be recovered approximately if the exact distribution of noise is not known?

Noise source $N[i]$ where $i=0..n, N[i] = 1$ with some unknown but low probability $p$, 0 otherwise; each $N[i]$ is independent of other $N[i]$

Filter $F[k]$ where $k=0..m$, unknown, $m$ is much less than $n$

Signal $S[i]$ where $i=0..n$, $S = N * F$, where * is convolution

Given $S$, estimate $F$. $N$ and $p$ are also unknown.

To give this some practical background, the average period between noise pulses can be 20-40 samples, the FIR filter can have a few hundred nonzero coefficients (so the filtered signals resulting from each pulse overlap significantly), and the total signal is a few thousand samples.

EDIT: Some real-world data here: http://pastebin.com/N3LZU6vm  I think what you see at sample 780-810 is pretty close to the impulse response to a single impulse.  This is generated by noise fed into a complex resonator, whether the noise is impulsive as described above is not actually known but is physically plausible.

Some synthetic data here: http://pastebin.com/qTYLP5eA  This is generated from the following python code:

noise = numpy.where( numpy.random.rand(5000) &lt; 1/20.0, numpy.ones(5000), numpy.zeros(5000))
filter = numpy.sin( numpy.arange(200) * 2 * math.pi / 11 ) * numpy.exp( - numpy.arange(200) / 30.0 )
signal = numpy.convolve( noise, filter )
for i in range(len(signal)): print int(signal[i] * (2**15)/ 3.0)


This creates a filter wavelet from a sinewave with period 11 samples and an exponentially decaying envelope, and noise where each sample has an independent 5% probability of being an impulse.

EDIT:  I think what I am looking for is very similar to what is involved in building an autoregressive model (AR).  There are slight differences in the formulation, for example AR is based on an IIR with lagged signal plus a noise term rather than FIR with lagged noise, and the noise is assumed to be Gaussian rather than impulses with Poisson distribution.  However I think AR proves that estimating a filter from filtered noise is possible, and in fact AR may just work for my problem.
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