reconstruction of distribution function, deconvolution

Hi,

I'm working in the field of chemical engineering and I'd like to reconstruct processes wich are distorted by the residence time distribution of the reactor. I simplified the problem in order to develope a strategy to work out the solution: First I generate two gaussian functions A and B. The function A and B are convoluted to obtain the distorted function C. Assuming function B and C are given I'd like to reconstruct function A. I tried to use [q,r] = deconvol(B,C) but it's not working ("First coefficient of A must be non-zero."). My idea now is to predict a arbitrary gauss function and to iteratively change the parameter and convolute with B until the resulting function (D) fits well with the function C. But I have no clue how to to it!

Below you'll find the code I have wrtitten so far. Who knows a way to reconstruct the initally given function A??? I would be very pleased to receive your hints or suggestions!

Nico

%%test convolution / deconvolution
clear all; close all;
x = linspace(0,10000,10000);

meanA = 250;                        % expectation value of curve A
sigmaA = 5;                         % standard deviation of curve A

meanB = 450;                        % expectation value of curve B
sigmaB = 6;                         % standard deviation of curve A

gaussA = gauss(x,meanA,sigmaA);    % gauss curve A
gaussB = gauss(x,meanB,sigmaB);    % gauss curve B

gaussC = conv(gaussA,gaussB);      % convolution of curve A and B         
gives distorted curve C
gaussC = gaussC(1:size(x,2));       % resize gaussC to inital points of data

figure(1);
plot(x,gaussA,'b',x,gaussB,'r',x,gaussC,'m'); 
legend('gauss (A)', 'gauss (B)', 'distorted gauss (C)'); 
xlim ([0 1000]); ylim([0 3]);

Answer 1

[Matt J]

You could use the equations

 meanC=meanA+meanB
 varC=varA+varB

to solve for meanA and varA, which determine A completely.

[Nicolai Zientek]

Thanks a lot for your answer. This will solve the problem discribed above. But I'd like to apply the problem to different kind of distribution functions e.g. logistic and gauss distributions. So it would be good to find a general solution.

So, how can I iteratively change the parameter of a gauss function, do a convolution and optimeze the parameter that the resulting function fits the distorted gauss function?

Thanks for further discussions, hints and suggestions!

Nico

Answer 2

[Matt J]

You could try fminsearch, or if you have the Optimization Toolbox, there's FMINUNC or FSOLVE.

[Nicolai Zientek]

Hi Matt,

thank you very much for your suggestions. I think this is the right way to solve the problem. I already used the excel solver wich works quite similiar but now I want to implement it into my Matlab code. Unfortunatly I'm quite new in working with Matlab. I checked the help for fminsearch but I don't know how to apply it to my problem. Could you please outline the approach? Thank you very much.

Nico

[Matt J]

Well, you've written most of it already. Basically, you need a function which measures the difference between the given gaussC and a hypothetical one generated according to the code you've posted. Then minimize this function using FMINSEARCH:

 initialguess=[250;5];
 estimate=fminsearch(@(params) fun(params,gaussC_given),initialguess);

 function Error=fun(parameters,gaussC_given) 
  x = linspace(0,10000,10000);
  meanA = parameters(1);                        % expectation value of curve A
  sigmaA = parameters(2);                         % standard deviation of curve A
  meanB = 450;                        % expectation value of curve B
  sigmaB = 6;                         % standard deviation of curve A
  gaussA = gauss(x,meanA,sigmaA);    % gauss curve A
  gaussB = gauss(x,meanB,sigmaB);    % gauss curve B
  gaussC = conv(gaussA,gaussB);      % convolution of curve A and B         
  gives distorted curve C
  gaussC_hypothetical = gaussC(1:size(x,2));       % resize gaussC to inital points of data

 Error=norm(gaussC_hypothetical - gaussC_given)^2;

end

[Matt J]

Rather than generating B from scratch each time you call fun(), you could also pre-compute gaussB and pass that to fun() as fixed data, just as we did above with gaussC_given:

 estimate=fminsearch(@(params) fun(params,gaussC_given, gaussB_given),initialguess);

[Matt J]

Nicolai Zientek commented:

Hello!

I have implemented the error function into the code but I still don't get how to optimize the parameter via minimization of the error value with fminsearch... Would be nice if a experienced user could give me a hint. Thank you very much!

This is the code I have written so far, I highlited the part where the parameters should be optimized...

Further down you'll find the code of the function that calculates the error value.

