classdef SystemIdentificationTlbx
    % Collection of Methods (Functions) for measuring the response (Impulse Response and the
    % Transfer function) of a LTI system
    
    properties
        IRcoefficients
        FRF
        time_vector
        frequency_vector
        fs
        N_fft
        
        % USEFUL PROPERTIES FOR PROCESSING
        % These properties are included to facilitate the access to
        % the data from external script. Example: in order to plot the
        % magnitude of the FRF use
        % >> plot(ObjectName.FRFMagnitude_dB)
        % instead of
        % >> plot(20*log10(abs(ObjectName.FRF)))
        
        COH; %Coherence of the estimation
        FRFMagnitude_dB % The 20*log10(abs(FRF))
        FRFPhase_rads % angle(FRF)
        FRFPhase_degs % rad2deg(angle(FRF))
        IR_MinimumPhasePart % The minimum phase part of the IR (see Get_Minimum_Phase method)
        IRpeakvector % The vector of the peaks of the IR
        
        % TO BE IMPLEMENTED
        % Group Delay Calculation [Gd,F] = grpdelay(SystemRespObj.IRcoefficients,1,SystemRespObj.N_fft,SystemRespObj.fs);
    end
    
    %%
    methods%Private Methods
        function obj = SystemResponseToolbox(InputSignal, OutputSignal, N_fft, fs, EstimatorType, varargin)
            % Returns a SYSTEMRESPONSE object (see systemresponse.m for
            % details of the properties and methods) after the estimation
            % of the System Response (performed via the HEstimator method)
            
            % Initialise the VARARGIN options
            AlignmentBoolean = 1; %perform the estimation with the alignment of the IN/OUT signals and then shift the IR back
            GetTheMinimumPhase = 0; %calculates the Minimum Phase component of the IR
            
            if (~isempty(varargin))
                for c = 1:length(varargin)
                    switch varargin{c}
                        case {'NoInputOutputAlignment'}
                            AlignmentBoolean = 0;
                            
                        otherwise
                            error(['Invalid optional argument, ', ...
                                varargin{c}]);
                    end % switch
                end % for each varargin entry
            end % if varargin is not empty
            
            if nargin < 3
                obj.N_fft = 2^nextpow2(length(h));
            else
                obj.N_fft = N_fft;
            end
            
            [h, H, coherence, f, t] = SystemResponseToolbox.HEstimator(InputSignal, OutputSignal, N_fft, fs, EstimatorType, AlignmentBoolean);
            %obj = SystemResponse(h, H, t, f, fs, N_fft);% Creates a new SYSTEMRESPONSE object with the new estimated data
            
            obj.IRcoefficients = h;%the coefficients of the Impulse response
            
            if GetTheMinimumPhase
                obj.IR_MinimumPhasePart = Get_Minimum_Phase(h, N_fft); %Calculates the Minimu Phase component of the Impulse response
            end
            
            obj.time_vector = transpose(t);
            
            obj.FRF = H; %The frequency response function (FRF)
            obj.FRFMagnitude_dB = ConvertIndB(H);
            obj.FRFPhase_rads = angle(H);
            obj.FRFPhase_degs = rad2deg(angle(H));
            
            obj.fs = fs;
            obj.frequency_vector = transpose(f);%(0:obj.N_fft/2)/fs;
            obj.COH = coherence;
        end
        
    end%%Private Methods
    
    methods(Static)
        
        
        function peakvector = findIRpeaks(coefficients, numberofpeaks)
            
            delta = 10^-3;
            
            [maxtab, mintab] = peakdet(coefficients, delta);
            
            peak_values = maxtab(:,2);
            peak_indices = maxtab(:,1);
            
            for peak_idx = 1:numberofpeaks
                
                [~, current_peak_idx] = max(peak_values);
                
                peakvector(peak_idx) = peak_indices(current_peak_idx);
                %set the maximum to a very low value so it will be excluded from the research
                peak_values(current_peak_idx) = 10^-20;
            end
        end%findIRpeaks        %%
        function Plot_Estimation_Results(obj)
            % Plotting the results of the HEstimator
            
            %Preset = get(handles.PreSetAxisChk,'value');
            figure('name','System Response Estimation');
            subplot(2,2,1); obj.PlotIR(obj); title('Estimated IR'); grid on;
            subplot(2,2,2); semilogx(obj.frequency_vector, obj.COH); title('Coherence'); xlim([50 (obj.fs)/2]); grid on;
            subplot(2,2,3); obj.Plot_Magnitude_FRF(obj); title('Estimated Magnitude FRF, dB');   xlim([50 (obj.fs)/2]); grid on;
            subplot(2,2,4); obj.Plot_Phase_FRF(obj); title('Estimated Phase FRF, rads');  xlim([50 (obj.fs)/2]); grid on;
        end%ploteverything
        
        %% Methods for the IR
        function PlotIR(obj, x_axis_range)
            
            plot(obj.IRcoefficients);
            
            if nargin < 2
                xlim([0 length(obj.IRcoefficients) - 1]);
            else
                xlim([x_axis_range(1) x_axis_range(2)]);
            end
            
            xlabel('Samples'); ylabel('Impulse Response Amplitude');
        end%plotIR
        
        
        
        function PlotIRWithPeaks(obj, numberofpeaks, c0, x_axis_range)
            obj.IRpeakvector = findIRpeaks(obj.IRcoefficients, numberofpeaks);
            peak_values = obj.IRcoefficients(obj.IRpeakvector);
            x_coordinates = obj.IRpeakvector;
            y_coordinates = peak_values;
            
            plot(obj.IRcoefficients);
            for idx_peak = 2:numberofpeaks
                time = (x_coordinates(idx_peak) - x_coordinates(1))/obj.fs;
                space = c0*time;
                hold all; plot([x_coordinates(1) x_coordinates(idx_peak)], [y_coordinates(idx_peak) y_coordinates(idx_peak)]);
                text(x_coordinates(idx_peak), y_coordinates(idx_peak),[num2str(roundn(space,-2)) 'm']);
            end
            xlim([x_axis_range(1) x_axis_range(end)]);
            xlabel('Samples'); ylabel('Impulse Response Amplitude');
        end%plotIRpeaks
        
        %% Methods for the FRF
        function Plot_Magnitude_FRF(obj)
            semilogx(obj.frequency_vector, obj.FRFMagnitude_dB);
            xlim([50 obj.fs/2]);
            xlabel('Frequency, Hz'); ylabel('Magnitude FRF');
        end%Plot_Magnitude_FRF
        
        function Plot_Phase_FRF(obj, UnwrappedBoolean)
            if nargin < 2
                UnwrappedBoolean = 0;
            end
            
            string1 = 'semilogx(obj.frequency_vector, ';
            
            if UnwrappedBoolean
                string2 = 'unwrap(obj.FRFPhase_rads)';
            else
                string2 = 'obj.FRFPhase_rads';
                
            end
            
            string3 = ' );';% closes the semilogx function in string1
            eval([string1 string2 string3]);
            xlim([50 obj.fs/2]);
            xlabel('Frequency, Hz'); ylabel('Phase FRF, rads');
        end%Plot_Phase_FRF
        
        function Plots(MethodStr, h, H, fs, COH, FigureNameTextStr)
            %PLOTS Summary of this function goes here
            % MethodStr can be 'SineSweep' 'HEstimator'
            if nargin < 6
            else
                figure('Name', [FigureNameTextStr ' WhiteNoise']);
            end
            
            N_fft = length(h);
            
            switch MethodStr
                case 'HEstimator'
                    NumberOfPlots = 4; % IR, Magnitude FRF, Phase FRF, COH
                case  'SineSweep'
                    NumberOfPlots = 3; % Same as above but it does not include Coherence
            end
            
            
            samplevect = (0:(N_fft - 1));
            timevect = samplevect/fs;
            freqvect = linspace(0, fs/2, N_fft/2 + 1);
            
            % IR
            subplot(NumberOfPlots, 1, 1); 
            plot(samplevect, h); title('Estimated IR'); xlabel('Samples'); ylabel('Amplitude');
            xlim([samplevect(1), samplevect(end)]);
            
            % Code to plot a double x-axis (samples and time)
