Problem using ifft and fft

Dear all,

I am facing a problem in the use of the fft and ifft. What i am try to do is not particularli complicated. I would like to extract the vector g coming from the deconvolution of two exponential functions namely f and h where

f = exp(-exp(x-log(tau0))*beta) ; h = exp(-exp(x));

The idea is to deconvolute them using the properties of the Fourier transform since I cannot find an easier deconvolution algorithm. However when i apply this the results are complitely wrong.

clear all


beta    = 0.2;
tau0    = 1;

%%%  definition xaxis log scale x = ln(t)
x = (-15:0.01:10);

%%% definition pure exponential decay and stretched exponential decay
f  = exp(-exp((x-log(tau0))*beta));
h  = exp(-exp(x));
%%% test with derivatives
df = diff(f);
dh = diff(h);

%%% deconvolution of f and h to find g

F = fftshift(fft(f));
H = fftshift(fft(h));

G = ifft(ifftshift(F./H));

figure(1)
    yyaxis left
    plot(x,h);    
    ylabel('h(x)')
    hold on
    yyaxis right
    ylabel('f(x)')
    plot(x,f);

    
    

Answer 1

[Matt J]

Make sure you are familiar with the concepts here,

https://www.mathworks.com/help/signal/ug/linear-and-circular-convolution.html

[Nicola Modolo]

Dear Matt,

Thank you for the suggestion. After reading the link I anderstood some of my mistakes and i updated the code as below. The test using two gaussian curves works fine. both for the convolution and deconvolution. However, When i try to deconvolute f and h the result is quite weird in particular referring to the first point of the resulting vector g. Am I missing something? thank you very much for your help and time.

Cordially,

Nicola

clear all


beta    = 0.3;
tau0    = 1;

%%%  definition x axis
step = 0.001;
x = (-30:step:20);
xconv = (-30:step/2:20);

%%%  test to convolute 2 gaussian functions
gauss1 = gaussmf(x,[1 0]);
gauss2 = gaussmf(x,[1 -5]);

%%%  definition of the vectors before fft and fft
sig1 = [gauss1 zeros(1,length(gauss2)-1)];
GAUSS1 = fft(sig1)*step;
sig2 = [gauss2 zeros(1,length(gauss1)-1)];
GAUSS2 = fft(sig2)*step;

%%%  convolution
CONV   = GAUSS2.*GAUSS1;
conv   = ifft(CONV)/step;

%%%  deconvolution test
DECONV = CONV./GAUSS1;
deconv = ifft(DECONV)/step;
deconv = deconv(1:length(gauss1));

%%% definition of pure exponential and strecthed exponential
f  = exp(-exp((xconv-log(tau0))*beta));
h  = exp(-exp(x));

%%% test derivatives
df = diff(f);
dh = diff(h);

%%%  definition of the vectors before fft and fft
lf = [f zeros(1, length(h)-1)];
lh = [h zeros(1, length(f)-1)];
F = fft(lf)*step;
H = fft(lh)*step;

%%% deconvolution of F and H and resize of g
g = ifft(F./H)/step;
g = g(1:length(h));

figure(1)
    yyaxis left
    plot(x,gauss1,x,gauss2,x,deconv);    
    ylabel('h(x)')
    hold on
    yyaxis right
    ylabel('f(x)')
    plot(xconv,conv);
    
figure(2)
    yyaxis left
    plot(xconv,f,x,h);    
    ylabel('h(x)')
    hold on
    yyaxis right
    ylabel('f(x)')
    plot(x,g);    

[Matt J]

Try this,

beta    = 0.3;
tau0    = 1;
step = 0.001;
N=1e5;

x = step*( (0:N-1) -ceil((N-1)/2)  );
f  = exp(-exp((x-log(tau0))*beta));
h  = exp(-exp(x));

%%% Numerical differentiation (right hand derivatives)
df = [diff(f),0]/step;
dh = [diff(h),0]/step;

%%Fourier transformation of the derivatives
dF = fft(ifftshift(df))*step;
dH = fft(ifftshift(dh))*step;

%%% deconvolution of dF and dH
g = ifft(divSpectrum(dF,dH),'symmetric')/step;
g = fftshift(g);

plot(x, conv(dh,g,'same')*step, x(1:500:end),df(1:500:end),'x');
legend('conv(dh,g)', 'df')

function Q=divSpectrum(dF,dH)

     tol=1e-5;

     supp=abs(dF)>tol & abs(dH)>tol;

     Q=dF./dH;

     Q(~supp)=0;


end