prog='TScorr';
%----------------------------------------------------------------------------------------
%                TIME SERIES CORRELATION - relation between two time signals
%========================================================================================
%
% DESCRIPTION:      This program has been developed to view and analyse a time series of
% ------------      a measured signal. The main features are:
%
%                   - viewing the time steries as a graph.
%                   - presentation of data points as histogram (discrete probability density).
%                   - calculation of power spectral density.
%                   - calculation and displaying the mean values and standard deviations of sub-series.
%                   - Identification of maximum and minimum in sub-series.
%
%                   Furthermore, some signal processing can be performed on the original data,
%                   namely detrending and high or low-pass filtering. The difference between the
%                   initial and the altered signal can be visualised in the graph of the time
%                   series and the spectral density the spectral density. 
%
% Program ID:
% -----------
% File name       : THA4.m
% Author          : Holger Koss (hko)
% Development Log : 2009-04-28  hko   Basic structure of the program
%                   2012-05-04  hko   Adoption to a general tool to first analysis and
%                                     treatment of a time history
%                   2012-10-29  hko   Adoption for course 11374
%
% Necessary files  : "TimeHistory"  - Ascii file containing in the first column the
%                                     time axis and in the second column the time history
%                                     of the investigated quantity. Both columns are in
%                                     model scale.
%
%                   "cpcent00.dat"  - File with 18 time series of pressure coefficients 
%                                     measured in a wind tunnel test on a model low-rise
%                                     building.
%========================================================================================
 close all
 clear all
%========================================================================================
% 1) READING OF INPUT DATA:
%    ----------------------
%    This version of the program is prepared to read a file with a matrix 
%    of 18 time series of measured pressure coefficients. In this case the
%    time series needs to be generated to plot the series.
%

% 1.1) Reading Data from Pressure Measurements on Low-rise Building:
%      -------------------------------------------------------------
%      Here, we read a file with 18 time series of pressure coefficients
%      measured on a wind tunnel model of a low-rise building. Each of these 
%      time series is a signal of a stochastic process. The input file does
%      not contain the time axis. To plot the process the time axis needs
%      to be generated separately.
%

    
FileName = 'cpcent00.dat';             % Name of input file
NFFTcase = 5*1024;                     % Filter depth of the fft-routine N*1024
SignalA = 1;                           % Number of 1st signal to be compared (1-18)
SignalB = 18;                           % Number of 2nd signal to be compared (1-18)
fid      = fopen(FileName,'r');
Fsamp    = 1600;                       % Sample frequency in 
qhmwk    = fscanf(fid,'%e',[1 1]);     % velocity pressure [kN/m^2]
cp       = fscanf(fid,'%e',[18 inf]);  % Matrix with pressure coefficient time series
Series   = cp';                        % Transposed matrix where each column is a signal
[m1 n1]  = size(Series);               % Number of rows (m1) and columns (n1)
status   = fclose(fid);
    
DT       = 1/Fsamp;                    % Calculation of time step
for i=1:m1                             % Generation of time axis with m1 steps
    TAx(i) = (i-1)*DT;
end
    
XA = Series(:,SignalA);                % Saving selected data to input vector
XB = Series(:,SignalB);                % Saving selected data to input vector


% NOTE: At this point in the program you should have following information available:
%
%       Fsamp    = sample frequency in [Hz]
%       DT       = time step between data points [s] = 1/Fsamp
%       SignalNo = Number of signal that has been chosen to be analysed - saved as X0(i)
%       m1       = number of data points (time steps)
%       NFFTcase = FFT filter length, determine by trial, shall not exceed m1
%       TAx(i)   = Vector with values for time axis (length = m1)
%       X0(i)    = Vector with data of stochastic process (length = m1)
%
%       Disregarding what data you want to analyse, just make sure that after reading the
%       input file you define the parameters and vectors listed above.

%========================================================================================
% 2) ANALYSIS SETTING:
%    -----------------
%    The stochastic process is now saved as X2(t) on which all subsequent
%    analysis will be peformed.
%
%    2.1 Statistical Parameter of 1st Time History (XA):
%    ---------------------------------------------------

XAmean    = mean(XA);              % mean value of the parent time history
XAstd     = std(XA);               % standard deviation of the parent time history
XAvar     = XAstd^2;               % corresponding variance
XAmax     = max(XA);               % maximum peak value occuring in parent time history
XAmin     = min(XA);               % corresponding minimum peak value
Tend      = TAx(m1);               % duration of (both!) time histories

%    2.2 Statistical Parameter of 2nd Time History (XB):
%    ---------------------------------------------------

XBmean    = mean(XB);              % mean value of the parent time history
XBstd     = std(XB);               % standard deviation of the parent time history
XBvar     = XBstd^2;               % corresponding variance
XBmax     = max(XB);               % maximum peak value occuring in parent time history
XBmin     = min(XB);               % corresponding minimum peak value

%    2.3 Statistical Parameter of 2nd Time History (XB):
%    ---------------------------------------------------

ABcorr = corr(XA,XB);

% R = corrcoef(X) returns a matrix R of correlation coefficients calculated
% from an input matrix X whose rows are observations and whose columns are variables
Mcorr  = corr(Series);   


disp(['Data read from file:',FileName])
fprintf(1,' \n');
fprintf(1,' Correlation cofficient: %g\n',ABcorr);


%========================================================================================
% 3) GRAPHICAL DISPLAY OF THE EXTREME VALUE ANALYSIS:
%----------------------------------------------------
% Display Definitions
% -------------------

scrsz = get(0,'ScreenSize');

