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

• This is done all the time in electrophysiology by a technique that is usually called reverse correlation (although it's a bit of a misnomer). Whether you can reconstruct F (the kernel) depends on the characteristic of N. What is the support of F? I can work out the details for you another time. No time now. In fact, this looks like a physiology problem. One cell fires randomly following Poisson. It synapses into another cell. You're recording from both cells. You're looking for the EPSP/IPSP kernel. – thang Jan 23 '13 at 19:49
• @thang: thank you, fascinating! and no, this is from a different field. From looking at the cases where pulses are far apart, support of F is 100-200 samples (a few times longer than average spacing between pulses), the envelope of F seems to be approximately exponentially decaying with time constant ~50 samples. – Alex I Jan 23 '13 at 20:17
• You have an interesting question. I need to think about it some more and try a couple things out. Do you happen to have any sample datasets that are representative? – Jason R Jan 24 '13 at 13:09
• @JasonR: I added some datasets. Let me know if you need more. Thanks! P.S. Bounty coming up :) – Alex I Jan 24 '13 at 20:20
• can you post whatever statistical properties of n that you know? if you know it's Poisson, post the estimated parameters. otherwise, post sample n. is it generated by an ergodic source? what is the mean separation between the pulses? or what is the separation distribution (exponential?)? – thang Jan 25 '13 at 2:05

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.

function [Fhat,q] = blindDeconvImpulseNoise(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
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

• If you can provide links to the images, then someone else can edit your post to add them in. I'm interested to see your results. – Jason R Jan 27 '13 at 18:59
• Sadly they thought of that too: new users can't link to things either – Gruntled Jan 27 '13 at 20:17
• @Gruntled: This is amazing! Probably the most epic stackexchange answer I have ever seen. Thank you. I'll have to read up on the algorithm and get your code running, but I think this does the job nicely. – Alex I Jan 27 '13 at 23:21
• Great to see someone could pull together a suitable answer. @AlexI I would love to know how this goes, so a follow up post would be greatly appreciated. I have a kind of similar issue I am working on. In my case, I'm trying to locate the time separation between impulses very precisely. – user2718 Jan 28 '13 at 13:35
• This is phenomenal! Ok @AlexI, I owe you some rep points :p I am very surprised that EM can converge with so little information given that there are so many solutions that would generate the same data points. Remark then that you may be able to impose additional regularization constraints by using MAP - en.wikipedia.org/wiki/Maximum_a_posteriori_estimation – thang Jan 28 '13 at 15:14

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

• thank you, this is interesting. the problem is that N is unknown also... N is a random signal, one can to some extent guess where the nonzero elements (pulses) are from the highest peaks in S but some of those can be caused from superposition of the impulse responses from several earlier peaks as well. the solution you have is great for known N and unknown F, but what we have here is independeltly random unknown N as well as unknown F. – Alex I Jan 24 '13 at 7:38
• then you can't recover anything. more unknown than variables. even if n were known, you can't guarantee to be able to recover f all the time. if you don't known n, then all bets are off. – thang Jan 24 '13 at 9:52

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

• the paper you quoted is absolutely unrelated as you pointed out. it is actually on extending the sampling theorem to irregular sampling. that aside, there is also a set of other papers on system identification using random input, but they all assume you know the input, and they control it to make well behaving Toeplitz matrix. 1 point out of every 30 samples means the signal has a 30 sample period. this actually represents a low freq limit on recoverability of your filter and actually is reflected by rank of the Toeplitz of this signal. – thang Jan 24 '13 at 15:28
• The general problem is to recover a signal from random views of that signal. I would think that the research on system ID from random inputs would also be relevant here. Why not post a link to an example or two that at least loosely fit the problem at hand. – user2718 Jan 24 '13 at 15:57
• Did you read the paper? it is to reconstruct the analog signal from random sampling of the signal. standard sampling is periodic. the research on system ID all end up with equation (2). they control the input to make well behaving matrix N. i suspect that OP has no control over the input. otherwise he would know what it is, which he doesn't. this is where the hangup is. otherwise this problem is just standard system id. standard system id assumes you know n and s (sometimes have control over n). what i posted is actually from standard id. in this problem, you don't know n. – thang Jan 24 '13 at 15:59
• Yes I did read that paper and several others on recovering signals from random samples. I well understand where the hang up is. You won't find a solution using deterministic methods. There may not be a ready made textbook solution that you can apply to this problem. Probably why it was posted in the first place. Since no one has offered up a solution, I am simply making a suggestion to broaden your horizons. It could spark an idea. – user2718 Jan 24 '13 at 16:15
• I see. I just wanted to remark that his issue is not from sampling. Random sampling has nothing to do with it, although it is a cool idea, and has a lot of implications in other areas. The specific problem that is presented here is not solvable and provably so. Basically, you're asked to find the solution to an under-ranked equation. It's mathematically not solvable. Usually what this means is that there is additional information from context (specific application) that OP is not disclosing. – thang Jan 24 '13 at 16:21