% Example ID81.m
% Comprehensive example on identification
% Copyright Babu Joseph, 6/96 Washington University, St.Louis, 
%
% The purpose of this program is to illustrate the various steps involved
% in identifying a SISO process
% We will look at the following steps
%
% Step 1. Preidentication test and data analysis
% Step 2. Design of PRBS input signal
% Step 3. Prefiltering of the data
% Step 4. Model fitting and validation
% Step 5. Accuarcy of model
% Step 6. Accuarcy of model parameters
% Step 7. Comparison of model step response with actual plant response
%
% We will use data given in file called id81pt.dat for pretest data
% and id81.dat for actual identification test data
%
clg;clc; hold off;echo on;
%
% load prelim test data
%
load ch16ex5_prelim_test.mat;
time= t;u=pwr;y=temp;
ndata=length(u);
idplot([ y u]);
subplot(2,1,1);ylabel('Measurement, Volts');
subplot (2,1,2); ylabel('Power');xlabel('Time, seconds');%

pause % to see the plot of the input and ouput data to get an idea of the response
%

% Note that the change in u was made about midpoint in the data
% Since u was constant before this point an examination of y before the step change
% was made will give us an idea of the noise content in the signal
%
% First extract the data before change in u was made
u1=u(1:ndata/2);y1=y(1:ndata/2);umean=mean(u1);ymean=mean(y1);
%
t_samp=1; % define sample time
% first look at the frequency content of output signal to determine process and noise characterisitcs
%

n1=length(y); % no of points in the data
z1=fft( y-ymean) ;% compute fourier transform
% This gives a complex number the magnitude of which is the power spectrum
% The frequency ranges from 2pi/(t_samp*n1) to 2pi/(t_samp)
% z1 is symmetric about the nyquist frequency pi/t_samp
% Hence only the frequency below the nyquist frequency is relevant
%
freq1=0:2*pi/(t_samp*n1): pi/t_samp;  % the frequency axis
amp=abs(z1(1:n1/2+1))*2/n1; % power spectrum normalized so that unit amplitude
%                            in the freq domain corresponds to unit amplitude
%                            in the time domain
%
loglog(freq1,amp,'-');
xlabel('frequency, rad/sec');ylabel('amplitude');title('Fourier Transform of the noise in the input signal');
pause % to see the frequency response

%

% Now examine just the noise;

n1=length(y1); % no of points in the data
z1=fft( y1-ymean) ;% compute fourier transform
% This gives a complex number the magnitude of which is the power spectrum
% The frequency ranges from 2pi/(t_samp*n1) to 2pi/(t_samp)
% z1 is symmetric about the nyquist frequency pi/t_samp
% Hence only the frequency below the nyquist frequency is relevant
%
freq1=0:2*pi/(t_samp*n1): pi/t_samp;  % the frequency axis
amp=abs(z1(1:n1/2+1))*2/n1; % power spectrum normalized so that unit amplitude
%                            in the freq domain corresponds to unit amplitude
%                            in the time domain
%
loglog(freq1,amp,'-');
xlabel('frequency, rad/sec');ylabel('amplitude');
pause % to see the frequency response
%
% From this we can see that the noise is increasing with frequency
% with a peak at w= .3 rad/sec with a peak amplitude of about .08
%  We can apply various low pass filters
% to the data to see if we can recover the step response

% 
% let us filter the data and see the response
%
% assume sample time=1 
clg; hold off;
yf=zeros(ndata,1);

wb=0.5;wbs=num2str(wb); % Desired bandwidth of filter
yf=butterworth_filter(y-ymean,wb,t_samp); % This is an mfile that I created to do second order filtering
 % see plot of filtered data
%
plot(yf);title(['filtered data with cutoff-freq= ',wbs]);
xlabel('Time, sec');ylabel('Deviation in output, volts');

pause;
clg;
%
% The step occured at 200 sec and the new steady state is reached in about 50 seconds
% which shows a time constant of about 20 seconds
% From this we can guess that the signal to noise ratio is about 3.
% This information is used in the design of the input test signal for the 
% identification test.
% From results of example id.4.1 we have the following data
% wmin=.03, wmax=.20 
% We can use an input amplitude of 1  and sample time of 1 sec
%
wmin=.03;wmax=.20;amp=.2,t_samp=1;

utest=prbs(wmin,wmax,t_samp,amp);
% we will use two cycles of the prbs 
% the first cycle will be used for training ( parameter estimation)
% and the second for model testing
 ut=[utest' utest']';
