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
                options.polrot
            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,"polrot",options.polrot);


            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 = randi(1000)
                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(X,'Signal') && isa(y,'Signal')
                Sum = X;
                Sum.signal = X.signal + y.signal;
            elseif isa(X,'Signal') && isnumeric(y)
                Sum = X;
                Sum.signal = X.signal + y;
            elseif isnumeric(X) && isa(y,'Signal')
                Sum = y;
                Sum.signal = X + y.signal;
            end

        end

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

            if isa(X,'Signal') && isa(y,'Signal')
                Diff = X;
                Diff.signal = X.signal - y.signal;
            elseif isa(X,'Signal') && isnumeric(y)
                Diff = X;
                Diff.signal = X.signal - y;
            elseif isnumeric(X) && isa(y,'Signal')
                Diff = y;
                Diff.signal = X - y.signal;
            end

        end

        function Product = times(X,y)

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

        end

        function Product = mtimes(X,y)

            if (isa(X,'Signal') && isnumeric(y) && isscalar(y)) || ...
                    (isnumeric(X) && isscalar(X) && isa(y,'Signal'))
                Product = times(X,y);
            else
                error('Signal:mtimes:UnsupportedOperands', ...
                    'Use element-wise .* for Signal multiplication, or scalar * Signal for scaling.');
            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;
            SignalCopy  = [];
            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

                [p, q] = rat(options.fs_out / options.fs_in);
                obj.signal  = resample(obj.signal,p,q,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 = 2025
                options.displayname = "";
                options.color = [];
                options.linestyle = '-';
                options.HandleVisibility (1,1) string {mustBeMember(options.HandleVisibility, ["on","off"])} = "on";
                options.normalizeToNyquist = 0;
                options.normalizeToSamplingRate = 0;
                options.addDCoffset = 0;
                options.normalizeToDC = 0;
                options.normalizeTo0dB = 0;
                options.show_onesided = false;
                options.max_num_lines = []; % Leave empty or omit to disable line rotation
                options.fft_length = [];
                % --- NEW options ---
                options.useWavelengthAxis (1,1) logical = false  % plot x-axis in wavelength
                options.lambda0_nm (1,1) double = 1310           % center wavelength [nm]
            end

            if isempty(options.fft_length)
                options.fft_length = 2^(nextpow2(length(obj.signal))-9);
            end

           
            useSamplingRateAxis = options.normalizeToSamplingRate ~= 0;
            useRadPerSampleAxis = options.normalizeToNyquist ~= 0 && ~useSamplingRateAxis;

            if ~useRadPerSampleAxis && ~useSamplingRateAxis
                [p_lin,f_Hz] = pwelch(obj.signal, hanning(options.fft_length), ...
                    options.fft_length/2, options.fft_length, ...
                    obj.fs, "centered", "power", "mean");
                f_GHz = f_Hz*1e-9; % keep frequency vector for frequency axis
            else
                [p_lin,f_rad] = pwelch(obj.signal, hanning(options.fft_length), ...
                    options.fft_length/2, options.fft_length, ...
                    "centered", "power", "mean");
                % In normalized modes, pwelch returns rad/sample centered on 0.
                % Divide by 2*pi for the f/fs axis where Nyquist is 0.5.
            end

            p_lin = movmean(p_lin,10);

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

            % --- If requested, build wavelength axis from frequency offset ---
            if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis
                c = physconst('LightSpeed');         % [m/s]
                lambda0_m = options.lambda0_nm*1e-9; % center wavelength [m]
                f_c = c / lambda0_m;                 % carrier frequency [Hz]

                % exact mapping
                f_abs    = f_c + f_Hz;               % absolute frequency [Hz]
                lambda_m = c ./ f_abs;               % wavelength [m]
                lambda_nm = lambda_m * 1e9;          % wavelength [nm]

                % assign axis
                x_vec   = lambda_nm(:);
                x_label = "Wavelength [nm]";
                dc_axis = f_Hz;

                % Sort to ensure axis is ascending
                [x_vec, sortIdx] = sort(x_vec, 'ascend');
                p_dbm = p_dbm(sortIdx, :);
                dc_axis = dc_axis(sortIdx);
            else
                % Frequency or normalized axes
                if ~useRadPerSampleAxis && ~useSamplingRateAxis
                    x_vec   = f_GHz;
                    x_label = "Frequency in GHz";
                    dc_axis = f_GHz;
                elseif useSamplingRateAxis
                    x_vec   = f_rad ./ (2*pi);
                    x_label = "Normalized Frequency f/fs";
                    dc_axis = x_vec;
                else
                    x_vec   = f_rad;     % normalized frequency in rad/sample
                    x_label = "Normalized Frequency [rad/sample]";
                    dc_axis = x_vec;
                end
            end

            if options.show_onesided
                keep_idx = dc_axis >= 0;
                x_vec = x_vec(keep_idx);
                dc_axis = dc_axis(keep_idx);
                p_dbm = p_dbm(keep_idx, :);
            end

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

            p_dbm = p_dbm+options.addDCoffset;

            if options.normalizeToDC
                [~,min_idx]=min(abs(dc_axis));
                pow_at_dc = p_dbm(min_idx);
                p_dbm = p_dbm-pow_at_dc;
            end
            % p_dbm = movmean(p_dbm,10);

            for s = 1:min(size(p_dbm))
                if isempty(options.color)
                    plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1, 'LineStyle', options.linestyle, 'HandleVisibility', options.HandleVisibility);
                else
                    plot(x_vec, p_dbm(:,s), 'DisplayName', options.displayname, 'LineWidth', 1, 'Color', options.color,'LineStyle',options.linestyle, 'HandleVisibility', options.HandleVisibility);
                end
            end

            % Limit number of lines if requested
            if ~isempty(options.max_num_lines) && options.max_num_lines > 0
                allLines = findall(ax, 'Type', 'Line');
                if length(allLines) > options.max_num_lines
                    numToRemove = length(allLines) - options.max_num_lines;
                    delete(allLines(1:numToRemove));
                end
            end

            % Axis labels and limits
            xlabel(x_label);

            if options.useWavelengthAxis && ~useRadPerSampleAxis && ~useSamplingRateAxis
                xlim([min(x_vec) max(x_vec)]);
            else
                if ~useRadPerSampleAxis && ~useSamplingRateAxis
                    % Keep your existing freq handling (you can fine-tune as needed)
                    % xlim([-128 128]); % example for 256 GSa/s if desired
                    xlim([min(x_vec) max(x_vec)]);
                elseif useSamplingRateAxis
                    if options.show_onesided
                        xlim([0, 0.5]);
                    else
                        xlim([-0.5, 0.5]);
                    end
                else
                    if options.show_onesided
                        xlim([0, pi]);
                    else
                        xlim([-pi, pi]);
                    end
                end
            end

            ylabel(ylab);

            % --- Y-Axis scaling (auto with margin) ---
            y_min = min(p_dbm(:));
            y_max = max(p_dbm(:));

            % Add 5% dynamic range margin on both sides
            y_range = y_max - y_min;
            if y_range == 0
                y_range = 10; % fallback if flat
            end

            y_margin = 0.05 * y_range;
            ylim([y_min - y_margin, y_max + y_margin]);

            % Set ticks automatically, avoid overpopulation
            try
            yticks(round(linspace(y_min, y_max, min(10, max(4, ceil(y_range/10))))));
            end
            grid on;
            
            % Add legend if not already present
            if isempty(get(gca, 'Legend'))
                legend;
            end


        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 = sum(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

            pkpos = sort(pkpos);

            if isempty(pks)
                warning(['Error in findpeaks, ususally the seuqnece is too short. No Peaks detected']);
                return
            end

            if max(p) < 0.3 || median(w) > 15 
                %median(w) > 15 part means “reject if the detected correlation peaks are too broad.” That can be sensible: a true sync peak should often be sharp.
                warning(['Error in findpeaks, ususally the seuqnece is too short. max(p) = ',num2str(max(p)),'; median(w)=',num2str(median(w)),'']);
                return
            end

            sequenceFound = 1;

            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

                if ~isempty(S)
                    obj = S{1};
                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 = "";
                options.mode = 1; %1= histogram method; 2= intuitive "line based" eye
            end

            mode = options.mode;

            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

                xlabel('Samples','Interpreter','latex')
                ylabel('Amplitude of Signal','Interpreter','latex');
                xlim([0 histpoints_horizontal])

            elseif mode == 1
                % generate eye diagram using histogram

                finite_eye = eye_mat(isfinite(eye_mat));
                if isempty(finite_eye)
                    finite_eye = sig(isfinite(sig));
                end
                amp_min = min(finite_eye);
                amp_max = max(finite_eye);
                amp_center = (amp_max + amp_min) / 2;
                amp_span = amp_max - amp_min;
                if amp_span == 0
                    amp_span = max(abs(amp_center),1);
                end
                amp_margin = 0.08 * amp_span;
                maxA = amp_center + amp_span/2 + amp_margin;
                minA = amp_center - amp_span/2 - amp_margin;
                if ~isa(obj,'Opticalsignal') && minA < 0 && maxA > 0
                    targetStep = max(abs([minA maxA])) / 2;
                    if targetStep > 0
                        stepMagnitude = 10^floor(log10(targetStep));
                        normalizedStep = targetStep / stepMagnitude;
                        if normalizedStep <= 1
                            tickStep = stepMagnitude;
                        elseif normalizedStep <= 2
                            tickStep = 2 * stepMagnitude;
                        elseif normalizedStep <= 5
                            tickStep = 5 * stepMagnitude;
                        else
                            tickStep = 10 * stepMagnitude;
                        end
                        axisLimit = 2 * tickStep;
                        maxA = axisLimit;
                        minA = -axisLimit;
                    end
                end

                % maxA = 0.12;
                % minA = -0.08;

                difference= maxA-minA;

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

                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
                cm=flip(cbrewer2("RdYlBu",4096));
                % cm=flip(cbrewer2("RdBu",4096));
                % cm=flip(cbrewer2("Blues",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");
                    yTickValues = linspace(maxA.*1e3,minA.*1e3,5);
                    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");
                    yTickValues = linspace(maxA,minA,5);
                    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");
                    yTickValues = linspace(maxA,minA,5);
                    min_ = min(obj.signal(100:end-100));
                    max_ = abs(max(obj.signal(100:end-100)));
                end
                xlabel('Time in ps')



                % add information

                if 0

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

                        minPeakDistance = max(10, floor(histpoints / (2*M)));
                        minPeakProminence = max(3, 0.05 * max(hist_interest_smoth));
                        [pk,loc] = findpeaks(hist_interest_smoth, ...
                            "MinPeakDistance",minPeakDistance, ...
                            "NPeaks",M, ...
                            "MinPeakProminence",minPeakProminence, ...
                            "SortStr","descend");
                        loc = sort(loc);

                        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

                    grid off

                end


                yTickPositions = linspace(1,histpoints,numel(yTickValues));
                yticks(yTickPositions);
                yticklabels(sprintfc('%.2f', yTickValues));

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

                % 
            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