%             Ax1 = gca;
%             Ax2 = axes('Position', get(Ax1,'Position'), 'XAxisLocation','top');
%             plot(timevect, h, 'color','k','parent', Ax2);
%             xlim([timevect(1), timevect(end)]);
             grid on;
            
             % Magnitude
            subplot(NumberOfPlots, 1, 2); semilogx(freqvect, LogDB(abs(H))); 
            title('Estimated Magnitude FRF'); xlabel('Freq, Hz'); ylabel('dB');
            grid on; xlim([50, fs/2]);
            
            % Phase
            subplot(NumberOfPlots, 1 ,3); 
            semilogx(freqvect, unwrap(angle((H)))); 
            title('Estimated FRF Unwrapped Phase'); xlabel('Freq, Hz');xlim([50, fs/2]);  grid on;
            
            if strcmp(MethodStr, 'HEstimator') % then plot the COH
            
            subplot(NumberOfPlots, 1, 4); semilogx(freqvect, COH); title('Coherence'); xlabel('Freq, Hz');
            grid on; xlim([50, fs/2]);
            end
        end
        
        
        
    end%Static Methods
end






function PlotStr = ReturnPlotString(AX, ThingToPlotAsAString, ColorString, Xlinlog,Unwrap,DegRad)
% Example ThingToPlotAsAString = 'abs(H)', or '20*log10(abs(H))', or 'angle(H)', or 'COH'

% Example ColorStr = 'r', or 'k'...
if nargin < 4
    
end

if nargin < 3
    %UnWrap = 0;
end

% Choosing the MATLAB functions to be used in the presentation of the results
% All the strings will be evaluated with the EVAL function

if Xlinlog == 1 %if linear scale
    plotting_function = 'plot'; %use function PLOT
else
    plotting_function = 'semilogx'; %use SEMILOGX
end

if Unwrap == 1
    unwrapping_function = 'unwrap';
else
    unwrapping_function = '';
end

if DegRad == 1
    degrees_function = 'rad2deg';
else
    degrees_function = '';
end

if isempty(ColorString)
    ColorStrTemp = '';
else
    ColorStrTemp = [' ''Color'', '' ' ColorString ' '','];
end

% Constructing the final string

PlotStr = [plotting_function '(' ...
    AX ', f,', ...
    degrees_function '(' unwrapping_function '(' ...
    ThingToPlotAsAString ')), ' ColorStrTemp  ' ''LineWidth'', 2);' ...
    ];
end