% get_response.m

%

% This MATLAB function uses Fourier analysis to estimate

% the response function of a system from the provided

% input and output data.

%

% K. Beers. MIT ChE. 8/14/2002

function [r,r_FT,iflag] = ...

    get_response(dt,t0,x,y,iplot);

iflag = 0;

% We first obtain the discrete Fourier transforms

% of the input and output signals

x_FT = fft(x);

y_FT = fft(y);

% We now compute the Fourier transform of the

% response function.  If x_FT is non-zero for

% a particular frequency, we set r_FT = 0 for

% this frequency.

N = length(x_FT);

r_FT = complex(zeros(1,N),zeros(1,N));

cutoff = sqrt(eps);

for m=1:N

    if(abs(x_FT(m) < cutoff))

        r_FT(m) = 0;

    else

        r_FT(m) = y_FT(m)/x_FT(m);

    end

end

% apply inverse Fourier transform to get

% time-space response function

r = real(ifft(r_FT));

% if desired, plot results

if(exist('iplot'))

    if(iplot)

        figure;

        t_val = linspace(0,(N-1)*dt,N);

        t_plot_min = 0;

        t_plot_max = t_val(N);

        % input signal

        subplot(2,2,1);

        plot(t_val,x);

        buffer = 0.05*(max(x)-min(x));

        axis([t_plot_min t_plot_max ...

                (min(x)-buffer) (max(x)+buffer)]);

        xlabel('t');

        ylabel('x(r) - input signal');

        title('estimated response function');

        % output function

        subplot(2,2,2);

        plot(t_val,y);

        xlabel('t');

        ylabel('y(t) - output signal');

        buffer = 0.05*(max(y)-min(y));

        axis([t_plot_min t_plot_max ...

                (min(y)-buffer) (max(y)+buffer)]);

        % reponse function (real-time)

        subplot(2,2,3);

        plot(t_val,r);

        xlabel('t');

        ylabel('r(t) - estimated response');

        buffer = 0.05*(max(r) - min(r));

        axis([t_plot_min t_plot_max ...

                (min(r) - buffer) (max(r) + buffer)]);

    end

end

iflag = 1;

return;