classdef Signal
    %SIGNAL Summary of this class goes here
    %   Detailed explanation goes here

    properties
        signal
        logbook
        fs
    end

    methods
        function obj = Signal(signal,options)
            %SIGNAL Construct an instance of this class
            %   Detailed explanation goes here
            arguments
                signal
                options.fs = [];
            end

            obj.signal = signal;
            obj.signal = obj.signal;
            obj.fs = options.fs;

            SignalType  = [];
            TimeStamp   = [];
            Length      = [];
            SignalPower = [];
            Nase        = [];
            Description = [];

            obj.logbook = table(SignalType,TimeStamp,Length,SignalPower,Nase,Description);


        end

        %% CONVERT TO INFORMATIONSIGNAL
        function [i_sig, varargout] = Informationsignal(obj,options)
            arguments
                obj
                options.fs
                options.logbook
            end

            if isa(obj,'Electricalsignal')

                %convert to information
                varargout{1} = obj.fs;
                i_sig = Informationsignal(obj.signal,"fs",options.fs,"logbook",options.logbook);

            elseif isa(obj,'Opticalsignal')

                error("Cannot convert from optical- to informationsignal. Use O/E conversion first.");

            end

        end

        %% CONVERT TO Electricalsignal
        function [e_sig, varargout] = Electricalsignal(obj,options)

            arguments
                obj
                options.fs
                options.logbook
            end

            obj.logbook = options.logbook;

            if isa(obj,'Opticalsignal')

                %convert to electrical
                varargout{1} = obj.nase;
                varargout{2} = obj.lambda;
                e_sig = Electricalsignal(obj.signal,"fs",obj.fs,"logbook",obj.logbook);

            elseif isa(obj,'Informationsignal')

                try
                    % specify fs at varargin{1}
                    e_sig = Electricalsignal(obj.signal,"fs",options.fs,"logbook",options.logbook);
                catch
                    error("Signal Conversion failed [I -> E] ");
                end

            end

        end


        %% CONVERT TO Opticalsignal
        function o_sig = Opticalsignal(obj, options)

            arguments
                obj
                options.fs
                options.logbook
                options.nase
                options.lambda
            end

            fn = fieldnames(options);

            for l = 1:numel(fn)
                try
                    obj.(fn{l}) = options.(fn{l});
                end
            end

            if isa(obj,'Electricalsignal')

                %convert to optical

                o_sig = Opticalsignal(obj.signal,"fs",obj.fs,"lambda",options.lambda,"logbook",obj.logbook,"nase",options.nase);


            elseif isa(obj,'Informationsignal')

                error("Cannot convert from information- to opticalsignal. Use E/O conversion first.");

            end

        end

        function plot(obj, options)
            % signal to plot: obj.signal
            % fsamp : obj.fs (e.g. 92e9 => 92 GHz)
            % length: length(obj.signal)

            arguments
                obj
                options.fignum
                options.displayname = [];
                options.timeframe = 0;
            end

            figure(options.fignum);  % If figure does not exist, create new figure

            % 2) Plot into the figure handle found or created in one
            t = (0:length(obj.signal)-1) / obj.fs;  % time vector

            if options.timeframe ~= 0
                %only show a certain timeframe of signal
                t = t(t<options.timeframe);
            end
            % 2 a) Actual plot (hold on, displayname??)

            dn = options.displayname;

            if isa(obj,'Opticalsignal')
                sig = abs(obj.signal).^2;
            else
                sig = obj.signal;
            end

            hold on;
            plot(t, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);

            % 2 c)
            % - xlabel if not already here: time in readable format (1 ms and not 1e-3 s)
            % - ylabel amplitude
            if isempty(get(gca, 'XLabel').String)
                xlabel('Time (mu s)');
            end

            if isempty(get(gca, 'YLabel').String)
                ylabel('Amplitude');
            end

            % Convert time axis to milliseconds for readability
            xticks = get(gca, 'XTick');
            set(gca, 'XTick', xticks, 'XTickLabel', xticks * 1e6);

            % Add legend if not already present
            if isempty(get(gca, 'Legend'))
                legend;
            end

            hold off;

        end

        %% Add signals from one signal to another, the first object will sustain
        function Sum = plus(X,y)

            if  isa(y,'Signal')
                Sum = X;
                Sum.signal = X.signal + y.signal;
            elseif isnumeric(y)
                Sum = X;
                Sum.signal = X.signal + y;
            end

        end

        function Product = times(X,y)

            if  isa(y,'Signal')
                Product = X;
                Product.signal = X.signal .* y.signal;
            elseif isnumeric(y)
                Product = X;
                Product.signal = X.signal .* y;
            end

        end

        %% Display length
        function return_length = length(obj)
            %METHOD1 Summary of this method goes here
            %   Detailed explanation goes here
            return_length = length(obj.signal);
        end

        %% Write Logbook Entry
        function obj = logbookentry(obj,varargin)

            if nargin > 1
                Description = varargin{1};
            else
                Description = "";
            end

            SignalType  = [string(class(obj))];
            TimeStamp   = [(datetime('now','TimeZone','local','Format','HH:mm:ss'))];
            Length      = num2str(obj.length, ['%' sprintf('.%df', 0)]);%[obj.length];
            SignalPower = [obj.power];
            Nase        = [0];

            cell = {SignalType , TimeStamp , Length , SignalPower(1) , Nase, Description};

            obj.logbook = [obj.logbook;cell];

        end

        %% Resample Signal
        function obj = resample(obj,options)

            arguments
                obj Signal
                options.fs_in double = obj.fs
                options.fs_out double
                options.n double = 10;
                options.beta double = 5;
            end

            if options.fs_in ~= obj.fs
                warning('The signals fs is different from the given fs_in while it should be the same.');
            end

            obj.signal  = resample(obj.signal,options.fs_out,options.fs_in,options.n,options.beta);

            desc = ['resample signal from ', num2str(options.fs_in*1e-9), ' GHz to ', num2str(options.fs_out*1e-9), ' GHz' ];

            obj = obj.logbookentry(desc);

            obj.fs = options.fs_out;

        end

        %%
        function spectrum(obj,options)

            arguments
                obj
                options.fignum
                options.displayname = "";
            end

            % spectrum_plot(obj.signal,options.fsamp,options.figurename,options.displayname);

            N = 2^(nextpow2(length(obj.signal))-8);

            [p_lin,w] = pwelch(obj.signal,hanning(N),N/2,N,obj.fs,"centered","power","mean");
            p_dbm = 10*log10(p_lin)+30; %dB to dBm in case of "power"

            figure(options.fignum);  % If figure does not exist, create new figure
            hold on
            plot(w.*1e-9,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
            xlabel("Frequency in GHz");
            %ylabel("Power/frequency (dB/Hz)");
            ylabel("Power (dBm)");
            xlim([-obj.fs/2 obj.fs/2].*1e-9)
            edgetick = 2^(nextpow2(obj.fs*1e-9));
            % xticks([-edgetick:16:edgetick]);
            xlim([-244, 244])
            ylim([-120,10]);
            yticks([-200:10:10]);
            legend

        end

        %% Power of signal
        function pow = power(obj,options)

            arguments
                obj
                options.unit power_notation = power_notation.dBm
            end

            pow = mean(abs(obj.signal).^2);

            switch options.unit
                case power_notation.dBm
                    if isa(obj,'Electricalsignal')
                        pow = pow / 50;
                    end
                    pow = 10*log10(pow)+30; %dbm

                case power_notation.mW
                    pow = pow .* 1e3; %mW

                case power_notation.W
                    %pow = pow % Watt

            end

        end

        %% Peak Power of Signal
        function pow_pk = power_peak(obj,options)

            arguments
                obj
                options.unit power_notation = power_notation.dBm
            end

            pow_pk = max(abs(obj.signal).^2); % dBm

            switch options.unit
                case power_notation.dBm
                    pow_pk = pow2db(pow_pk)+30; %dbm
                case power_notation.mW
                    pow_pk = pow_pk .* 1e3; %mW
                case power_notation.W
                    pow_pk = pow_pk; % Watt
            end

        end

        %% PAPR of signal
        function papr = papr_lin(obj)
            %PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power)
            % divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor.
            % papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2;

            papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W); %linear

        end

        %% PAPR of signal
        function papr_db = papr_db(obj)
            %PAPR The peak-to-average power ratio (PAPR) is the peak amplitude squared (giving the peak power)
            % divided by the RMS value squared (giving the average power).[1] It is the square of the crest factor.
            % papr = max(abs(timesignal))^2 / rms(timesignal)^2; ODER papr = peak2rms(sig)^2;

            papr = obj.power_peak("unit",power_notation.W) / obj.power("unit",power_notation.W);
            papr_db = 10*log10(papr);

            average_power = obj.power;
            peak_power = obj.power_peak;
            papr_db = peak_power - average_power;  %db

        end

        %%
        function obj = normalize(obj,options)

            arguments
                obj Signal
                options.mode normalization_mode = normalization_mode.rms
            end

            switch options.mode
                case normalization_mode.rms
                    obj.signal = obj.signal/sqrt(mean(abs(obj.signal).^2,"all"));
                case normalization_mode.oneone
                    obj.signal = obj.signal - min(obj.signal);
                    obj.signal = obj.signal/max(abs(obj.signal));
                    obj.signal = (2*obj.signal) - 1;
            end

        end

        %%
        function [obj] = delay(obj,delay,options)

            arguments
                obj Signal
                delay double = 0
                options.mode delay_mode = delay_mode.samples
            end

            if options.mode == delay_mode.samples

                obj.signal=delayseq(obj.signal,delay);

            elseif options.mode == delay_mode.time

                obj.signal=delayseq(obj.signal,delay,obj.fs);

            end

        end

        %%
        function [obj,S] = tsynch(obj,options)
            % time sync and cut
            arguments
                obj Signal
                options.reference Signal
                options.fs_ref = 0;
            end

            %normalize the signal
            a = obj.normalize("mode","oneone").signal;

            %resample the reference
            q = obj.fs/options.fs_ref;
            b = options.reference.resample("fs_in",options.fs_ref,"fs_out",obj.fs).normalize("mode","oneone").signal;

            %estimate delay between signals
            [co,lags] = xcorr(a,b);
            [~,pos] = max(co);
            D = lags(pos);

            %estimate start pos of signal
            maxpeaknum = floor(length(a)/length(b));
            [pks,pkpos] = findpeaks(co./max(co),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum);
            shifts = lags(pkpos);
            
            %Cut occurences of ref signal from signal
            S = {};
            for c = shifts
                sig = obj.delay(-c,'mode','samples');
                sig.signal = sig.signal(1:length(b));
                S{end+1,1} = sig;
            end
            
            %plot all synced signals and the ref signal
            figure;hold on;
            for i = 1:size(S,1)
                plot(S{i}.normalize('mode','oneone').signal(1000:1100),'LineWidth',0.1,'Color',[0.2157    0.4941    0.7216]);
                plot(b(1000:1100),'LineWidth',1);
            end
            
            %return the sinal with the highest correlation...
            [~,idx]=max(pks);
            obj.signal = S{idx}.signal;


        end

        function er = extinctionratio(obj,fsym,M)
            histpoints = 1024;                                              %% verticale resolution
            histpoints = floor(histpoints/2)*2+1;                           %% to have the eye digram centered around one point make the vertical resolution uneven
            histpoints_horizontal = 512;                                    %% horizontal resolution
            hist_data=zeros(histpoints,histpoints_horizontal );             %% initilize eye diagram

            if isa(obj,'Opticalsignal')
                sig = abs(obj.signal).^2;
            elseif isa(obj,'Electricalsignal')
                sig = obj.signal;
            else
                sig = obj.signal;
            end

            x = (sig);                                               %% make input signal rea)l

            x = resample(x,fsym*histpoints_horizontal/2,obj.fs);         %% up sample to original fsym rate

            if mod(length(x),2)==1                                          %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol
                x = x(1:end-1);
            end

            eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols)

            maxA = max(sig(100:end-100));
            minA = min(sig(100:end-100));

            difference= maxA-minA;

            data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1;

            for n=1:size(data_ind_y,1)
                nn=histcounts(data_ind_y(n,:),1:histpoints+1);
                hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-(
            end

            plot_data = 20*log10(hist_data);
            plot_data(plot_data==-Inf) = 0;

            maxall = 0;

            for l = 1:size(plot_data,2)
                [maxpk_,pos_] = max(plot_data(:,l));
                if maxpk_ > maxall
                    maxall = maxpk_;
                    posxall = l;
                    posyall = pos_;
                end
            end

            hist_interest = plot_data(:,posxall);
            hist_interest_smoth = smooth(hist_interest,20);
            [pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);

            for i = 1:numel(loc)
                ppeak(i) = maxA - (difference/histpoints*loc(i));
            end



            if isa(obj,'Opticalsignal')
                er=10*log10(ppeak(1)/ppeak(end));
            elseif isa(obj,'Electricalsignal')
                if  mean([ppeak(1),ppeak(end)]) < 1e-2
                    disp("No Extiction Ration for Bipolar Electrical Signal. Calculating Outer OMA instead...")
                    er=max(ppeak)-min(ppeak);
                else
                    er=10*log10(ppeak(1)/ppeak(end));
                end
            else
                er=10*log10(ppeak(1)/ppeak(end));
            end

            if 0
                findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
            end


        end

        function eye(obj,fsym,M)

            mode = 1;

            histpoints = 1024;                                              %% verticale resolution
            histpoints = floor(histpoints/2)*2+1;                           %% to have the eye digram centered around one point make the vertical resolution uneven
            histpoints_horizontal = 512;                                    %% horizontal resolution
            hist_data=zeros(histpoints,histpoints_horizontal );             %% initilize eye diagram

            if isa(obj,'Opticalsignal')
                sig = abs(obj.signal).^2;
            elseif isa(obj,'Electricalsignal')
                sig = obj.signal;
            else
                sig = obj.signal;
            end

            x = (sig);                                               %% make input signal rea)l

            x = resample(x,fsym*histpoints_horizontal/2,obj.fs);         %% up sample to original fsym rate

            if mod(length(x),2)==1                                          %% if the signal lenght is not divisible by 2 (symbols displayed in the eye diagram are 2) remove last symbol
                x = x(1:end-1);
            end

            eye_mat = reshape(x(1:end-mod(length(x),histpoints_horizontal)),histpoints_horizontal,floor(length(x)/histpoints_horizontal)); %% reshape signal into 256 rows each row has the histogram(eye data of all symbols)

            figure(922)
            clf
            if mode == 2
                % generate "intuitive eye diagram" by drawing lines on top over
                % each other; only draw 1000 lines, otherwise the plot is too
                % crowded


                col = cbrewer2('Set1',2);
                for n=1:1000
                    hold on
                    plot(eye_mat(:,n),'LineStyle',':','LineWidth',0.1,'Color',col(2,:));
                end
                ylabel('Amplitude of Signal');

            elseif mode == 1
                % generate eye diagram using histogram

                maxA = max(sig(100:end-100));
                minA = min(sig(100:end-100));

                % maxA = 0.0015;
                % minA = 0;



                difference= maxA-minA;

                data_ind_y=round((eye_mat-minA)/difference*(histpoints-1)) +1;

                for n=1:size(data_ind_y,1)
                    nn=histcounts(data_ind_y(n,:),1:histpoints+1);
                    hist_data(:,n)=flip(nn.'); %without flip, the eye is upside down :-(
                end

                plot_data = 20*log10(hist_data);
                plot_data(plot_data==-Inf) = 0;

                imagesc(plot_data);



                % beautify
                colormap(cbrewer2("Blues",4096));


                if isa(obj,'Opticalsignal')
                    title("Optical Eye")
                    ylabel("Power in mW");
                    y_tickstring = string(linspace(maxA.*1e3,minA.*1e3,16));
                    min_ = min(abs(obj.signal(100:end-100)).^2);
                    max_ = abs(max(obj.signal(100:end-100)).^2);
                elseif isa(obj,'Electricalsignal')
                    title("Electrical Eye")
                    ylabel("Voltage in V");
                    y_tickstring = string(linspace(maxA,minA,16));
                    min_ = min(obj.signal(100:end-100));
                    max_ = abs(max(obj.signal(100:end-100)));
                else
                    title("Digital Eye")
                    ylabel("Digital Signal Amplitude");
                    y_tickstring = string(linspace(maxA,minA,16));
                    min_ = min(obj.signal(100:end-100));
                    max_ = abs(max(obj.signal(100:end-100)));
                end
                xlabel('Time in ps')

                % add information

                if 1

                    pwr_dbm = round(obj.power,3);
                    pwr_lin = obj.power("unit",power_notation.W);

                    papr_ = obj.papr_lin;%round(papr(obj.signal.^2),3);

                    yline( histpoints-(pwr_lin - minA)/difference*histpoints );
                    yline( histpoints-(min_    - minA)/difference*histpoints );
                    yline( histpoints-(max_    - minA)/difference*histpoints );

                    maxall = 0;
                    for l = 1:size(plot_data,2)
                        [maxpk_,pos_] = max(plot_data(:,l));
                        if maxpk_ > maxall
                            maxall = maxpk_;
                            posxall = l;
                            posyall = pos_;
                        end
                    end

                    hold on
                    xline(posxall)



                    % Define properties
                    boxPosition = [0.15 0.86 0.2 0.05]; % Position for the first box [x y width height]
                    boxColor = [0.9 0.9 0.9]; % Light grey background color
                    boxEdgeColor = 'k'; % Black edge color
                    boxLineStyle = '--'; % Dashed line style
                    boxFontWeight = 'bold'; % Bold font

                    % Create first annotation box for Power
                    annotation('textbox', boxPosition, ...
                        'String', ['Power: ',num2str(pwr_dbm),' dBm'], ...
                        'BackgroundColor', boxColor, ...
                        'EdgeColor', boxEdgeColor, ...
                        'LineStyle', boxLineStyle, ...
                        'FontWeight', boxFontWeight, ...
                        'HorizontalAlignment', 'center');

                    % Adjust position for the second box (slightly to the right)
                    boxPosition = [0.37 0.86 0.2 0.05]; % Adjusted position

                    % Create second annotation box for PAPR
                    annotation('textbox', boxPosition, ...
                        'String', ['PAPR(lin):',num2str(papr_),''], ...
                        'BackgroundColor', boxColor, ...
                        'EdgeColor', boxEdgeColor, ...
                        'LineStyle', boxLineStyle, ...
                        'FontWeight', boxFontWeight, ...
                        'HorizontalAlignment', 'center');

                    try
                    hist_interest = plot_data(:,posxall);
                    hist_interest_smoth = smooth(hist_interest,20);
                    a = scatter(hist_interest_smoth+posxall,1:length(hist_interest_smoth),4,'.','MarkerEdgeColor','red');

                    [pk,loc] = findpeaks(hist_interest_smoth,"MinPeakDistance",40,"NPeaks",M,"MinPeakHeight",30);
                    scatter(posxall,loc,'red','Marker','x','LineWidth',2);
                    yline(loc,'Color','red','LineWidth',1,'LineStyle',':');

                    for i = 1:numel(loc)
                        ppeak(i) = maxA - (difference/histpoints*loc(i));
                    end

                    oma = false;
                    if isa(obj,'Opticalsignal')
                        er=10*log10(ppeak(1)/ppeak(end));
                    elseif isa(obj,'Electricalsignal')
                        if  mean([ppeak(1),ppeak(end)]) < 1e-2
                            oma = true;
                            er=max(ppeak)-min(ppeak);
                        else
                            er=10*log10(ppeak(1)/ppeak(end));
                        end
                    else
                        er=10*log10(ppeak(1)/ppeak(end));
                    end

                    % Adjust position for the third box (slightly to the right)
                    boxPosition = [0.59 0.86 0.2 0.05]; % Adjusted position

                    % Create third annotation box for Vmax
                    if ~oma
                        thirdboxstring = ['ER (db):',num2str(er),' dB'];
                    else
                        thirdboxstring = ['OMA outer:',num2str(er),' V'];
                    end

                    

                    annotation('textbox', boxPosition, ...
                        'String',thirdboxstring , ...
                        'BackgroundColor', boxColor, ...
                        'EdgeColor', boxEdgeColor, ...
                        'LineStyle', boxLineStyle, ...
                        'FontWeight', boxFontWeight, ...
                        'HorizontalAlignment', 'center');
                    end

                    yticks(linspace(0,histpoints,16));
                    y_tickstring = sprintfc('%.2f', y_tickstring);
                    yticklabels(y_tickstring);

                    xticks(linspace(0,histpoints_horizontal,8))
                    x_tickstring = sprintfc('%.2f', linspace(0, 2/fsym, 8) .* 1e12);
                    xticklabels(x_tickstring);
                end




            end

            % disp('h');
            %
            % fsig = obj.fs;
            % q = fsig/fsym;
            %
            % if q > 10 && isinteger(q)
            %     sig = (obj.signal);
            % else
            %     sig = (obj.resample("fs_in",fsig,"fs_out",fsym*30).signal);
            %     q = 10;
            % end
            %
            % figure()
            % clf
            % cursor = 200*q;
            %
            % for s = 1:700
            %     plot(sig(cursor-q:cursor+q),'Color','black','LineWidth',0.1,'LineStyle','-');
            %     hold on
            %     cursor=cursor+q;
            %     s=s+1;
            % end
            %
            % ylim([-3 3]);


        end

    end
end