Recover Filter Coefficients from Filtered Noise

Question

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) < 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.

Answer 1

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>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) < 1e-10
     % No change on this iteration: stop here
     break
  end
  last_sig = sigmasq_hat;

  if 24*60*60*(now - starttime) > 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

Answer 2

Ok, so for the discrete case, this is actually pretty simple. It comes out to a big matrix equation.

Basically, you have that:

(1) $s = n * f$

(notice that the conventional way to write this is $s = n*f$ but you can just as well transpose as convolution is commutative)

as you have stated. In fact, this is really simple :p I have used lower case letters to represent your variables. This is because they are sequences, represented as vectors. I am reserving uppercase for matrices.

Define $N$ to be (the transpose of) the Toeplitz matrix of $n$ given by:

$ N = \left[ \begin{array}{ccc} n_1 & 0 & . & 0 \\ n_2 & n_1 & . & 0 \\ . & . & . & . \\ n_m & n_{m-1} & . & n_1 \\ 0 & n_m & . & n_2 \\ . & . & . & . \\ 0 & 0 & 0 & n_m \\ \end{array} \right] $

So equation (1) becomes

(2) $s = Nf$

just a matrix multiplication. In this equation, you know $s$ and $N$ and want to find $f$.

You probably have a lot of input and output samples compared to the size of the filter, so the length of $s$ is a lot bigger than $f$.

From linear algebra, recall that if rank(N)>length(f), then you can completely reconstruct $f$ from $n$ and $s$. Otherwise, you can't. In this case, you can shuffle things around, but there's always going to be missing information.

So how do you solve for $f$?

It's just a simple over-determined matrix equation.

One way is to use least square by evaluating $ argmin_f||s-Nf||$. There are many ways to do this, and it is the topic of several chapters in numerical methods. I won't go into that.

EDIT:

as per discussion below, it seems I misread and n is also unkown. In this case it is impossible to recover f. This is a case of system identification. See here Get an input signal from LTI system output, find impulse response for a brief discussion about it.

This issue is really quite simple. I say c = axb (multiply, not convolution) and then tell you that c is 10. what is b? There is no way to tell.

Your specific scenario is actually worse than that. It is difficult to even make a connection between the distribution of s and n. This is primarily because output s generated by convolving n with f converts most distributions into something that looks normal (recall the central limit theorem), assuming that n is generated by an iid process.

If you know that $n$ is sparse and can guarantee that the pulses are always farther apart than the support of f, then just model this into the formula.

And in fact, you don't have to do anything. whenever $s>0$, the portion of $s$ in the vicinity that is $>0$ a shifted version of your $f$. Of course, there is no way to know the degree of shift.

If you know that there is a small probability of overlap between the pulses, then what you can do is generate $f_1$, ... $f_p$ for each time $s>0$. Most of them will look similar (plus noise), and the rest will look like a shifted sum of the majority. What you can do, then is, out of the $f_j$, take the ones with smallest support, shift them so the $L_2$ norm is minimal, and average the shifted versions them together. This gives you an estimated $f$ (within shift in time) from pulses that are far apart. You can similarly model 2, 3, 4, etc. overlapping pulses...

Answer 3

Here are a couple of ideas (not solutions)

First, imagine what would happen if you're noise was instead a periodic impulse at the average period between samples (lets say a pulse every 30 samples). You would end up with a periodic sequence. Can you try averaging your data, perhaps using 30 sample long windows in such a way that the result, given sufficient sample size approaches this periodic response. Once you have this "average" response signal, find the DFT over one period of the signal. Then interpolate the DFT points to prvide enough detail (points/bins) to exceed the estimated support of F (say 200 bins). Then take the IDFT of the result. The reslult will provide an approximation of the FIR "filter" coefficients.

Second you can look into research on the representation and recovery of signals using random sampling. This is a different problem, but the techniques may be useful for the problem you are trying to solve. Here is a link to an article: http://deepblue.lib.umich.edu/bitstream/handle/2027.42/22215/0000648.pdf;jsessionid=92ABBC582E3D6909A24775F88D6543D2?sequence=1