Reconstructed output mismatch for LTI system

Question

I have a system with measured input (u) and output (y).

I assume that this is an linear time-invariant (LTI) system and I want to find its impulse response (ir).

This is the Matlab code I use to find ir:

ir = real(ifft(fft(y)./fft(u)));

And this is how I reconstruct the system output:

Y = conv(ir, u);

However, as can be seen below, the reconstructed output differs greatly from the true output. What am I missing here?

The full matlab code including data can be found here

Answer 1

The problem is that the formula $$H(f) = \frac{Y(f)}{X(f)},$$ where $H(f)$ is the frequency response, and $X(f)$, $Y(f)$ are the input and output, is valid only for frquencies $f$ where $|X(f)| \neq 0$.

Furthermore, if $|X(f_0)| \approx 0$, you're going to run into numerical errors, since $|H(f_0)| \rightarrow \infty$.

Your input signal has a bunch of places where it's zero, or close to zero. I wrote a little function that takes care of those cases:

function h = impres(x,y)

X = fft(x);
Y = fft(y);

for i = 1:length(X)
    if abs(X(i)) > 10
        H(i) = Y(i)/X(i);
    else
        H(i) = 0;
    end
end

h = ifft(H);

end

While not perfect, and could be tweaked, I get:

where the red line is the estimated output.

One more thing to be aware of: the equivalent in time of the IDFT is the circular convolution, not the linear convolution. In this case there is not too much difference, but you want to be as correct as you can in your code.