[Math] How to perform deterministic deconvolution

Question

Here is my problem:

I have a random variable $A$ that is the sum of independent random variables $B$ and $C$, i.e. $A=B+C$. All three random variables are on the real domain. $B$ is a Gaussian with known parameters and I have the form of the p.d.f. $f_A$ for $A$ with good estimate of the parameters. $A$ is symmetric around its mean, which is equal to the mean of Gaussian $B$. Assume that variances of $A$ and $B$ (and $C$) are low enough though that their p.d.f.'s are very close to zero outside the interval $[\mu-a,\mu+a]$, where $a$ is reasonable.

Unfortunately, there is no nice closed-form expression for the characteristic function of $A$. I am trying to find the p.d.f. of $C$ (though having characteristic function $\phi_C$ of $C$ would also help.)

If I had a nice form for $A$'s characteristic function $\phi_A$, I would use the fact that $f_A=f_B*f_C$, which in turn means $\phi_A=\phi_B\phi_C$, and would simply divide $\phi_A$ by characteristic function of the Gaussian (since convolution in "time" domain is multiplication in "frequency" domain) to at least get $\phi_C$ to work with.

At this point I would be happy with a numerical approximation for p.d.f. $f_C$. Poking around the web, I found that one may be able to use "deterministic deconvolution" to remove the known Gaussian p.d.f. from the p.d.f. of $A$, but I couldn't find any references on how to actually do it (there a lot of literature on the use of deconvolution in seismography though.) I tried generating a vectors $f_A(\mathbf{x})$ and $f_B(\mathbf{x})$ where $\mathbf{x}$ is a large ordered list of equally-spaced numbers from $\mu-a$ to $\mu+a$, taking FFTs of both vectors, doing element-by-element division and then taking the inverse FFT of the result, which yielded garbage. Can somebody help?

Answer 1

Assuming you know exactly $f_B$ and that you measure exactly $f_A$, it's true that $f_C$ is in principle obtainable from $\phi_C= \phi_A/\phi_B$. The problem is, obviously, division by zero; this teorically prevents us to obtaining $\phi_C$ where $\phi_B=0$, and practically also in the regions where $\phi_B \approx 0$ (numerical instabilities). If we know that the "error" ($B$) is "small", that means that it has little variance (with respect to C) and hence its Fourier transform is wideband - this might alliviate the problem of regions with $\phi_B=0$; and because its a gaussian, we don't need to worry about zero crossings; but numerical noise (measure and computation) complicates things.

No foolproof recipes here, but some experiment to get the feeling - and at least reassure us that we are doing the fft right. Because you know the means, I'll assume that we have zero mean; that, together with the pdf being symmetric, means that the Fourier transform is real (nothing important, but it gets a little easier to graph).

In Octave:

N=100; % point in frecuency
dw=6/N; % delta freq.
w=[-3:dw:3-dw]; % frequencies

phib=e.^(- w.^2); % CF of a zero mean gaussian is a gaussian
phic=sinc(w*1.8); % suppose this is the CF of C
phia=phib .* phic; % this is the exact CF of A
plot([phia',phib',phic']); legend('CF(A)','CF(B)','CF(C)');

% compute the "exact" (except for roundin/discretization errors) densities
pa=abs(fftshift(fft(phia)))/N;  % the shift gets the center where it should be
pb=abs(fftshift(fft(phib)))/N;
pc=abs(fftshift(fft(phic)))/N;
plot([pa',pb',pc']);       legend('p(A)','p(B)','p(C)');

Lets now recover C from A,B and, to see the problem, let's add a little noise to both.

phian=phia+(rand(1,N)-0.5)*0.01; % noisy version of CF(A)
phibn=phib+(rand(1,N)-0.5)*0.01; % noisy version of CF(B)
phicn = phian ./ phibn;          % naive deconv
phicn(phicn>1)=1;                % first naive correction: truncate the CF below 1
phicn(phicn<-1)=-1;              % ... and over -1
% lets plot the noisy CF 
plot([phian',phibn',phicn']);  legend('CFn(A)','CFn(B)','CFn(C)');
% ouch... and the densities
pan=abs(fftshift(fft(phian)))/N;
pbn=abs(fftshift(fft(phibn)))/N;
pcn=abs(fftshift(fft(phicn)))/N;
plot([pan',pbn',pcn']);       legend('pn(A)','pn(B)','pn(C)')

Because I considered relatively high frequencies (in relation with the gaussian) there are some clear numerical problems. To alliviate it, a first, rather brutal -but not very bad- recipe: look at the frequencies where $\phi_B$ is below some threshold, and consider $\phi_C$ "unknown" for those values: for example, set it to zero.

phicn2(abs(phibn)<0.02)=0;
plot([phian',phibn',phicn2']); legend('CFn(A)','CFn(B)','CFn2(C)')

pcn2=abs(fftshift(fft(phicn2)))/N;
plot([pan',pbn',pcn2']); legend('pn(A)','pn(B)','pn2(C)')

That's better. From this, you could get more ideas. In this case, because we know $\phi_B$ is gaussian, I'd restrict the frequencies to the range where its values exceed some threshold.