Deconvolution of probability density functions via FFT

I am trying to implement a code that yields the probability density function , given knowledge of two other PDFs, and . The whole procedure is based on the paper Pernpeintner et al., A Method to Obtain Planar Mixture Fraction Statistics in Turbulent Flows Seeded With Tracer Particles (2011), but I am having trouble getting any meaningful results.I was hoping someone here could maybe point out what my mistake is.

The random variables f and ω are statistically independent and have the relationship .

The PDFs supposedly have the relationship:

There is a numerical method to invert this, but I thought I would try the analytical method first to get some verifiable results. For the analytical method one makes the substitution for all random variables. This yields a convolution integral:

Which should be solvable via Fourier transformation: and the relationship .

I attempted to implement this in Matlab and to verify the method with synthetic PDFs for which I know the answer. I generate as a Gaussian and as a Beta distribution and compute with the first equation above. The variable space for all of the random variables is 0-1 with 256 samples.

Then I change variables, use an FFT on and , divide it and perform the iFFT to get back. However, the resulting signal is nowhere near the original Gaussian (after changing variables back), see this figure:

Below I also added the Matlab code. I would appreciate it a lot if somebody could explain to me, what I am doing wrong or maybe what I conceptually do not get correct.

% Evaluation points 8-bit

centers = 0.5:1:255.5;

centers_norm = centers/256;

% Seeding PDF (omega)

seedPDF = pdf('Beta',centers_norm,5,2);

seedPDF = seedPDF./trapz(centers_norm,seedPDF);

% mix PDF (f)

mixPDF = pdf('Normal',centers,128,20);

mixPDF = mixPDF./trapz(centers,mixPDF)*2^8;

% Convolute seeding and mixing pdf to get intensity measurement (theta)

for i = 1:length(mixPDF)

mix_interp = interp1(centers_norm,mixPDF,centers_norm(i)./centers_norm,'linear');

mix_interp(isnan(mix_interp)) = 0;

miePDF(i) = trapz(centers_norm,seedPDF./centers_norm.*mix_interp);

end

% low-pass filter pdf

windowWidth = 5;

kernel = ones(windowWidth,1)/windowWidth;

seedPDF_lowp = filter(kernel,1,seedPDF);

miePDF_lowp = filter(kernel,1,miePDF);

% seedPDF_lowp = seedPDF;

% miePDF_lowp = miePDF;

% transform to hat quantities

variable_hat = log(centers);

seedPDF_hat = seedPDF_lowp.*centers_norm;

% seedPDF_hat = [seedPDF_hat fliplr(seedPDF_hat)];

miePDF_hat = miePDF_lowp.*centers_norm;

% miePDF_hat = [miePDF_hat fliplr(miePDF_hat)];

% Perform Fourier transformation and deconvolute

n = 2^(nextpow2(length(seedPDF_hat)));

W = 1;

seedPDF_FT = fft(W'.*seedPDF_hat,n);

miePDF_FT = fft(W'.*miePDF_hat,n);

testPDF_FT = fft(W'.*mixPDF.*centers_norm,n);

mixPDF_FT = miePDF_FT./(seedPDF_FT);

mixPDF_hat = ifft(mixPDF_FT,n);

mixPDF_analytic = mixPDF_hat./centers_norm;

plot(centers_norm,miePDF,'b')

hold on

plot(centers_norm,miePDF_lowp,'b-.')

plot(centers_norm,mixPDF,'r')

plot(centers_norm,seedPDF,'g')

plot(centers_norm,mixPDF_analytic,'r--')

trapz(centers_norm,centers_norm.*mixPDF)

trapz(centers_norm,centers_norm.*mixPDF_analytic)

[Aiqin Jiang]

Hi Dan, did you find the answer? I have similar question like you have here. if you find the answer, would you please share? thanks! Aiqin

[Dan Fries]

Not really, I got the numerical method to work very well, but the analytical approach is still an enigma. It might have something to do with the the Fourier transformation's sensitivity to noise or my inability to do the transformation into logarithmic space correctly.