classdef OSA_Advantest < handle
    % Advantest Q8384 Optical Spectrum Analyzer (GPIB)
    % User inputs in nanometers

    % EXAMPLE:

    % span_nm     = 70;
    % lambda_c_nm = 1310;
    % osa = OSA_Advantest( ...
    %     'span_nm',        span_nm, ...
    %     'lambda_c_nm',    lambda_c_nm, ...
    %     'resolution_nm',  0.1, ...
    %     'sampling_points', 10001 );
    %
    % osa.configure();
    % [lambda_nm, psd_dBm] = osa.measure();
    % figure()
    % plot(lambda_nm,psd_dBm);

    properties (Access = public)
        span_nm double = 5
        resolution_nm double = 0.05
        lambda_c_nm double = 1550
        sampling_points double = 1001
        timeout_s double = 10

        lambda_nm
        psd_dBm
    end

    properties (Constant)
        valid_resolutions_nm = [0.01 0.02 0.05 0.10 0.20 0.50]
        valid_sampling_points = [101 201 501 1001 2001 5001 10001];
    end

    properties (Access = protected)
        visaAddress string = "GPIB0::3::INSTR"
    end

    methods
        function obj = OSA_Advantest(options)
            arguments
                options.visaAddress string = "GPIB0::3::INSTR"
                options.span_nm double = 5
                options.resolution_nm double = 0.05
                options.lambda_c_nm double = 1550
                options.sampling_points double = 1001
                options.timeout_s double = 10
            end

            fn = fieldnames(options);
            for k = 1:numel(fn)
                obj.(fn{k}) = options.(fn{k});
            end

            assert(ismember(obj.resolution_nm, obj.valid_resolutions_nm), ...
                'Invalid resolution. Allowed: 0.01 | 0.02 | 0.05 | 0.10 | 0.20 | 0.50 nm');
            assert(ismember(obj.sampling_points, obj.valid_sampling_points), ...
                'Invalid # of Samp. Points. Allowed: 101 201 501 1001 2001 5001 10001');

            v = obj.connect();
            writeline(v,"*IDN?");
            disp(readline(v));
            delete(v);
        end

        function configure(obj)
            assert(ismember(obj.resolution_nm, obj.valid_resolutions_nm), ...
                'Invalid resolution. Allowed: 0.01 | 0.02 | 0.05 | 0.10 | 0.20 | 0.50 nm');
            assert(ismember(obj.sampling_points, obj.valid_sampling_points), ...
                'Invalid # of Samp. Points. Allowed: 101 201 501 1001 2001 5001 10001');

            v = obj.connect();

            writeline(v,"FMT0,HED0,SDL0");
            writeline(v,"BUZ0");
            writeline(v,"QUI0");

            writeline(v, sprintf("RES %gnm", obj.resolution_nm));
            writeline(v, sprintf("SPT %g",  obj.sampling_points));
            writeline(v, sprintf("CEN %gnm", obj.lambda_c_nm));
            writeline(v, sprintf("SPA %gnm", obj.span_nm));
            writeline(v,"SWE 0");

            delete(v);
        end

        function [osa_x_nm, osa_y_dBm] = measure(obj)
            v = obj.connect();

            writeline(v,"CSB");
            writeline(v,"MEA1");

            active = 1; cnt = 0;
            while active
                pause(1);
                writeline(v,"MEA?");
                active = str2double(readline(v));
                cnt = cnt + 1;
                if cnt > obj.timeout_s
                    delete(v);
                    error("OSA_Q8384:Timeout","Measurement timeout");
                end
            end

            writeline(v,"PKL");

            writeline(v,"ODN?");
            readline(v); % number of points (not strictly needed)

            writeline(v,"FMT0,HED0,SDL0");

            writeline(v,"OSD0");           % Y-axis
            osa_y_dBm = str2num(readline(v));

            writeline(v,"OSD1");           % X-axis
            osa_x_m = str2num(readline(v));
            osa_x_nm =  osa_x_m.*1e9;
            delete(v);

            obj.lambda_nm=osa_x_nm;
            obj.psd_dBm=osa_y_dBm;

        end

        % --- Replace / extend plot() in OSA_Q8384 class ---

        function plot(obj, options)
            arguments
                obj
                options.FigureNumber double = 1
                options.Hold logical = true
                options.Title string = "OSA Spectrum"
                options.LineWidth double = 0.8
                options.Grid logical = true
                options.Color (1,3) double = [0 0.4470 0.7410]
                options.DisplayName string = ""
            end

            if isempty(obj.psd_dBm) || isempty(obj.lambda_nm)
                error('Record spectrum first: run osa.measure()');
            end

            fig = figure(options.FigureNumber);
            fig.Color = 'w';

            ax = gca;
            if ~options.Hold
                cla(ax);
            else
                hold(ax,'on');
            end

            plot(ax, obj.lambda_nm, obj.psd_dBm, ...
                'LineWidth', options.LineWidth, ...
                'Color', options.Color, ...
                'DisplayName', options.DisplayName);

            xlabel(ax,'Wavelength [nm]');
            ylabel(ax,'Optical Power [dBm]');

            if options.Title ~= ""
                title(ax, options.Title);
            end

            xlim(ax,[min(obj.lambda_nm) max(obj.lambda_nm)]);
            yl = ylim(ax);
            ylim(ax,[yl(1)-1 yl(2)+1]);

            if options.Grid
                grid(ax,'on');
                grid(ax,'minor');
            end

            set(ax, ...
                'FontSize',11, ...
                'LineWidth',1, ...
                'Box','on');

            if options.DisplayName ~= ""
                legend(ax,'show','Location','best');
            end
        end

        function osnr = computeOSNR(obj, options)
            % OSNR computation based on ratio method (RBW-independent)
            % Assumes obj.lambda_nm [nm], obj.psd_dBm [dBm/RBW] exist

            arguments
                obj
                options.bw_signal_nm double = 0.1    % channel bandwidth (0.1 for CW?!?!)
                options.ref_bw_nm double = 0.1       % OSNR reference bandwidth
                options.carrier_frequency = NaN;
                options.noise_guard_nm = 1;
                options.noise_window_nm =1;
                options.debugplots = 0;
            end

            assert(~isempty(obj.lambda_nm) && ~isempty(obj.psd_dBm), ...
                'Run measure() before computing OSNR');

            obj.lambda_nm = obj.lambda_nm;
            psd    = obj.psd_dBm;

            % --- find carrier wavelength recorded in OSA ---
            % options.carrier_frequency [nm], options.noise_guard_nm, options.noise_window_nm

            if ~isnan(options.carrier_frequency)

                % local search window around expected CW carrier
                idx_search = obj.lambda_nm > (options.carrier_frequency - options.noise_guard_nm) & ...
                    obj.lambda_nm < (options.carrier_frequency + options.noise_guard_nm);

                [~, imax_local] = max(psd(idx_search));
                lambda_search = obj.lambda_nm(idx_search);
                lambda_c = lambda_search(imax_local);
                [~,imax] = find(obj.lambda_nm==lambda_c);

                % define local noise windows directly adjacent to signal (guarded)
                idx_noise_left  = obj.lambda_nm > (lambda_c - options.noise_guard_nm - options.noise_window_nm) & ...
                    obj.lambda_nm < (lambda_c - options.noise_guard_nm);

                idx_noise_right = obj.lambda_nm > (lambda_c + options.noise_guard_nm) & ...
                    obj.lambda_nm < (lambda_c + options.noise_guard_nm + options.noise_window_nm);

                idx_noise = idx_noise_left | idx_noise_right;

            else
                % fallback: global maximum (unsafe for pathological ASE cases)
                [~, imax] = max(psd);
                lambda_c = obj.lambda_nm(imax);

                % ASE from outer quarters of spectrum
                N = numel(obj.lambda_nm);
                idx_noise = false(size(obj.lambda_nm));
                idx_noise([3:round(N/4), round(3*N/4):N-1]) = true;
            end

            % --- extract signal window ---
            idx_signal = obj.lambda_nm > (lambda_c - options.bw_signal_nm/2) & ...
                obj.lambda_nm < (lambda_c + options.bw_signal_nm/2);

            if ~any(idx_signal)
                idx_signal(imax) = true;
            end

            cc_lambda = obj.lambda_nm(idx_signal);
            cc_psd    = psd(idx_signal);

            % --- estimate ASE from local noise (flat ASE assumption) ---
            noise_lambda = obj.lambda_nm(idx_noise);
            noise_psd    = psd(idx_noise);


            ase_psd = interp1(noise_lambda, noise_psd, cc_lambda, ...
                'linear', 'extrap');

            % --- convert to linear (mW / RBW) ---
            sig_lin = 10.^(cc_psd/10);
            ase_lin = 10.^(ase_psd/10);

            % --- integrate (RBW cancels) ---
            P_sig_lin = sum(sig_lin);
            P_ase_lin = sum(ase_lin);

            P_sig_dB = 10*log10(P_sig_lin);
            P_ase_dB = 10*log10(P_ase_lin);

            osnr.p_sig_db = P_sig_dB;
            osnr.p_ase_db = P_ase_dB;

            % --- OSNR direct from peak value to the noise estimate out of band ---
            osnr.direct_dB = psd(imax) - 10*log10(mean(ase_lin)) ;

            % --- OSNR in measurement bandwidth ---
            osnr.measured_dB = P_sig_dB - P_ase_dB;

            % --- normalize to reference bandwidth (usually 0.1 nm) ---
            osnr.ref_dB = osnr.measured_dB + ...
                10*log10(options.bw_signal_nm / options.ref_bw_nm);

            % --- ASE-corrected OSNR (optional definition) ---
            osnr.corrected_dB = ...
                10*log10(10^(P_sig_dB/10) - 10^(P_ase_dB/10)) - P_ase_dB;

            osnr.corrected_ref_dB = osnr.corrected_dB + ...
                10*log10(options.bw_signal_nm / options.ref_bw_nm);

            % metadata
            osnr.lambda_c_nm = lambda_c;
            osnr.bw_signal_nm = options.bw_signal_nm;
            osnr.ref_bw_nm = options.ref_bw_nm;

            if options.debugplots
                figure(100)
                xline(lambda_c,'LineWidth',2,'DisplayName','Max. OSA Value','Color',[0.8,0.8,0.8])
                obj.plot("DisplayName",'OSNR','FigureNumber',100);
                yline(osnr.p_ase_db,'LineWidth',2,'DisplayName','Signal Estimation');
                yline(osnr.p_sig_db,'LineWidth',2,'DisplayName','Noise Estimation');
                scatter(obj.lambda_nm(idx_noise),obj.psd_dBm(idx_noise),'Marker','.','LineWidth',4,'DisplayName','Noise Samples','MarkerEdgeColor','red');
                scatter(obj.lambda_nm(idx_signal),obj.psd_dBm(idx_signal),'Marker','.','LineWidth',4,'DisplayName','Signal Samples','MarkerEdgeColor','green');
                scatter(obj.lambda_nm(idx_signal),ase_psd,'Marker','.','LineWidth',4,'DisplayName','Signal Samples','MarkerEdgeColor','magenta');
                ylim([min(obj.psd_dBm), max([obj.psd_dBm,osnr.p_sig_db+5])]);

                figure(101);hold on
                xline(lambda_c,'LineWidth',2,'DisplayName','Max. OSA Value','Color',[0.8,0.8,0.8])
                xline(lambda_c - options.bw_signal_nm/2,'LineWidth',2,'DisplayName','CW - 0.5 BW','Color',[0.6,0.6,0.8])
                xline(lambda_c + options.bw_signal_nm/2,'LineWidth',2,'DisplayName','CW + 0.5 BW','Color',[0.6,0.8,0.8])
                obj.plot("DisplayName",'OSNR','FigureNumber',101);
                yline(osnr.p_ase_db,'LineWidth',2,'DisplayName','Signal Estimation');
                yline(osnr.p_sig_db,'LineWidth',2,'DisplayName','Noise Estimation');
                scatter(obj.lambda_nm(idx_signal),obj.psd_dBm(idx_signal),'Marker','.','LineWidth',4,'DisplayName','Signal Samples','MarkerEdgeColor','green');
                scatter(obj.lambda_nm(idx_signal),ase_psd,'Marker','.','LineWidth',4,'DisplayName','Signal Samples','MarkerEdgeColor','magenta');
                xlim([lambda_c - options.bw_signal_nm, lambda_c + options.bw_signal_nm])
                ylim([min(obj.psd_dBm), max([obj.psd_dBm,osnr.p_sig_db+5])]);
            end

        end

    end

    methods (Access = protected)
        function v = connect(obj)
            v = visadev(obj.visaAddress);
        end
    end
end