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

3 次查看(过去 30 天)
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!

回答(4 个)

Image Analyst
Image Analyst 2013-2-18
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 2013-2-19
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 2013-2-19
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 个评论
Image Analyst
Image Analyst 2013-2-20
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.

请先登录,再进行评论。


Ákos
Ákos 2013-2-20
编辑:Ákos 2013-2-20
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 个评论
Image Analyst
Image Analyst 2013-2-20
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.

请先登录,再进行评论。

Community Treasure Hunt

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

Start Hunting!

Translated by