% get_response_example.m

%

% This MATLAB script demonstrates the use of

% Fourier analysis to estimate the response

% function of a damped-spring system.

%

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

% First, we generate the exact response function

% that we will use to synthesize the output signal.

% We will also want to compare this exact response

% function to that estimated from only the input

% and output data.  We use the script file

% damped_spring_response.m that was written

% previously.

damped_spring_response;

i_r_zero = 1;  % position of t = 0 in response signal

% Now, we generate as an input signal a square pulse

% of duration 1.

x = zeros(size(t_val));

j_step = max(find(t_val < 1));

x(1:j_step) = 1;

% We use convolution to synthesize the output data from

% the known response function.

iplot = 1;

y = compute_convolution(dt,t0,x,r,i_r_zero,iplot);

% we then bring the input vector up to the same length

% by padding with zeros

num_pad = length(y) - length(x);

x = [x zeros(1,num_pad)];

% From this point on, we assume that we only have knowledge

% of the measured input x(t) and the measured output y(t),

% and that we do not know yet the response function.

% We now call the routine get_response() to estimate the

% response function from the input and output data.

iplot = 1;

[r_est,r_est_FT] = get_response(dt,t0,x,y,iplot);

% We now graph the known and estimate response functions

% in Fourier and time space.

figure;

% actual response function

subplot(2,1,1);

N_plot = length(t_val);

plot(t_val,r);

xlabel('t');

ylabel('r(t) - actual');

title('Actual and estimated response functions');

% estiamted response function

subplot(2,1,2);

plot(t_val,r_est(1:N_plot));

xlabel('t');

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

% We now consider the effect of noise on the input and output

% signals on our estimation of the response function.

% We add to the input and output a Gaussian-distributed,

% uncorrelated noise.

noise_mag = 0.1;

x_noisy = x + noise_mag*randn(size(x));

y_noisy = y + noise_mag*randn(size(y));

[r_noisy,r_noisy_FT] = get_response(dt,t0, ...

    x_noisy,y_noisy,iplot);

% We now graph the known and estimate response functions

% in Fourier and time space.

figure;

% actual response function

subplot(2,1,1);

plot(t_val,r);

xlabel('t');

ylabel('r(t) - actual');

title('Actual and estimated response functions (with noise)');

% estimated response function

subplot(2,1,2);

plot(t_val,r_noisy(1:N_plot));

xlabel('t');

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

% We can obtain better estimates of the response function

% if we first use a filter to remove the noise from the

% input and output signals.

omega_crit = pi/dt;

omega_filter_max = 0.975*omega_crit;

x_filtered = filter_signal(dt,x,omega_filter_max);

y_filtered = filter_signal(dt,y,omega_filter_max);

[r_filtered,r_filtered_FT] = get_response(dt,t0, ...

    x_filtered,y_filtered,iplot);

% plot results when using filter on input and output signals

figure;

% actual response function

subplot(2,1,1);

plot(t_val,r);

xlabel('t');

ylabel('r(t) - actual');

title('Actual and estimated response functions (with filter)');

% estimated response function

subplot(2,1,2);

plot(t_val,r_filtered(1:N_plot));

xlabel('t');

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