How can I determinate the delay and the time constant from script?

1 view (last 30 days)
Hi!
I have an input and an output vector, for example:
in=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1]
out=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 2; 12; 17; 20; 21; 22; 22; 22; 22; 22; 22; 22; 22; 22; 22; 22; 22]
If you want to see it in a figure, type the following code to the Command Window:
in=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1];
out=[0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 2; 12; 17; 20; 21; 22; 22; 22; 22; 22; 22; 22; 22; 22; 22; 22; 22];
x=(1:1:30);
plot(x,in,'red')
hold
plot(x,out,'blue')
hold
I would like to determinate the delay and the time constant. I have System Identification Toolbox, and I can determinate the delay in GUI mode (Td=2.864), but I would like to determinate it in a matlab script. Which functions have to use, and how? Maybe I have to use the delayest function, but it doesn't works in this way:
DATA=iddata(out,in,0.1)
delayest(DATA)
The result is 0.
So I think this is not the right way. Can somebody write me an example to this two vector?
Thanks in advance!

Answers (4)

Image Analyst
Image Analyst on 18 Feb 2013
What about just thresholding and using find()? Let's say that you say a signal starts when its amplitude exceeds some threshold:
thresholdValue = 0.5;
inStartTime = find(in > thresholdValue, 1, 'first')
outStartTime = find(out > thresholdValue, 1, 'first')
lagTime = outStartTime - inStartTime

Ákos
Ákos on 19 Feb 2013
It is not too simple, because the problem, that there is a big noise on the signal, the initial value is oscillate with random amplitudes between 0 and even 0.3 * average of final value, and the final value is oscillate too.
But I found the way, how can I use the functions of System Identification Toolbox right:
DATA=iddata(out,in,0.1);
model=arx(DATA,[1,1,1]);
type = idproc('P1D');
result=procest(DATA,type);
fprintf('Td= %f sec, Tp1= %f sec', result.Td, result.Tp1);
But really thank you for trying to help.

Rajiv Singh
Rajiv Singh on 19 Feb 2013
Try impulse response estimation based approaches, assuming you are able to find a good FIR model for data (which you can ascertain using COMPARE)
m=impulseest(DATA); h = impulseplot(m); showConfidence(h,3)
The first response value that is significantly out of 3 sd region is at t = 0.3 which indicates a 3 sample (0.3 sec) delay. See:
  1 Comment
Image Analyst
Image Analyst on 20 Feb 2013
This is not an answer. It is supposed to be a comment to someone but we don't know who because you posted it as an independent answer. Please copy it to a comment where it belongs.

Sign in to comment.


Ákos
Ákos on 20 Feb 2013
Edited: Ákos on 20 Feb 2013
Thank you, it's a helpful answer, but according to the compare funcion, the DATA, and the impulseest(DATA) does not fit well.
Here is my measured datas:
If you Import the data file as a matrix, and you use the following code, you can see my system (Input in the 5th column (this is a duty cycle [%]), output in the 3rd column):
in=data(:,5);
out=data(:,3);
x=0.1:0.1:length(in)/10;
plot(x,out,'blue');
hold
plot(x,100*in,'red');
hold
Blue line shows the output, and the red line shows the input(*100). This is a nonlinear system, and i would like to determinate the delay and the time constant.
You can choose one part of the full system with this code (you can change the value of variable 'which'):
which=5;
dc_begin_vector=[1 638 1925 3237 4115 5071 6043 7065 7985 9326 10512 11802 14400 15547 16920 18249 19447 20715 21987 23319 24667 25925 27095];
sample_out_full=data(:,3);
sample_in_full=data(:,5);
sample_out=sample_out_full((dc_begin_vector(which)-100):((dc_begin_vector(which+1))-1));
sample_in=sample_in_full((dc_begin_vector(which)-100):((dc_begin_vector(which+1))-1));
sample_out=sample_out-min(sample_out);
sample_in=sample_in-min(sample_in);
x_axis=0.1:0.1:(length(sample_out)/10);
plot(x_axis,sample_out,'blue')
hold
plot(x_axis,100*sample_in,'red')
hold
If you can determinate the delay and the time constant to this one part, it's good for me.
  1 Comment
Image Analyst
Image Analyst on 20 Feb 2013
This is not an answer. It is supposed to be a comment to someone but we don't know who because you posted it as an independent answer. Please copy it to a comment where it belongs.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!