Problem using ifft and fft

2 次查看（过去 30 天）
Nicola Modolo
Nicola Modolo 2020-12-15
编辑： Matt J 2020-12-17
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);

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2020-12-15
编辑：Matt J 2020-12-15
Make sure you are familiar with the concepts here,
https://www.mathworks.com/help/signal/ug/linear-and-circular-convolution.html
  2 个评论
Nicola Modolo
Nicola Modolo 2020-12-16
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
Matt J 2020-12-17
编辑：Matt J 2020-12-17
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Special Functions 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by