%====================================================================================
figure('Name','Time History A','Position',[5 0.50*scrsz(4) 0.7*scrsz(3) 0.45*scrsz(4)])

% Setting axis range

Y1max = XAmin+(XAmax-XAmin)*1.2;
Y1min = XAmin-(XAmax-XAmin)*0.02;
DY1   = Y1max-XAmax;

Ex   = floor(log10(Tend));
X1min = 0.;
X1max = ceil(Tend/(10^Ex))*10^Ex;

plot(TAx,XA,'-b');
hold on

plot3([TAx(1),TAx(m1)],[XAmean,XAmean],[1,1],'--m');
hold on
plot3([TAx(1),TAx(m1)],[XAmean+XAstd,XAmean+XAstd],[1,1],'--g');
hold on
plot3([TAx(1),TAx(m1)],[XAmean-XAstd,XAmean-XAstd],[1,1],'--g');
hold on

text(Tend/20,Y1max-0.3*DY1,['Time History File is "',FileName,'"',' ; Signal No.:',num2str(SignalA)],'FontSize',9) 
text(Tend/20,Y1max-0.7*DY1,['Time History duration is ',num2str(Tend),' sec'],'FontSize',9) 

text(0.98*Tend,Y1max-0.3*DY1,['Mean value of XA(t): ',num2str(XAmean)],'FontSize',9,'HorizontalAlignment','right')
text(0.98*Tend,Y1max-0.7*DY1,['Standard deviation of XA(t): ',num2str(XAstd)],'FontSize',9,'HorizontalAlignment','right')

xlabel('time [s]');
ylabel('ordinate of variable XA(t)');
title('1st TIME HISTORY in COMPARISON');

axis([X1min X1max Y1min Y1max]);

eval(['print -dtiff -zbuffer TimeHistA']);

%====================================================================================
figure('Name','Time History B','Position',[5 35 0.7*scrsz(3) 0.45*scrsz(4)])

Y2max = XBmin+(XBmax-XBmin)*1.2;
Y2min = XBmin-(XBmax-XBmin)*0.02;
DY2   = Y2max-XBmax;

Ex   = floor(log10(Tend));
X2min = 0.;
X2max = ceil(Tend/(10^Ex))*10^Ex;

plot(TAx,XB,'-b');
hold on

plot3([TAx(1),TAx(m1)],[XBmean,XBmean],[1,1],'--m');
hold on
plot3([TAx(1),TAx(m1)],[XBmean+XBstd,XBmean+XBstd],[1,1],'--g');
hold on
plot3([TAx(1),TAx(m1)],[XBmean-XBstd,XBmean-XBstd],[1,1],'--g');
hold on

text(Tend/20,Y2max-0.3*DY2,['Time History File is "',FileName,'"',' ; Signal No.:',num2str(SignalB)],'FontSize',9) 
text(Tend/20,Y2max-0.7*DY2,['Time History duration is ',num2str(Tend),' sec'],'FontSize',9) 

text(0.98*Tend,Y2max-0.3*DY2,['Mean value of XB(t): ',num2str(XBmean)],'FontSize',9,'HorizontalAlignment','right')
text(0.98*Tend,Y2max-0.7*DY2,['Standard deviation of XB(t): ',num2str(XBstd)],'FontSize',9,'HorizontalAlignment','right')

xlabel('time [s]');
ylabel('ordinate of variable XB(t)');
title('2nd TIME HISTORY in COMPARISON');

axis([X1min X1max Y2min Y2max]);

eval(['print -dtiff -zbuffer TimeHistB']);

%====================================================================================
figure('Name','Correlation','Position',[0.715*scrsz(3) 0.50*scrsz(4) 0.29*scrsz(3) 0.45*scrsz(4)])


plot(XA,XB,'.b','MarkerSize',3)
hold on

%
% Calculation of coordinate values for graph design

X1 = min(XA);
X2 = max(XA);
Y1 = min(XB);
Y2 = max(XB);
DX = abs((X2-X1)/20);
DY = abs(Y2-Y1)/20;

%
% Calculating the line coordinates for full correlation (corr=1) with
% refernce in intersection point of both mean values.
%

plot ([X1 X2],[XBmean-(XAmean-X1) XBmean+(X2-XAmean)],'-k');
hold on
%plot ([X1 X2],[XBmean+(XAmean-X1) XBmean-(X2-XAmean)],'-.k');
plot ([XAmean+(XBmean-Y1) XAmean-(Y2-XBmean)],[Y1 Y2],'-.k');
hold on

title('CORRELATION');
xlabel('signal A');
ylabel('signal B');

text(X1+DX,max(XB)-DY,['Corr = ',num2str(ABcorr)],'FontSize',9) 

%
% Plotting the MEAN value and STANDARD DEVIATION of 1st signal XA
%
plot3([XAmean,XAmean],[Y1,Y2],[1,1],'--m');
hold on
plot3([XAmean+XAstd,XAmean+XAstd],[Y1,Y2],[1,1],'--g');
hold on
plot3([XAmean-XAstd,XAmean-XAstd],[Y1,Y2],[1,1],'--g');
hold on

%
% Plotting the MEAN value and STANDARD DEVIATION of 2nd signal XB
%
plot3([X1,X2],[XBmean,XBmean],[1,1],'--m');
hold on
plot3([X1,X2],[XBmean+XBstd,XBmean+XBstd],[1,1],'--g');
hold on
plot3([X1,X2],[XBmean-XBstd,XBmean-XBstd],[1,1],'--g');
hold on


axis([X1 X2 Y1 Y2]);
axis equal;

eval(['print -dtiff -zbuffer Correlation']);




