prog='THA5'; %---------------------------------------------------------------------------------------- % TIME HISTORY ANALYSIS - Detrend, Filter and Display %======================================================================================== % % 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 : THA5.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 % 2012-11-07 hko Reading additional types of files % 2013-11-11 hko Different formats of PSD % % 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. % % "18Signals.dat" - File with 18 time series of pressure coefficients % measured in a wind tunnel test on a model low-rise % building (based on cpcent.00). % % "cpcent.00" - Fragmented file with cp time series from 12 storm % events, each lead by the velocity pressure [kPa] % in the wind tunnel test at eaves height. % (Exercise: "Characteristic Wind Load Statistic on Low-rise % Buildings") % % "BendTS.txt" - Time series of the base bending moment [N] of a % high-rise building. % (Exercise: "Dynamic Response of a High-rise Building % to Wind Loading") %======================================================================================== close all clear all %======================================================================================== % 1) READING OF INPUT DATA: % ---------------------- % This version of the program is prepared to read two different input files. % A simple file with two columns, time axis (t) and variable X(t), and 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 = Read simple input file % 2 = Read matrix of cp-series % 3 = Read observation 10min-mean data Reykjavik harbour 2002 (hole year) % 4 = Reading on hour of anemometer data from Oeresund Bridge (30Hz) % 5 = Reading fragmented CpCent data files with wind tunnel speeds inbetween blocks % 6 = custom sequence Iread = 3; % 1.1) Reading simple data set: % ---------------------------------------------------------- % File with two columns: time axis and data of variable X(t). The % values in this file have been recorded at top of the pylon of the % Great Belt Bridge (TimeHistory). if Iread == 1 % FileName = 'Data.txt'; % Name of input file (10min wind Great Belt Bridge) % FileName = 'TimeHistory.dat'; % Name of input file (10min wind Great Belt Bridge) FileName = 'BendTS.txt'; % Name of input file (Base bending moment CDJ) % FileName = 'test1.txt'; % Name of input file (Base bending moment CDJ) % FileName = 'Tow2X5.273'; % Name of input file (Base bending moment CDJ) fid = fopen(FileName,'r'); Series = fscanf(fid,'%e',[2 inf]); % Matrix with time series Series = Series'; % Transposed matrix [m1 n1] = size(Series); % Number of rows (m1) and columns (n1) status = fclose(fid); disp(['Data read from file:',FileName]) TAx = Series(:,1); % Time series X0 = Series(:,2); % Saving selected data to variable vector SignalNo = 1; % Default value since there is only one signal to analyse NFFTcase = 5*1024; % Filter depth of the fft-routine N*1024 % NFFTcase = 512; % Filter depth of the fft-routine N*1024 Nwindow = 8; % Number of windows (default = 8) for pwelch SFD calculation DT = (TAx(10)-TAx(1))/9; % Calculation of time step Fsamp = 1/DT; % Calculation of sample frequency Nbin = 30; % Number of bins to generate histogram Nsub = 30; % Number of sub-series end % 1.2) 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. % if Iread == 2 %FileName = '18Signals.dat'; % Name of input file FileName = 'cpR6400.00'; % Name of input file fid = fopen(FileName,'r'); % Echo print on screen 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); disp(['Data read from file:',FileName]) SignalNo = 12; % Number of signal to be analysed (1-18) X0 = Series(:,SignalNo); % Saving selected data to input vector NFFTcase = 2*1024; % Filter depth of the fft-routine N*1024 Nwindow = 8; % Number of windows (default = 8) for pwelch SFD calculation Fsamp = 6400; % Sample frequency in DT = 1/Fsamp; % Calculation of time step Nbin = 100; % Number of bins to generate histogram Nsub = 11; % Number of sub-series for i=1:m1 % Generation of time axis with m1 steps TAx(i) = (i-1)*DT; end end % 1.3) Reading Data from Wind Records at Reykjavik Harbour, entire year 2002: % ---------------------------------------------------------------------- % The data are provided by a private weather station located in the harbour of Reykjavik. % The file contains amogst other 10 minutes mean and gust wind speeds continuously % recorded throughout the year 2002. The file format is given in the list below. % % Column Content % ----------------------------------------------------- % 1 Hours % 2 Minutes % 3 Day % 4 Month % 5 Year % 6 10 Minutes mean wind speed [m/s] % 7 Gust wind speed [m/s] % 8 Wind Direction [deg] % 9 Standard deviation of wind direction [deg] % 10 Atmospheric pressure [mbar] % 11 Air temperature [degC] % if Iread == 3 FileName = 'Reykjavik2002.txt'; % Name of the 1st time history file fid1 = fopen(FileName,'r'); Series = fscanf(fid1,'%g',[11 inf]); Series = Series'; [m1 n1] = size(Series); status = fclose(fid1); SignalNo = 6; % Number of signal to be analysed (1-18) X0 = Series(:,SignalNo); % Saving selected data to input vector NFFTcase = 10*1024; % Filter depth of the fft-routine N*1024 Nwindow = 8; % Number of windows (default = 8) for pwelch SFD calculation Fsamp = 0.001667; % Sample frequency in [Hz] DT = 600; % 10min mean data Nbin = 50; % Number of bins to generate histogram Nsub = 12; % Number of sub-series % Calculating a continous time index: % ----------------------------------- % (time index is generated in 10-minutes steps as the smallest time unit available) DayMon = [31 28 31 30 31 30 31 31 30 31 30 31]; % days per month for time axis for i=1:m1 if Series(i,4)==1; Month=0 ; end if Series(i,4)==2; Month=44640 ; end if Series(i,4)==3; Month=84960 ; end if Series(i,4)==4; Month=129600; end if Series(i,4)==5; Month=172800; end if Series(i,4)==6; Month=217440; end if Series(i,4)==7; Month=260640; end if Series(i,4)==8; Month=305280; end if Series(i,4)==9; Month=349920; end if Series(i,4)==10; Month=393120; end if Series(i,4)==11; Month=437760; end if Series(i,4)==12; Month=480960; end TAx(i) = Series(i,1)*60+Series(i,2)+(Series(i,3)-1)*24*60+Month; % time axis in minutes end end % 1.4) Reading record file from Oeresund Bridge Monitoring: % ---------------------------------------------------------------------- % The data are recorded on the Oresund Bridge and cover a period of one hour. % The input file has one sub array "wind" with 8 columns of following structure: % % Anemometer A: located at second shortest cable % 1st U component % 2nd V component % 3rd W component % 4th Azimuth % % Anemometer B: located at midspan % 5st U component % 6nd V component % 7rd W component % 8th Azimuth % % With V direction coinciding with geographical north, U coincide with Eastand W % is vertical direction Azimuth is the angle that define the position of the wind % flow along the hoorizontal plane. It is calculated from 0 to 360 clockwise along % the horizontal plane. 0 is the Geographical north. % Each columns is 1 hour long (108000) sampled at 30Hz. if Iread == 4 % FileName = 'wind_RawData_20120320_162217.mat'; % Starting 4.22pm % FileName = 'wind_RawData_20120320_172250.mat'; % Starting 5.22pm % FileName = 'wind_RawData_20120320_182323.mat'; % Starting 6.22pm FileName = '3 records of consecutive hours'; % Starting 6.22pm load 'wind_RawData_20120320_162217.mat'; Series1 = wind; [m1 n] = size(Series1); load 'wind_RawData_20120320_172250.mat'; Series2 = wind; [m1 n] = size(Series2); load 'wind_RawData_20120320_182323.mat'; Series3 = wind; [m1 n] = size(Series3); m3 = 3*m1; % All three files are of the same length! X0a1 = zeros(m1,1); X0a2 = zeros(m1,1); X0a3 = zeros(m1,1); X0b1 = zeros(m1,1); X0b2 = zeros(m1,1); X0b3 = zeros(m1,1); X0a = zeros(m3,1); % vector for combined series X0b = zeros(m3,1); % vector for combined series % Horizontal resulting conponent for anemometer A: for i=1:m1 X0a1(i) = sqrt(Series1(i,1)^2+Series1(i,2)^2); X0a2(i) = sqrt(Series2(i,1)^2+Series2(i,2)^2); X0a3(i) = sqrt(Series3(i,1)^2+Series3(i,2)^2); X0a(i) = X0a1(i); X0a(i+m1) = X0a2(i); X0a(i+2*m1) = X0a3(i); end % Horizontal resulting conponent for anemometer B: for i=1:m1 X0b1(i)= sqrt(Series1(i,5)^2+Series1(i,6)^2); X0b2(i)= sqrt(Series2(i,5)^2+Series2(i,6)^2); X0b3(i)= sqrt(Series3(i,5)^2+Series3(i,6)^2); X0b(i) = X0b1(i); X0b(i+m1) = X0b2(i); X0b(i+2*m1) = X0b3(i); end % X0 = X0a; % Three hours wind speeds at second shortest cable X0 = X0b; % Three hours wind speed at midspan % X0 = X0a1; % First hour wind speed at short cable % X0 = X0a2; % Second hour wind speed at short cable % X0 = X0a3; % Third hour wind speed at short cable % X0 = X0b1; % First hour wind speed at midspan % X0 = X0b2; % Second hour wind speed at midspan % X0 = X0b3; % Third hour wind speed at midspan % X0 = Series1(:,4); % One hour wind direction m1 = m3; % Record length: m1 for one hour % m3 for three hours SignalNo = 9999; % Combination of different signals NFFTcase = 20*1024; % Filter depth of the fft-routine N*1024 Nwindow = 8; % Number of windows (default = 8) for pwelch SFD calculation Fsamp = 30; % Sample frequency in [Hz] DT = 1/Fsamp; % Time step [s] Nbin = 100; % Number of bins to generate histogram Nsub = 3; % Number of sub-series for i=1:m1 % Generation of time axis with m1 steps TAx(i) = (i-1)*DT; end end % 1.5) Reading orginal formated data files for CpCent time series: % ---------------------------------------------------------------------- % The input file consists of time series from 18 signals. The time % series have a record length 4096 steps sampled at 1600Hz in a wind % tunnel test. The mean wind speed at which the test has been % performed is given as mean velocity pressure [kPa] at the start of the % block with 4096x18 data entries. All-in-all 12 of the above % described data sets, i.e. mean velocity pressure and signal time % series, are contained in the input file. The data reading accounts % for the special structure of the input file. % if Iread == 5 %FileName = 'cpR6400.00'; % Name of input file FileName = 'cpcent.00'; % Name of input file %FileName = 'cpcent.01'; % Name of input file Nstorm = 12; % number of sub-series Ntap = 18; % No of taps length = 1*4096; % number of time steps per storm SignalNo = 5; % Number of signal to be analysed (1-18) disp(['Data read from file:',FileName]) Series = zeros(Nstorm*length,Ntap); fid = fopen(FileName,'r'); % Echo print on screen index=0; for istorm = 1:12 f1 = (istorm-1)*length+1; f2 = istorm*length; qhmwk(istorm) = fscanf(fid,'%e',[1 1]); % velocity pressure [kN/m^2] cp = fscanf(fid,'%e',[18 length]); Series(f1:f2,:) = cp'; % saving the 12 data sets as continuous times series fprintf(1,'Storm Number considered: %g %g\n',istorm,qhmwk(istorm)) end status = fclose(fid); [m1 n1] = size(Series); % Number of rows (m1) and columns (n1) X0 = Series(:,SignalNo); % Saving selected data to input vector NFFTcase = 10*1024; % Filter depth of the fft-routine N*1024 Nwindow = 8; % Number of windows (default = 8) for pwelch SFD calculation Fsamp = 6400; % Sample frequency in DT = 1/Fsamp; % Calculation of time step Nbin = 100; % Number of bins to generate histogram Nsub = 1; % Number of sub-series for i=1:m1 % Generation of time axis with m1 steps TAx(i) = (i-1)*DT; end end if Iread == 6 FileName = 'test1.txt'; % Name of input file (Base bending moment CDJ) fid = fopen(FileName,'r'); Series = fscanf(fid,'%e',[12 inf]); % Matrix with time series Series = Series'; % Transposed matrix [m1 n1] = size(Series); % Number of rows (m1) and columns (n1) status = fclose(fid); disp(['Data read from file:',FileName]) TAx = Series(:,1); % Time series X0 = Series(:,2); % Saving selected data to variable vector SignalNo = 1; % Default value since there is only one signal to analyse NFFTcase = 5*1024; % Filter depth of the fft-routine N*1024 % NFFTcase = 512; % Filter depth of the fft-routine N*1024 Nwindow = 8; % Number of windows (default = 8) for pwelch SFD calculation Fsamp = 10000; % Sample frequency in DT = 1/Fsamp; % Calculation of time step Nbin = 100; % Number of bins to generate histogram Nsub = 12; % Number of sub-series for i=1:m1 % Generation of time axis with m1 steps TAx(i) = (i-1)*DT; end end %======================================================================================== % NOTE: At this point in the program you should have following information available: % % SignalNo = Number of signal that has been chosen to be analysed - saved as X0(i) % X0(i) = Vector with data of stochastic process (length = m1) % TAx(i) = Vector with values for time axis (length = m1) % m1 = number of data points (time steps) % NFFTcase = FFT filter length, determine by trial, shall not exceed m1 % Nwindow = Number of windows (default = 8) for pwelch SFD calculation % Fsamp = sample frequency in [Hz] % DT = time step between data points [s] = 1/Fsamp % Nbin = Number of bins to generate histogram % Nsub = Number of sub-series in which the signal can be divided % to calculate sub-mean and rms-values % % 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: % ----------------- % 2.2) Parameter Definitions: % The setting of the parameter switches different options for the analysis on % and off. The detailed setting for the different actions are defined in % section 1.4. <0> = no action % <1> = action activated % % The program allows for following data modification and analysis: % % 1. DETREND % 2. STATSITSICS part 1 % 3. SPECTRAL DENSITY part 1 % 4. DIGITAL FILTERING % 5. STATSITSICS part 2 % 6. SPECTRAL DENSITY part 2 % DoAct1 = 0; % Linear detrending the time series (no break points) DoAct2 = 1; % Identify, display and safe sub-series maxima and minima in vector DoAct3 = 0; % Digital filtering DoAct4 = 0; % Saving modified data in external file "Data2.dat" % Adjustments for display % Displ = 1 displaying modified data, if applied, and related results (set as default) % Displ = 2 displaying both initial and modified data and results for comparison Displ = 2; % Display parameter % 2.2) Parameter Definitions: % ---------------------- % General: % -------- pi = 4*atan(1.); % Circular constant % Spectral Density: % ----------------- % The spectral density is calculated using Welch's method in following format: % % [Pxx,f] = pwelch(x,window,noverlap,nfft,fs) % % Description (from the Matlab Desktop Help): % ------------------------------------------- % [Pxx,w] = pwelch(x) estimates the power spectral density Pxx of the input signal % vector x using Welch's method. Welch's method splits the data into overlapping segments, % computes modified periodograms of the overlapping segments, and averages the resulting % periodograms to produce the power spectral density estimate. % % - The vector x is segmented into eight sections of equal length, each with 50% overlap. % - Any remaining (trailing) entries in x that cannot be included in the eight segments of % equal length are discarded. % - Each segment is windowed with a Hamming window (see hamming) that is the same length % - as the segment. % % [Pxx,f] = pwelch(x,window,noverlap,nfft,fs) uses the sampling frequency fs specified in % hertz (Hz) to compute the PSD vector (Pxx) and the corresponding vector of frequencies (f). % In this case, the units for the frequency vector are in Hz. The spectral density produced % is calculated in units of power per Hz. If you specify fs as the empty vector [], the % sampling frequency defaults to 1 Hz. Nw = Nwindow; % Number of windows (default = 8) window = floor(m1*2/(1+Nw)); % length of windows assumed 50% overlap (still automatic default) nfft = NFFTcase; % Filter depth (should not exceed numer of time steps!) fs = Fsamp; % Sample frequency [Hz] % Digital filtering of the signal: % -------------------------------- % Fn = order of the filter (using standard 6th-order Butterworth filter) % Ftype = 'high' for a highpass digital filter with cutoff frequency CutOff % (high-frequency signals pass) % Ftype = 'low' for a lowpass digital filter with cutoff frequency CutOff % (low-frequency signals pass) % CutOff = Frequency [Hz] below/above which the frequency content will be filtered out Fn = 6; Ftype ='low'; CutOff = 10; %======================================================================================== % 3) DATA PROCESSING: % ---------------- % 3.1 Detrending unmodified X0(t): % -------------------------------- % The function "detrend" removes the mean value or linear trend from a vector or matrix. % % detrend(x,'linear') - removing of a linear trend % detrend(x,'linear',bp) - removing of linear trend between break points % % A breakpoint between two segments is defined as the data point that the two segments share. % The break points are given in vector "bp". if DoAct1 == 1 X1 = detrend(X0,'linear'); % no breakpoints defined else X1 = X0; end % 3.2 Digital Filtering of X1(t): % ------------------------------- % Calculation of unfiltered spectrum for later comparison X1d = detrend(X1,'constant'); % Removing mean value before performing FFT [Pxx,f] = periodogram(X1d,[],nfft,fs); % Periodogram [Sxx,f] = pwelch(X1d,window,[],nfft,fs); % Spectral density using Welch method Df = (f(10)-f(1))/9; % Frequency resolution of calculated spectrum XvarG = trapz(f(:),Sxx(:,1)); % Area underneath calcultaed SDF-curve (geometrical variance) Xvar = (std(X1d))^2; % statistical variance SNx1(:,1) = Sxx(:,1).*f./XvarG; % Ordinate normalised Spectral Curves SNx1(:,2) = Sxx(:,1)./XvarG; % Normalizing Spectral Curves to Unit Area (SNxx=Sxx/XvarG) SNx1(:,3) = Sxx(:,1)./XvarG*Xvar; % Spectrum with statistical variance as area SNx1(:,4) = Sxx(:,1); % Spectrum as calculated with pwelch (geometric variance underneath) % Performance of digital filtering if DoAct3 == 1 Nyquist = Fsamp/2; [b,a] = butter(Fn,CutOff/Nyquist,Ftype); Xf = filtfilt(b,a,X1); X2 = Xf; else X2 = X1; end % The stochastic process is now saved as X2(t) on which all subsequent % analysis will be peformed. % % 3.3 Statistical Parameter of Time History (X2): % ----------------------------------------------- % Echo print on print is default Xmean = mean(X2); % mean value of the parent time history Xstd = std(X2); % standard deviation of the parent time history Xvar = Xstd^2; % corresponding variance Xmax = max(X2); % maximum peak value occuring in parent time history Xmin = min(X2); % corresponding minimum peak value Tend = TAx(m1); % duration of parent time history % 3.4 Statistics on Sub-Series: % ----------------------------- % if DoAct2==0; Nsub=1; end Lsub = floor(m1/Nsub); % Length of sub-series in number of time steps Smean = zeros(Nsub,1); Srms = zeros(Nsub,1); Sx = zeros(Nsub,1); St = zeros(Nsub,1); Smin = zeros(Nsub,1); Smax = zeros(Nsub,1); SXmin = zeros(Nsub,1); SXmax = zeros(Nsub,1); for k=1:Nsub Smin(k) = +10E10; %Value to compare with data points in sub-series to identify sub-minimum. Smax(k) = -10E10; %Value to compare with data points in sub-series to identify sub-maximum. k1 = 1+(k-1)*Lsub; % Step number where sub-series starts k2 = k*Lsub; % Step number where sub-series ends Smean(k)= mean(X2(k1:k2)); % mean value of sub-series Srms(k) = std(X2(k1:k2)); % standard deviation of sub-series if DoAct2==1 for j=k1:k2 if X2(j)<=Smin(k) ; Smin(k)=X2(j) ; SXmin(k) = TAx(j) ; end if X2(j)>=Smax(k) ; Smax(k)=X2(j) ; SXmax(k) = TAx(j) ; end end end Sx(k) = k; % Number of sub-series St(k) = TAx(k*(Lsub-1)); % Location of sub-series boundaries on time axis end % 3.5 Calculating a histogram and probability density on detrended data: % ---------------------------------------------------------------------- RangeX = (Xmax-Xmin)*1.03; % Expanding range about 3% BinLow = Xmin-0.03*(Xmax-Xmin)/2; % Lower start point for bin grid DBin = RangeX/Nbin; % Bin width % Generating vector with Nbin+1 bin boundaries: BIN = zeros(Nbin+1,1); BIN(1) = BinLow; for i=1:Nbin BIN(i+1)=BinLow+i*DBin; end % Counting data points per bin: BinCount=zeros(Nbin,1); for i=1:Nbin BinL=BIN(i); BinU=BIN(i+1); for j=1:m1 if(X2(j)>BinL)&&(X2(j)<=BinU) BinCount(i)=BinCount(i)+1; end end end % Conversion to relative bin frequency RelFreq=zeros(Nbin,1); for i=1:Nbin RelFreq(i)=BinCount(i)/(m1*DBin); end Dx = DBin/10; i=0; for x = Xmin:Dx:Xmax i=i+1; pdf(i,1) = x; pdf(i,2) = 1/(Xstd*sqrt(2*pi))*exp(-0.5*((x-Xmean)/Xstd)^2); end % % Printing basic parameter of the analysis % ---------------------------------------- fprintf(1,'TIME HISTORY of Variable X\n'); fprintf(1,' Number of time steps in time history: %10.0f [-]\n',m1); fprintf(1,' Duration of parent time history : %10.2f [s]\n',Tend); fprintf(1,' Number of sub-series : %10.4g [s]\n',Nsub); fprintf(1,' Duration of sub-series : %10.2f [s]\n',Tend/Nsub); fprintf(1,' Time step width DT : %10.4g [s]\n',DT); fprintf(1,' Sample frequency (if [T]=s) : %10.4g [Hz]\n',Fsamp); fprintf(1,' Mean value of X(t) : %10.4g [x]\n',Xmean); fprintf(1,' Standard deviation of X(t) : %10.4g [x]\n',Xstd); fprintf(1,' Corresponding variance : %10.4g [x^2]\n',Xvar); fprintf(1,' Maximum peak value in X(t) : %10.4g [x]\n',Xmax); fprintf(1,' Minimum peak value in X(t) : %10.4g [x]\n',Xmin); fprintf(1,' \n'); fprintf(1,'SPECTRAL DENSITY parameters for Sxx\n'); fprintf(1,' Number of overlapping sub-windows : %10.0f [-]\n',Nw); fprintf(1,' Sub-window length : %10.0f [-]\n',window); fprintf(1,' Filter depth of fft-routine : %10.0f [-]\n',nfft); fprintf(1,' \n'); fprintf(1,' \n'); % 3.6 Spectral Density of X(t): % ----------------------------- X2d = detrend(X2,'constant'); % Removing mean value before performing FFT [Pxx,f] = periodogram(X2d,[],nfft,fs); % Periodogram [Sxx,f] = pwelch(X2d,window,[],nfft,fs); XvarG = trapz(f(:),Sxx(:,1)); % Area underneath calcultaed SDF-curve (geometrical variance) SNx2(:,1) = Sxx(:,1).*f./XvarG; % Ordinate normalised Spectral Curves SNx2(:,2) = Sxx(:,1)./XvarG; % Normalizing Spectral Curves to Unit Area (SNxx=Sxx/XvarG) SNx2(:,3) = Sxx(:,1)./XvarG*Xvar; % Spectrum with statistical variance as area SNx2(:,4) = Sxx(:,1); % Spectrum as calculated with pwelch (geometric variance underneath) % % 3.7 Saving modified data to external file: % ------------------------------------------ if DoAct4 ==1 fid1 = fopen('Data2.dat','w'); for i=1:m1 fprintf(fid1,' %g %g\n',TAx(i),X2(i)); end status = fclose(fid1); end %======================================================================================== % 4) GRAPHICAL DISPLAY OF THE EXTREME VALUE ANALYSIS: %---------------------------------------------------- % Display Definitions % ------------------- scrsz = get(0,'ScreenSize'); %==================================================================================== figure('Name','Time History','Position',[5 0.50*scrsz(4) 0.7*scrsz(3) 0.45*scrsz(4)]) % Setting axis range Y1max = Xmin+(Xmax-Xmin)*1.2; Y1min = Xmin-(Xmax-Xmin)*0.02; DY = Y1max-Xmax; X1min = TAx(1); X1max = TAx(m1); if Displ == 2 plot(TAx,X1,'-c'); hold on end if DoAct2==1 for j=1:Nsub plot(SXmax(j),Smax(j),'o','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',5);hold on plot(SXmin(j),Smin(j),'o','MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',5);hold on end end plot(TAx,X2); hold on plot3([TAx(1),TAx(m1)],[Xmean,Xmean],[1,1],'--m'); hold on plot3([TAx(1),TAx(m1)],[Xmean+Xstd,Xmean+Xstd],[1,1],'--g'); hold on plot3([TAx(1),TAx(m1)],[Xmean-Xstd,Xmean-Xstd],[1,1],'--g'); hold on plot([0 0],[Y1min Y1max],':k'); for j=1:Nsub plot([St(j) St(j)],[Y1min Y1max],':k'); end text(Tend/20,Y1max-0.3*DY,['Time History File is "',FileName,'"',' ; Signal No.:',num2str(SignalNo)],'FontSize',9) text(Tend/20,Y1max-0.7*DY,['Time History duration is ',num2str(Tend),' sec'],'FontSize',9) text(0.98*Tend,Y1max-0.3*DY,['Mean value of X(t): ',num2str(Xmean)],'FontSize',9,'HorizontalAlignment','right') text(0.98*Tend,Y1max-0.7*DY,['Standard deviation of X(t): ',num2str(Xstd)],'FontSize',9,'HorizontalAlignment','right') xlabel('time [s]'); ylabel('ordinate of variable X(t)'); title('TIME HISTORY ANALYSIS'); axis([X1min X1max Y1min Y1max]); %eval(['print -dtiff -zbuffer TimeHist']); %==================================================================================== figure('Name','Histogram','Position',[0.715*scrsz(3) 0.50*scrsz(4) 0.29*scrsz(3) 0.45*scrsz(4)]) XX=1.05*max(RelFreq); for i=1:Nbin plot([0 RelFreq(i) RelFreq(i) 0],[BIN(i) BIN(i) BIN(i+1) BIN(i+1)],'k'); % area([0 RelFreq(i) RelFreq(i) 0],[BIN(i) BIN(i) BIN(i+1) BIN(i+1)],'FaceColor',[.7 0 0]); hold on end plot([0 XX],[Xmean Xmean],'--m'); hold on plot([0 XX],[Xmean+Xstd Xmean+Xstd],'--g'); hold on plot([0 XX],[Xmean-Xstd Xmean-Xstd],'--g'); hold on plot(pdf(:,2),pdf(:,1),'-b','LineWidth',2); hold on title('HISTOGRAM'); xlabel('relative frequency'); ylabel('Data value'); axis([0 XX Y1min Y1max]); %eval(['print -dtiff -zbuffer Histogram']); %==================================================================================== figure('Name','Spectral Density','Position',[5 35 0.33*scrsz(3) 0.395*scrsz(4)]) if Displ == 2 loglog(f,SNx1(:,1),'.','MarkerSize',5,'Color','c'); hold on end loglog(f,SNx2(:,1),'.','MarkerSize',5,'Color','b'); hold on title('SPECTRAL DENSITY of X(t)'); xlabel('frequency [Hz]'); ylabel('Normalised Spectrum S_x_x(f) * f / \sigma_x^2'); grid on %eval(['print -dtiff -zbuffer SDF']); %==================================================================================== figure('Name','Spectral Density 2','Position',[5 35 0.33*scrsz(3) 0.395*scrsz(4)]) loglog(f,SNx2(:,2),'.','MarkerSize',5,'Color','b'); hold on title('SPECTRAL DENSITY of X(t)'); xlabel('frequency [Hz]'); ylabel('Normalised Spectrum S_x_x(f) / \sigma_x^2'); grid on %eval(['print -dtiff -zbuffer SDF2']); %==================================================================================== figure('Name','Spectral Density 3','Position',[5 35 0.33*scrsz(3) 0.395*scrsz(4)]) loglog(f,SNx2(:,3),'.','MarkerSize',5,'Color','b'); hold on title('SPECTRAL DENSITY of X(t)'); xlabel('frequency [Hz]'); ylabel('Spectrum S_x_x(f)'); grid on %eval(['print -dtiff -zbuffer SDF3']); %==================================================================================== figure('Name','Periodogram','Position',[5 35 0.33*scrsz(3) 0.395*scrsz(4)]) loglog(f,Pxx,'.','MarkerSize',5,'Color','b'); hold on title('PERIODOGRAM of X(t)'); xlabel('frequency [Hz]'); ylabel('Spectrum S_x_x(f)'); grid on %eval(['print -dtiff -zbuffer SDF4']); %==================================================================================== figure('Name','SubSeries Parameters','Position',[0.34*scrsz(3) 35 0.33*scrsz(3) 0.395*scrsz(4)]) plot(Sx,Smean,'s','MarkerEdgeColor','k','MarkerFaceColor','m','MarkerSize',7); hold on plot(Sx,Srms,'o','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',7); hold on title('SUB-SERIES PARAMETERS'); xlabel('number of sub-series'); ylabel('Mean and rms value'); legend('mean','rms','Location','Best'); plot([0 Nsub],[Xmean Xmean],'--m'); hold on plot([0 Nsub],[Xstd Xstd],'--g'); hold on grid on %eval(['print -dtiff -zbuffer SubSeriesParam']);