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

    properties
        signal
        logbook
        fs

        gitSHA
        gitStatus
        gitPatch
    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;

            [~,obj.gitSHA]      = system('git rev-parse HEAD');
            [~,obj.gitStatus]   = system('git status --porcelain');
            % [~,obj.gitPatch]    = system('git diff');

            %%% Stuff for Logbook %%%
            SignalType  = [];
            TimeStamp   = [];
            Length      = [];
            SignalPower = [];
            Nase        = [];
            SignalCopy  = [];
            ModifierName = [];
            ModifierCopy= {};
            Description = [];

            obj.logbook = table(SignalType,TimeStamp,Length,SignalPower,Nase,SignalCopy,ModifierName, ModifierCopy, 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;
                options.clear = 0;
                options.color = [];
            end

            figure(options.fignum);  % If figure does not exist, create new figure
            if options.clear
                clf
            end

            % 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;
            if isempty(options.color)
                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1);
            else
                plot(t* 1e6, sig(1:length(t)), 'DisplayName', dn, 'LineWidth', 0.1, 'Marker', '.', 'LineStyle','none', 'MarkerSize', 0.1,'Color',options.color);
            end
            % 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

           
            % 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

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

            if  isa(y,'Signal')
                Diff = X;
                Diff.signal = X.signal - y.signal;
            elseif isnumeric(y)
                Diff = X;
                Diff.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 == 2
                Description = varargin{1};
                CallingModifier = evalin('caller','obj');
            elseif nargin == 3
                Description = varargin{1};
                CallingModifier = varargin{2};
            else
                Description = "";
                CallingModifier = evalin('caller','obj');
            end

            if isa(CallingModifier,"Signal")
                CallingModifierStruct = obj.objToStructFilteredRecursive(CallingModifier);
                ModifierCopy = {CallingModifierStruct};
            end

            ModifierCopy = {CallingModifier};

            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];
            SignalCopy  = obj.signal;
            ModifierName = class(CallingModifier);


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

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

        end

        function s = objToStructFilteredRecursive(~,obj)
            % Convert the object to a structure using 'struct' and catch warnings
            warnState = warning('off', 'MATLAB:structOnObject');
            s = struct(obj);  % Convert to struct
            warning(warnState);  % Restore previous warning state

            % Get all field names of the struct
            fields = fieldnames(s);

            % Loop over each field and handle filtering
            for i = 1:numel(fields)
                fieldData = s.(fields{i});

                if isstruct(fieldData)  % If the field is a struct, call recursively
                    s.(fields{i}) = obj.objToStructFilteredRecursive(fieldData);
                elseif numel(fieldData) > 1000  % Remove field if it has more than 1000 elements
                    %s = rmfield(s, fields{i});
                    s.(fields{i}) = [];
                elseif isa(fieldData,'table')
                    s = rmfield(s, fields{i});
                end
            end
        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

            if options.fs_in == options.fs_out

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

                obj = obj.logbookentry(desc,obj);

            else

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

                obj.fs = options.fs_out;

            end

        end

        %%
        function spectrum(obj,options)

            arguments
                obj
                options.fignum
                options.displayname = "";
                options.color = [];
                options.normalizeToNyquist = 0;
                options.normalizeTo0dB = 0;
                options.max_num_lines = []; % Leave empty or omit to disable line rotation
                options.fft_length = [];
            end

            if isempty(options.fft_length)
                options.fft_length = 2^(nextpow2(length(obj.signal))-9);
            end
            
            if options.normalizeToNyquist == 0
                [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,obj.fs,"centered","power","mean");
                w = w.*1e-9;
            else
                [p_lin,w] = pwelch(obj.signal,hanning(options.fft_length),options.fft_length/2,options.fft_length,"centered","power","mean");
            end

            if options.normalizeTo0dB
                p_lin = p_lin./ max(p_lin);
                p_dbm = 10*log10(p_lin); % normalized to 0 dB
                ylab = "normalized to 0 dB";
            else
                p_dbm = 10*log10(p_lin);
                ylab = "Power (dB/Hz)";
            end

            figure(options.fignum);
            ax = gca;
            hold on

            if isempty(options.color)
                hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1);
            else
                hLine = plot(w,p_dbm,'DisplayName',options.displayname,'LineWidth',1,'Color',options.color);
            end

            % If user wants to limit the number of lines, check and remove old lines
            if ~isempty(options.max_num_lines) && options.max_num_lines > 0
                allLines = findall(ax, 'Type', 'Line');
                if length(allLines) > options.max_num_lines
                    % Sort lines by creation order. Usually, the oldest lines appear first in allLines.
                    % If needed, you can sort by UserData or other criteria.
                    numToRemove = length(allLines) - options.max_num_lines;
                    delete(allLines(1:numToRemove));
                end
            end

            if options.normalizeToNyquist == 0
                xlabel("Frequency in GHz");
                edgetick = 2^(nextpow2(obj.fs*1e-9));
                xticks(-edgetick:16:edgetick);
                xlim([100*round(min(w)/100,1)-10, 100*round(max(w)/100,1)+10])
                xlim([-128 128]);%256GSa/s
            else
                xlabel("Normalized Frequency");
                xlim([-pi, pi]);
            end

            ylabel(ylab);

            try
            ylim([max(min(floor(min(p_dbm))-3, ax.YLim(1)),-40), min(max(ceil(max(p_dbm))+3, ax.YLim(2)),10)]);
            catch
            ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]);
            end
            ylim([floor(min(p_dbm))-3, ceil(max(p_dbm))+3]);
            yticks(-200:10:10);
            grid on; grid minor;
            legend
            % legend('Interpreter','none');

        end


        function move_it_spectrum(obj,options)
            arguments
                obj
                options.fignum
                options.displayname = "";
                options.color = [];
                options.normalizeToNyquist = 0;
                options.normalizeTo0dB = 0;
            end

            data_in = obj.signal;

            if size(data_in,1) > size(data_in,2)
                data_in = data_in';
            end

            for pol = 1:size(data_in,1)

                %compute FFT of input
                Data_in = fft( data_in(pol,:) );

                %psd = Data_in.*conj(Data_in);
                psd = Data_in;

                %Use only magnitude of FFT (which was complex)
                psd = abs(psd);

                %Shift the spectrum to yield
                psd = fftshift(psd);

                %divide by N
                psd = psd/length(data_in(pol,:));

                psd_plot = 20*log10(psd);

                % psd_plot = psd_plot - max(psd_plot);


                %smoothing
                % psd_smoothed = smooth(psd,1000);
                %
                % psd_smoothed = 10*log10(psd_smoothed);

                % psd_smoothed = psd_smoothed - max(psd_smoothed);

                carrier_power_time_dbm = 20*log10(   mean(abs(data_in)) .^2  )+30; % dB -> +30 -> dBm
                carrier_power_freq_dbm = max(psd_plot);


                % psd_plot = psd_plot - max(psd_plot);

                %% cspr
                c = mean(data_in).^2;
                s = mean(data_in.^2);
                cspr = 10*log10(c / s);

                testParseval = 1;

                if testParseval == 1

                    E_FreqDomain =1/length(psd) * sum((psd.*length(psd)).^2);

                    %test parseval
                    E_TimeDomain = sum( (data_in(pol,:).^2) );

                    if isequal(round(E_FreqDomain,1),round(E_TimeDomain,1))
                        disp('Parseval is right!');
                    else
                        % disp('Something is wrong here?!');
                    end
                end

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

                if 1
                    %Frequency Axis
                    freq_vec = linspace(-obj.fs/2,obj.fs/2,length(psd));
                    freq_vec = reshape(freq_vec,size(psd_plot));


                    if nargin == 4
                        p = plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5,'DisplayName',options.displayname);
                        % plot(freq_vec*1e-9,psd_smoothed','Linewidth',1,'Color',[0 0 0],'DisplayName',[char(varargin{2}),' smoothed']);

                    else
                        p = plot(freq_vec*1e-9,psd_plot,'Linewidth',0.5);
                        % plot(freq_vec*1e-9,psd_smoothed','Linewidth',1,'Color',[1 1 1],'LineStyle',':');

                    end

                    %xlim([freq_vec(1)/1e9-2 freq_vec(end)/1e9+2])
                    xlabel('frequency [GHz]')
                else
                    %Wavelength Axis
                    freq_vec = physconst('LightSpeed')*linspace(-obj.fs/2,obj.fs/2,length(psd))./((physconst('LightSpeed')/1310e-9)^2)*1e9;
                    freq_vec = freq_vec+1310;
                    if nargin == 4
                        plot(freq_vec',psd_plot,'Linewidth',0.5,'DisplayName',options.displayname)
                    else
                        plot(freq_vec,psd_plot','Linewidth',0.5);
                    end

                    %xlim([freq_vec(1)/1e9-2 freq_vec(end)/1e9+2])
                    xlabel('wavelength [nm]')
                end
                hold on


            end

            % xlim([-150 150])
            % ylim([-100,0]);
            ylabel('magnitude [dBm]')
            legend
            grid minor;


        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

        %% Normalize
        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

        %% Delay
        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);
                % obj.signal=circshift(obj.signal,delay);

            elseif options.mode == delay_mode.time

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

            end

        end

        %%
        function [obj,S,inverted,sequenceFound,sequenceStarts] = tsynch(obj,options)
            % time sync and cut
            arguments
                obj Signal
                options.reference Signal
                options.fs_ref = 0;
                options.debug_plots = 0;
            end
     
            
            
            S = {};
            inverted = -1;
            sequenceFound = 0;
            sequenceStarts = [];


            %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;

            max_occurences = floor(length(a)/length(b));

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

            %estimate start pos of signal
            maxpeaknum = floor(length(a)/length(b));
            try
                [pks,pkpos,w,p] = findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend');
            catch
                warning(['Error in findpeaks, ususally the seuqnece is too short. Max peak num: ', num2str(maxpeaknum)]);
                return
            end

            if mean(w) > 10 || mean(p) > 10
                return
            else
                sequenceFound = 1;
            end

            if options.debug_plots
                figure(121212);clf
                subplot(1,2,1);
                findpeaks(abs(co./max(co)),'MinPeakDistance',length(b)/2,'MinPeakHeight',0.2,'NPeaks',maxpeaknum,'SortStr','descend')
            end

            

            shifts = lags(pkpos);
            sequenceStarts = shifts;
            % shifts = shifts(shifts>=0);

            if numel(shifts) > 0
    
                %Cut occurences of ref signal from signal (only positive shifts)
                if all(sign(co(pkpos))==-1)
                    inverted = 1;
                end
               
                for c = shifts
                    sig = obj.delay(-c,'mode','samples');
                    sig.signal = sig.signal(1:length(b)) .* -inverted;
                    S{end+1,1} = sig;
                end
    
                %return/keep the sinal with the highest correlation (only within positive shifts)
                [~,idx]=max(pks);
                obj.signal = S{idx}.signal;
                %put signal with highest corr. to first index in S array
                swap = S{1};
                S{1} = S{idx};
                S{idx} = swap;
    
                for c = 1:numel(shifts)
                    S{c}.logbook = [];
                end

            else

                %do nothing when shifts are negative or there are none...

            end

            %plot all synced signals and the ref signal
            if options.debug_plots
                figure(121212);hold on;
                subplot(1,2,2);
                for i = 1:size(S,1)
                    hold on
                    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
                drawnow;
            end




        end






        %%
        function obj = filter(obj,a,b)

            lbdesc = ['Filtering signal with H = a: ',num2str(a),' / b: ',num2str(b)];
            obj  = obj.logbookentry(lbdesc,obj);

            obj.signal = filter(a,b,obj.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,options)

            arguments
                obj
                fsym
                M
                options.fignum = 100;
                options.displayname = "";
            end

            mode = 1;

            histpoints = 2048;                                              %% verticale resolution
            histpoints = floor(histpoints/2)*2+1;                           %% to have the eye digram centered around one point make the vertical resolution uneven
            histpoints_horizontal = 2048;                                    %% 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

            startpos = floor(0.1*length(sig));
            endpos = floor(0.9*length(sig));
            endpos = min(endpos,startpos+200000);
            x = sig(startpos:endpos);                                               %% 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(options.fignum)
            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))*2;
                minA = min(sig(100:end-100))*2;

                % 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

                ax = gca;

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

                imagesc(plot_data);

                % beautify
                cm=flip(cbrewer2("RdYlBu",4096));
                % cm=flip(cbrewer2("RdBu",4096));
                cm(1,:) = [1,1,1]; % set zeros to white => clean background
                colormap(cm);
                % colormap('turbo');
                % ax.CLim =  [0   50];
                if isa(obj,'Opticalsignal')
                    title(['Optical Eye ',options.displayname])
                    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 ',options.displayname])
                    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 ',options.displayname])
                    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)



                    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",10,"NPeaks",M,"MinPeakHeight",30,"MinPeakProminence",10);

                        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

                        plot_infos = 0;
                        if plot_infos

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



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

                    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