%%test convolution / deconvolution
clear all; close all;
x = linspace(0,10000,10000);        % x values 
meanA = 250;                        % expectation value of curve A
sigmaA = 5;                         % standard deviation of curve A
meanB = 450;                        % expectation value of curve B
sigmaB = 6;                         % standard deviation of curve B
gaussA = gauss(x,meanA,sigmaA);    % gauss curve A
gaussB = gauss(x,meanB,sigmaB);    % gauss curve B
gaussC = conv(gaussA,gaussB);       % gauss cuve C = convolution of gauss   curve A and B   
gaussC = gaussC(1:size(x,2));       % resize gaussC to inital points of data
%%plot data
figure(1);
plot(x,gaussA,'b',x,gaussB,'r',x,gaussC,'m'); 
legend('gauss (A)', 'gauss (B)', 'distorted gauss (C)'); 
xlim ([0 1000]); ylim([0 3]);
%%optimization
meanAn = 240;                       % initial guess for new gauss curve A
sigmaAn = 4.5;                        % initial guess for new gauss curve A
Error = GaussErr(x,meanAn,sigmaAn,gaussB,gaussC);
%{
hear comes the part where the parameters are optimized
%}
gauss1n = gauss(x,meanAn,sigmaAn);  % new gauss curve A with optimized   parameter
gauss3n = conv(gauss1n,gaussB);     % new gauss curve C via convolution of ne gauss A with gauss B
gauss3n  = gauss3n(1:size(x,2));    % resize gauss3n to inital points of data 
hold on
plot(x,gauss1n,'k.',x,gauss3n,'g.');

And here comes the error function:

function err = GaussErr(x,mean,sigma,gauss2,gauss3);
pred = gauss(x,mean,sigma);           % predicted gauss curve
pred = conv(pred,gauss2);             % convolution 
pred = pred(1:size(x,2));             % rezise to initial points of data
err = sum((pred(:)-gauss3(:)).^2);    % error value
end

[Matt J]

z0 = [mean0; sigma0]; %initial guess

meanANDsigma=fminsearch(@(z) GaussErr(x,z(1),z(2),gauss2,gauss3) , z0);

Answer 3

[Image Analyst]

Sounds like just a simple inverse filter would do the job. If you have the Image Processing Toolbox, type wiener into the help and see what pops up. The Wiener filter is basically an inverse filter taking into account the noise.

[Matt J]

That may recover the signal A, but you would still have to do some sort of post-fitting to recover its parameters meanA and sigmaA. It would also not take advantage of that parameterization to improve the accuracy of the reconstruction.

It would, however, be a good way of generating an initial guess.

[Nicolai Zientek]

Good idea, but unfortunatly I don't have the Image Processing Toolbox!

Answer 4

[Nicolai Zientek]

Hello!

I have implemented the error function into the code but I still don't get how to optimize the parameter via minimization of the error value with fminsearch... Would be nice if a experienced user could give me a hint. Thank you very much!

This is the code I have written so far, I highlited the part where the parameters should be optimized...

Further down you'll find the code of the function that calculates the error value.

%%test convolution / deconvolution
clear all; close all;
x = linspace(0,10000,10000);        % x values 

meanA = 250;                        % expectation value of curve A
sigmaA = 5;                         % standard deviation of curve A

meanB = 450;                        % expectation value of curve B
sigmaB = 6;                         % standard deviation of curve B

gaussA = gauss(x,meanA,sigmaA);    % gauss curve A
gaussB = gauss(x,meanB,sigmaB);    % gauss curve B

gaussC = conv(gaussA,gaussB);       % gauss cuve C = convolution of gauss   curve A and B   
gaussC = gaussC(1:size(x,2));       % resize gaussC to inital points of data

%%plot data

figure(1);
plot(x,gaussA,'b',x,gaussB,'r',x,gaussC,'m'); 
legend('gauss (A)', 'gauss (B)', 'distorted gauss (C)'); 
xlim ([0 1000]); ylim([0 3]);

%%optimization

meanAn = 240;                       % initial guess for new gauss curve A
sigmaAn = 4.5;                        % initial guess for new gauss curve A

Error = GaussErr(x,meanAn,sigmaAn,gaussB,gaussC);

%{
hear comes the part where the parameters are optimized
%}

gauss1n = gauss(x,meanAn,sigmaAn);  % new gauss curve A with optimized   parameter
gauss3n = conv(gauss1n,gaussB);     % new gauss curve C via convolution of ne gauss A with gauss B
gauss3n  = gauss3n(1:size(x,2));    % resize gauss3n to inital points of data 

hold on
plot(x,gauss1n,'k.',x,gauss3n,'g.');

And here comes the error function:

function err = GaussErr(x,mean,sigma,gauss2,gauss3);

pred = gauss(x,mean,sigma);           % predicted gauss curve
pred = conv(pred,gauss2);             % convolution 
pred = pred(1:size(x,2));             % rezise to initial points of data
err = sum((pred(:)-gauss3(:)).^2);    % error value
end

[Matt J]

This seemed like it should be a Comment to my Answer, so I copied it (and responded to it) there.