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']);