classdef AWG
    %AWG Summary of this class goes here
    %   Detailed explanation goes here

    properties(Access=public)

        kover %oversampling factor e.g. 16
        repetitions %repeat the signal to generate a longer sequence?
        fdac %needed
        normalize %want to normalize at first? either 0 or 1
        bit_resolution %bit res. of quantizer (e.g. 5 bit)
        dac_min
        dac_max

        skew_active = 0;
        awg_skew ; %skew vector for each output channel

        lpf_active = 0;
        lpf_type ;
        f_cutoff ;
        lowpass ;

        signal_length;
        H_lpf;


    end

    methods (Access=public)
        function obj = AWG(options)
            %AWG Construct an instance of this class
            %   Detailed explanation goes here

            arguments
                options.preset          = []      ;    
                options.kover           = 16;
                options.repetitions     = 1;
                options.normalize       = 1;
                options.fdac            = 92e9;
                options.bit_resolution  = 5.5
                options.dac_min         = -0.5;
                options.dac_max         = 0.5;
                options.skew_active     = 0;
                options.awg_skew        = 0; 
                options.lpf_active      = 0;
                options.lpf_type        = 0;
                options.f_cutoff        = 32e9;
                options.lowpass         = 1;
            end

            if isempty(options.preset)
                obj.skew_active = 0;
                obj.lpf_active = 0;

            elseif options.preset == "M8196A"
                % M8196A (92GBd) https://www.keysight.com/us/en/product/M8196A/92-gsa-s-arbitrary-waveform-generators.html
                obj.dac_max         = 0.5;
                obj.dac_min         = -.5;
                obj.lowpass         = 1; %LP
                obj.f_cutoff        = 32e9;
            elseif options.preset == "M8199B"
                %https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf
            end

            obj.kover = options.kover; %oversampling factor e.g. 16
            obj.repetitions = options.repetitions; %repeat the signal to generate a longer sequence?
            obj.fdac = options.fdac;
            obj.normalize = options.normalize;%want to normalize at first? either 0 or 1
            obj.bit_resolution = options.bit_resolution;%bit res. of quantizer (e.g. 5 bit)
            obj.dac_min = options.dac_min;
            obj.dac_max = options.dac_max;

            obj.lpf_active = options.lpf_active;

            obj.skew_active = options.skew_active;
            obj.awg_skew = options.awg_skew;

            

        end

        function elec_out = process_channel(obj,data_in)
            %METHOD1 Summary of this method goes here
            %   Detailed explanation goes here
            arguments(Input)
                obj
                data_in 
            end
            
            arguments(Output)
                elec_out 
            end

            obj.signal_length = length(data_in);

            if obj.normalize
                % 5.1. Normalize the signal to 1 Vpp and set the amplitude of the signal
                data_in = data_in/(max(data_in)-min(data_in));
            else
                % 5.1. Cut the Signal at -1 and 1 and scale to amplitude
                data_in(data_in > 1) = 1;
                data_in(data_in < -1) = -1;
            end


            % 1. Quantize the signal
            if obj.bit_resolution>0
                elec_out = obj.quantization(data_in) ;
            else
                elec_out = data_in;
            end

            % 2. Sample and hold + repeat (data_out: 1xsignal length)
            elec_out = repmat(elec_out,obj.repetitions,obj.kover);
            elec_out = reshape(elec_out',[],1);

            % 3. Add skew
            if obj.skew_active
                elec_out = obj.skew(elec_out);
            end         

            % 4. Apply LPF on the signal
            if obj.lpf_active
                lpf = Filter('filtdegree',4,"f_cutoff",obj.f_cutoff,"fsamp",obj.kover*obj.fdac,"filterType",filtertypes.bessel_inp);
                elec_out = lpf.process(elec_out);
%                 obj.H_lpf = obj.buildFilter(1);
%                 data_out = obj.lpf(data_out) ;
            end

        end

    end

    methods (Access=private)

        function quant_out = quantization(obj,x_in)

            steps = round(2^obj.bit_resolution) + rem(round(2^obj.bit_resolution),2) ;

            if isreal(x_in)
                % shift signal and clip to quantizer intervall
                x_in =   min(max(x_in-obj.dac_min,0),obj.dac_max-obj.dac_min);

                % quantize signal
                x_in = round((steps-1)/(obj.dac_max-obj.dac_min)*x_in);

                % scale and shift back
                quant_out = x_in*(obj.dac_max-obj.dac_min)/(steps-1) +obj.dac_min;

            else

                % shift signal and clip to quantizer intervall
                x_in = min(max(real(x_in)-obj.dac_min,0),obj.dac_max-obj.dac_min) + ...
                    1i*min(max(imag(x_in)-obj.dac_min,0),obj.dac_max-obj.dac_min);

                % quantize signal and shift back signal
                x_in = round((steps-1)/(obj.max-obj.dac_min)*x_in);

                % scale and shift back
                quant_out = x_in*(obj.dac_max-obj.dac_min)/(steps-1) +(1+1i)*obj.dac_min;

            end

        end

        function skew_out = skew(obj,x_in)

            % Exact delay calculation

            % Calculate transfer function
            fsimu = obj.kover*obj.fdac;
            faxis = linspace(fsimu/2, fsimu/2, length(x_in)+1);
            faxis = fftshift(faxis(1:end-1));

            skew_out = ifft(fft(x_in) .* exp(-1i*2*pi*obj.awg_skew*faxis)' );

            skew_out = real(skew_out) ; % get rid of negligible imaginary part

        end

        function lpf_out = lpf(obj,x_in)
            lpf_out = ifft(obj.H_lpf.*fft(x_in));
        end

        function H = buildFilter(obj,filterType)
            
            filtdegree = 3;
            fsimu = obj.kover*obj.fdac;

            rp = 0.5; %passband ripple
            rs = 0.5; %stopband ripple

            switch filterType
                case 1

                    % Bessel filter, impulse invariant transformed

                    [B, A] = besself(filtdegree, 2*pi*obj.f_cutoff);
                    [B ,A] = impinvar(B,A,fsimu);

                case 2

                    % Bessel filter, impulse bilinear transformed

                    [Z, P, K] = besself(filtdegree, 2*pi*obj.f_cutoff);
                    [Z ,P, K] = bilinear(Z,P,K,fsimu);
                    [B ,A] = zp2tf(Z ,P ,K);

                case 3

                    % Butterworth filter
                    if obj.lowpass == 1 %lowpass
                        [B, A] = butter(filtdegree, obj.f_cutoff/(fsimu/2),'low');
                    else % highpass
                        [B, A] = butter(filtdegree, obj.f_cutoff/(fsimu/2),'high');
                    end

                case 4

                    % Chebyshev 1 filter

                    [B, A] = cheby1(filtdegree,rp, obj.f_cutoff/(fsimu/2));

                case 5

                    % Chebyshev 2 filter

                    [B, A] = cheby2(filtdegree,rs, obj.f_cutoff/(fsimu/2));

                case 6

                    % Elliptic filter

                    [B, A] = ellip(filtdegree,rp,rs,obj.f_cutoff/(fsimu/2));

                case 7

                    % Hamming filter

                    g=(filtdegree-1)/2;
                    wc=obj.f_cutoff/(fsimu/2);
                    B = wc*sinc(wc*(-g:g)).*hamming(filtdegree)';
                    A=1;

                case 8

                    % Raised Cosine filter

                    B = firrcos(filtdegree,obj.f_cutoff,para.df,fsimu);
                    A=1;

                case 9

                    % Sinc filter

                    g=(filtdegree-1)/2;
                    wc=obj.f_cutoff/(fsimu/2);
                    B = wc*sinc(wc*(-g:g));
                    A=1;

            end

             H = freqz(B, A, obj.repetitions*obj.kover*obj.signal_length,'whole');


        end

    end
end