% next we peturb the process with this signal to generate the actual identification data
%
load ch16ex5_training_run.mat;  % this loads the training data set
y=temp; u=pwr;ylen=length(y);
ymean=mean(y);y=y-ymean;umean=mean(u);u=u-umean;
idplot([ y u]);
subplot(2,1,1);ylabel('Measurement, Volts');
subplot (2,1,2); ylabel('Power');xlabel('Time, seconds');%


pause; % to see the plant response;
% Notice that there was an unknown disturbance of large magnitude entering
% the process around t=300 seconds. Hence we must remove this data before
% fitting a model
y=y+ymean;u=u+umean;% add back the mean to get original data
y=y(1:300);u=u(1:300);time=time(1:300);
ymean=mean(y);y=y-ymean;umean=mean(u);u=u-umean;uorig=u;yorig=y;

idplot([ y u]);
subplot(2,1,1);ylabel('Measurement, Volts');
subplot (2,1,2); ylabel('Power');xlabel('Time, seconds');%

pause; % to see the plant response after bad data was removed;
%
% next we will apply a prefilter to remove the high frequency noise
%
% From our input design we know that the noise is dominant at w > .5 rad/sec
% Hence apply a low pass filter with cutoff at w=.5
%
wb=.5;

yf=butterworth_filter(y,wb,t_samp);
uf=butterworth_filter(u,wb,t_samp);
idplot([yf uf]);
subplot(2,1,1);ylabel('Measurement, Volts');
subplot (2,1,2); ylabel('Power');xlabel('Time, seconds');%

pause;%  to see filtered data 
%
% Now fit a model using commands from identification toolbox
%
ndata=length(uf);
z=[yf uf];
th1=arx(z,[1 1 1]); % fit second order model with delay of 2
[num1,den1]=th2tf(th1);sysd1=tf(num1,den1,1)
[s1,t1]=step(sysd1); % get step response of first order system
th2=arx(z,[2 2 1]); % fit second order model with delay of 2
[num2,den2]=th2tf(th2);sysd2=tf(num2,den2,1)
[s2,t2]=step(sysd2); % get step response of second order system
th3=arx(z,[3 3 1]); % fit second order model with delay of 2
[num3,den3]=th2tf(th3);sysd3=tf(num3,den3,1)
[s3,t3]=step(sysd3); % get step response of third order system

th1=sett(th1,t_samp); % specify sampling interval
th2=sett(th2,t_samp); % specify sampling interval
th3=sett(th3,t_samp); % specify sampling interval


;pause

% to see how well model fits the validation data
%
load ch16ex5_validation.mat; % this loads the validation data

%
% The plant data  is too noisy
% let us filter out high frequency noise and compare
yt=temp;ut=pwr;
ytmean=mean(yt);yt=yt-ytmean;utmean=mean(ut);ut=ut-utmean;

wb=.5; ytf=butterworth_filter(yt,wb,t_samp);
utf=butterworth_filter(ut,wb,t_samp);
idplot([ytf utf]);
subplot(2,1,1);ylabel('Measurement, Volts');
subplot (2,1,2); ylabel('Power');xlabel('Time, seconds');%

pause; % to see filtered validation data
ymodel1=idsim(utf,th1); % compute model response
ymodel2=idsim(utf,th2); % compute model response
ymodel3=idsim(utf,th3); % compute model response

plot(ytf); hold on; plot(ymodel1); hold off; % compare model with plant
ylabel('Ouput, volts'); xlabel('Time, sec');
pause % to see comparison of model and plant for validation set
% let us also compute the least square error for the test data
%
lse1=norm(ytf-ymodel1)^2
lse2=norm(ytf-ymodel2)^2
lse3=norm(ytf-ymodel3)^2

pause;
% Normally we will pick a model with the smallest lse for the validation data
%
% Now let us try fitting using POLYID software provided with this book
%
time=time(1:300);
[hm, hh, hl, K, theta, tau]=polyid (yorig, uorig, time);
% compute model output
ymodelp=conv(utf,hm*t_samp);ymodelp=ymodelp(1:length(utf));
plot(ytf); hold on; plot(ymodelp); hold off; % compare model with plant
lse=norm(ytf-ymodelp)^2
ustep=ones(length(hm),1);
sp=conv(ustep,hm*t_samp);sp=sp(1:length(hm));

plot(sp);hold on;plot(t1,s1,'.');plot(t2,s2,'-');plot(t3,s3,'+');
ylabel('Step response ouput, volts'); xlabel('Time, sec');

[lse1 lse2 lse3 lse]
d2c(sysd1)
d2c(sysd2)
d2c(sysd3)

pause;
close all