Impulse response and signal reconstruction

Hello, I've been given two data sets to help characterise a single-input-single output physical system for some sensors.

As attached, I've filtered and plotted the input and output data for each data set. I'd like to be able to generate (I presume) an impulse response function that will allow me to back out the input signal from a given transient output signal.

I've tried a few ways to do this. Here is one example below. The result from estimating the input data from the generated impulse responses is clearly very poor. You can see this for both data sets A and B (i.e. short and long pulses).

Does anyone know if this is the right method? I'd have thought generating a transfer function and then impulse response from discrete data would've been straightforward but I'm clearly missing something.

Cheers!

% LOAD
load('data_A.mat'); 
time = data.time;
input_A = data.input;
output_A = data.output;

% TEST PARAMETERS
samplefreq = 5e6; 
npoints_pretrigger = 1000; 

% clean data
d_A = output_A - mean(output_A(1:npoints_pretrigger)); 


nx=1:10;  
sys_id  = n4sid(iddata(input_A,output_A,1/samplefreq),nx);  % nx ends up being selected as 1 or 2
[A, B, C, D] = ssdata(sys_id);
[num, den] = ss2tf(A, B, C, D);
h3 = impulse(tf(num, den), time_A); % Generate impulse response using transfer function

output_estimation = fftfilt(h3, output_A)   % recover input data

% another check to compare the system estimation to the input signal 
compare(iddata(input_A,X_A,1/samplefreq), sys_id)

Answer 1

[Maneet Kaur Bagga]

Hi,

As per my understanding you are aiming to characterize a system using the system identification to obtain an impulse response function that can be used to reconstruct the input signal from the output signal. Please refer to the below suggestions as a possible workaround.

Increase the range of "nx" to test different model orders and potentially use model selection criteria to choose the best order.

% Instead of nx=1:10, try a range and use a criterion for selection
nx_range = 1:10;
sys_id_models = cell(length(nx_range), 1);
fit_percent = zeros(length(nx_range), 1);

for i = 1:length(nx_range)
    sys_id_models{i} = n4sid(iddata(input_A, output_A, 1/samplefreq), nx_range(i));
    fit_percent(i) = sys_id_models{i}.Report.Fit.FitPercent;  % Get the fit percentage
end

% Select the model with the best fit
[~, best_idx] = max(fit_percent);
sys_id = sys_id_models{best_idx};

If overfitting is suspected, you may use "regularization" to prevent it.

% Use regularization (requires the System Identification Toolbox)
opt = n4sidOptions('Regularization', 'auto');
sys_id = n4sid(iddata(input_A, output_A, 1/samplefreq), 'best', opt);

Perform frequency domain analysis to check if the system behaves differently at various frequencies.

% Convert to frequency response data and plot Bode plot
frd_sys = idfrd(sys_id, linspace(0, samplefreq/2, 4096), 'FrequencyUnit', 'Hz');
bode(frd_sys);

If the impulse response is estimated well, use deconvolution directly to recover the input.

% Deconvolution to estimate the input signal
input_estimation = deconv(output_A, h3);

Introduce a non-parametric impulse response estimation as an alternative to the parametric method.

% Non-parametric impulse response estimation
h = impulseest(iddata(output_A, input_A, 1/samplefreq));
impulse(h, time);

Please refer to the MATLAB Documentation for further understanding:

impulseest: https://in.mathworks.com/help/signal/gs/impulse-step-and-ramp-functions.html

Hope this helps!