Asked by Keith Hernandez
on 1 Nov 2018

I am trying to deconvolve two signals using fft(). The issue I'm having is with the pulse width. If the pulse width is an even number of ones with zero padding, the fft() yields a zero in the resulting array. Therefore, when I try to divide my Y_fft by H_fft, I get NaN as the result.

The code below zeros pads a pulse to some power of 2 because my signal I want deconvolve will be forced to this length. I have played around with finding the zeros and replacing them 1e-15 or replacing them with the average of the surrounding values. Both are quick fixes with no real math to prove why to do this. Therefore, my deconvolved signal is reconstructed but with some jitter. I would greatly appreciate feedback on the correct way to handle this. At the very least, some understand as to why even number of pulses yield zeros after the fft().

Note: The position of the zeros seem to also be at location of powers of 2 and shifted.

if true

clc

clear

close all

%original signal

x = [1 2 3 4 5 6 5 4 3 2 1];

lx = length(x);

%odd pulse case

pulses_odd = 3;

h1 = ones(1,pulses_odd);

%creation of convolved signal

y1 = conv(x,h1);

ly1 = length(y1);

y1 = [y1 zeros(1, pow2(nextpow2(ly1)) - ly1)];

ly1 = length(y1);

%zero padding impulse function

h1 = [ones(1,pulses_odd) zeros(1,ly1-pulses_odd)];

Y_fft1 = fft(y1);

H_fft1 = fft(h1);

x1_decon = Y_fft1./H_fft1;

x1_decon = ifft(x1_decon);

%even pulse case

pulses_even = 4;

h2 = ones(1,pulses_even);

%creation of convolved signal

y2 = conv(x,h2);

ly2 = length(y2);

y2 = [y2 zeros(1, pow2(nextpow2(ly2)) - ly2)];

ly2 = length(y2);

%zero padding impulse function

h2 = [ones(1,pulses_even) zeros(1,ly2-pulses_even)];

Y_fft2 = fft(y2);

H_fft2 = fft(h2);

x2_decon = Y_fft2./H_fft2;

x2_decon = ifft(x2_decon);

plot(x,'linewidth',2)

hold on

plot(x1_decon,'d')

plot(x2_decon,'s')

hold off

end

Answer by David Goodmanson
on 1 Nov 2018

Accepted Answer

Hi Keith,

This is happening because the array length of 16 (after you zerofill) and a pulse length of 4 have a factor in common. For an array length of 16, a wave of period 4 fits (4 periods give 16) so that wave can contribute to the fft.

For the constant-height pulse of length 4, a wave of period 4 exactly averages out over the width of the pulse and the fft amplitude due to the nonzero part of the pulse is zero. It doesn't matter what is going on in the 12 points off the end of the pulse because the pulse height is zero there and can't contribute to the fft. So you get zero.

The same argument shows that you get zero for a wave of period 2. For a pulse length of 3 or 5 or anything that does not have a factor in common with 16, a similar argument shows that zero cannot happen.

One way to fix this, similar to what you are doing, is to introduce a small element into the pulse function that doesn't affect the convolution to any notable extent but does not allow a zero in the fft. The annotated line below does that. If the length of the pulse is an even number, including that line gives a good result and commenting it out gives NaNs.

%original signal

x = [1 2 3 4 5 6 5 4 3 2 1];

pulses = 4;

h = ones(1,pulses);

h(1) = h(1) + 1e-8; % assure nonzero fft

%creation of convolved signal

y = conv(x,h);

N = pow2(nextpow2(length(y)))

y(end+1:N) = 0 % zerofill

h(end+1:N) = 0; % zerofill

z = fft(y);

zh = fft(h);

x_decon = ifft(z./zh);

figure(1)

plot(x,'linewidth',2)

hold on

plot((x_decon),'d')

hold off

Keith Hernandez
on 2 Nov 2018

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.