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