ifft returns NaN when plotting the impulse response function

23 visualizaciones (últimos 30 días)
粤轩 杨
粤轩 杨 el 10 de Abr. de 2023
Comentada: 粤轩 杨 el 10 de Abr. de 2023
Hello, I seem to be having issues using MATLAB's fft and ifft functions.
An impulse response function h(t) has the following formula: inj(t) * h(t) = AIF(t). We know that the graphs of inj(t) and AIF(t) are as followed. I wrote the following code to do the deconvolution but h(t) in the output graph is zero.
I realized that the h returns from ifft is NaN but I don't know how to correct the code. Thanks!
load("AIF_1.mat");
load("inj_1.mat");
inj_1 = inj(201:2400);
inj1_FFT = fft(inj_1);
AIF_1 = AIF(201:2400);
AIF1_FFT = fft(AIF_1);
h_FFT = AIF1_FFT ./ inj1_FFT;
h_FFT(isnan(h_FFT)==1) = 0;
h = ifft(h_FFT);
X = zeros(1,200)
ht = [X,h];
plot(time,real(ht));title('h(t)');

Respuesta aceptada

Paul
Paul el 10 de Abr. de 2023
Hi 粤轩 杨,
It looks like h_FFT also has a few values that are inf, in addition to the NaNs
load("AIF_1.mat");
load("inj_1.mat");
inj_1 = inj(201:2400);
inj1_FFT = fft(inj_1);
AIF_1 = AIF(201:2400);
AIF1_FFT = fft(AIF_1);
h_FFT = AIF1_FFT ./ inj1_FFT;
sum(isnan(h_FFT))
ans = 3
sum(isinf(h_FFT))
ans = 2
% quick correction to get things to run, not sure if this is really what
% should be done
inj1_FFT(inj1_FFT == 0) = eps;
h_FFT = AIF1_FFT ./ inj1_FFT;
sum(isnan(h_FFT))
ans = 0
sum(isinf(h_FFT))
ans = 0
h = ifft(h_FFT);
X = zeros(1,200);
ht = [X,h];
plot(time,real(ht));title('h(t)');
  1 comentario
粤轩 杨
粤轩 杨 el 10 de Abr. de 2023
Hi Paul,
The output graph is a little strange. But you finging of inf values in h_FFT made me realize that I didn't intercept the inj correctly.
After changing to
inj_1 = inj(200:2400)
the output graph is good!
Thank you for your inspiration!

Iniciar sesión para comentar.

Más respuestas (1)

粤轩 杨
粤轩 杨 el 10 de Abr. de 2023
There are two related attachments for your information. Thanks!
  1 comentario
粤轩 杨
粤轩 杨 el 10 de Abr. de 2023
Just as mentioned above, the final code should be
load("AIF_1.mat");
load("inj_1.mat");
inj_1 = inj(200:2400);
inj1_FFT = fft(inj_1);
AIF_1 = AIF(200:2400);
AIF1_FFT = fft(AIF_1);
h_FFT = AIF1_FFT ./ inj1_FFT;
h = ifft(h_FFT);
X = zeros(1,199);
ht = [X,h];
plot(time,real(ht));title('h(t)');
And the output graph looks like

Iniciar sesión para comentar.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by