Student project- problems using Matlab's filter command

J

Thread Starter

Javier Fontalvo

Dear community:

I'm new here, I hope I'm not incurring in any fault. For a project at the university I need to do system identification in Matlab using some given data. I recieve data in four columns: time, imput with removed offset, actual input and actual output. I need to remove the offset from the actual output, and using the half of both input and output with no offset data, find the values of the model parameters using the Linear Least Squares method. Up to here I think I made things right. To check how good is the estimation, I have to check the response of the calculated system using the filter command.

Here comes my problem: I only obtain as output an impulse (from zero to infinite) which does nothing to do with the input signal. I suspect the mistake is related with the filter command, but I cannot see where it is. I would really appreciate if you can help me with some advice to find the problem.

Code:<pre>
clc
clear
load Student_data.mat
%%
vtime=Log:),1);% takes time data from Log
per=(vtime(2)-vtime(1))/1000 %calculates the period
time=[0:per:per*length(vtime)-per]; %creates a new time vector starting from zero
not_off_input=Log:),2); %takes input with no offset data from Log
off_input=Log:),3);%takes input with offset data from Log
off_output=Log:),4);%takes output with offset data from Log
plot(time,off_input); hold on%
plot(time, off_output);hold on%
%%
avv=mean(off_output(1:168))%calculates the average of the output data till the 168th element
not_off_output=off_output-avv;% substracts the mean to obtain the output data with no offset
figure
plot(time,not_off_input);hold on
plot(time, not_off_output);hold on
%%
col1=not_off_output(5:(length(not_off_output)/2)-1); % first column of the leas squares parameters matrix, from the 5th to the half-matrix element
col2=not_off_output(4:(length(not_off_output)/2)-2);
col3=not_off_input(5:(length(not_off_output)/2)-1);
col4=not_off_input(4:(length(not_off_output)/2)-2);
psi=[col1 col2 col3 col4];
y=not_off_output(6:length(not_off_output)/2);%otputs matrix
theta=y\psi;% matrix of errors
a1=-theta(1);
a2=-theta(2);
b1=theta(3);
b2=theta(4);
%%
figure
out_fil=filter([0 b1 b2],[1 a1 a2],not_off_input);% use of filter command to simulate the data
plot(time,out_fil);hold on
clc
clear</pre>
Thanks very much in advance for your advice.

Best regards,
Javier
 
Dear Javier,

Can you explain what are you trying to do from the line "col1=not_off_output(5:..." onwards?

Gerrit.
 
J

Javier Fontalvo

> Can you explain what are you trying to do from the line
> "col1=not_off_output(5:..." onwards?

Hi Gerrit, thanks for your response.

I want to use the iteration of equations:
output(1)=-a1*output(0)-a2*output(-1)+b1*input(0)+b2*input(-1)
Then output(2)=-a1*output(1)-a2*output(0)+b1*input(1)+b2*input(0) and so on, in form of the matrix equation y=theta*psi, to find the parameters a1, a2, b1 and b2.

I took the sixth element of the output as output(1) to avoid estimate the values output(0) and output(-1). That is why I create the columns of theta starting from 5 (output(0)) and 4 (output (-1)).

As a lecturer's requirement, I have to use only the first half of de input and output data to use the least squares method. That is why the vector y goes till <i>length(output/input)/2</i> and the columns of theta go till <i>length(output/input)/2-1</i> and <i>length(output/input)/2-2</i>

Then, once I've built the vector y and the matrix theta using half of the data, I use the operator to calculate the vector psi, which contains the parameters -a1, -a2, b1 and b2. And then I change the sign of first two to obtain a1 and a2.

Finally, I use the filter command to simulate the whole data with the parameters obtained with the least squares method. In this point I'm not so sure how the filter command works, but looking at the command's help it is sort of passing the data through a discrete equation that uses these parameters. I hope it is better explained now, as my English is still not so good. Thanks again for your help,

Regards,

Javier
 
Hello Javier,

OK, your explanation helps a lot to interpret your code. First, there are two secondary issues (I think) for you to consider:

Issue 1. Why do you calculate the average of the output data up to index 168? This is in line:

avv=mean(off_output(1:168))%...

Why not index 'length(off_output)' or half that?

Issue 2:

> I took the sixth element of the output as output(1) to avoid estimate the
> values output(0) and output(-1). That is why I create the columns of theta
> starting from 5 (output(0)) and 4 (output (-1)).

I still do not understand why you start 'col2' and 'col4' at index 4 instead of index 1 (and 'y' at index 6 instead of 3, etc.). However, any change here should not significantly affect the outcome of your estimation.

<snip>

Issue 3:
> Then, once I've built the vector y and the matrix theta using half of the data,
> I use the operator to calculate the vector psi, which contains the
> parameters -a1, -a2, b1 and b2. And then I change the sign of first two to obtain
> a1 and a2.

<snip>

You are confusing theta and psi here, compared to the code. You also miss out the operator itself. I think that the first sentence as you intended it should be: "Then, once I've built the vector y and the matrix psi using half of the data, I use the operator \ to calculate the vector theta..."

The key of your problem seems to be here, in the calculation of theta: "theta=y\psi;%..." You need to have a look at the help for this operator (search for 'mldivide' in the MATLAB help) and check whether the above is correct or whether it should be "theta=psi\y;%...".

Your use of the filter function looks correct at first sight, but we can look into that if necessary once you have checked the theta calculation. Please let us know how you get on.

Regards,
Gerrit.
 
J

Javier Fontalvo

Dear Gerrit,

As you wrote in your answer, the mistake was in the use of the \ operator, I was putting the matrixes inverted. Once I changed the order, I got the correct simulation :)

I really appreciate your help with this problem, thanks for your time. I hope to participate frequently in this excellent forum.

Best Regards,

Javier
 
Top