classdef ChannelFreqResp < handle
    % Linear pre-compensation of Channel effects. This module has two modes:

    % A) In acquire mode it sends a real OFDM singal with all subcarriers assigned in
    % order to acquire later the system's frequency response after the
    % signal passed the system...
    % ChannelFreqResp.SendOFDM

    % B) In the "not acquire" mode it loads the frequency response generated by the freq_res module,
    % calculates the inverse frequency response with the right length matched
    % with the signal length and distorts the signal before tranmission
    % give channel indices for measure mode

    %TYPICAL FLOW:

    % referencesignal = ChannelFreqResp.buildOFDM()

    % referencesignal -> SYSTEM -> measuredsignal

    % ChannelFreqResp.estimate(measuredsignal)

    properties(Access=public)
        Nacq
        Navg
        Ncp

        f_ref %fs of transmitted dmt sig
        f_observe %fs of the observed signal

        symlen
        seqlen

        refsig
        H
        H_all
        H_apply

        Nfft
        df
        faxis
    end

    methods (Access=public)

        function obj = ChannelFreqResp(options)
            %NAME Construct an instance of this class
            %   Detailed explanation goes here

            arguments
                options.Nacq = 1024;
                options.Navg = 64;
                options.Ncp = 63;
                options.f_ref;
                options.f_observe;
            end

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

            obj.Nfft = 2*obj.Nacq + 1 ;
            obj.df = obj.f_ref/obj.Nfft ; % calculate frequency grid
            obj.faxis = (0:obj.Nfft-1)*obj.df ;

        end

        function output = buildOFDM(obj)
            % generate OFDM or DMT signal that is send over the channel

            obj.symlen    = (obj.Nacq*2+obj.Ncp+1);

            obj.seqlen    = obj.symlen*obj.Navg;

            obj.refsig    = obj.genRNDDMT(15)';

            obj.refsig    = Informationsignal(obj.refsig,"fs",obj.f_ref);

            output = obj.refsig;

        end

        function output = estimate(obj, data_in, options)
            arguments
                obj
                data_in
                options.save logical = false
                options.savePath = "";
                options.fileName = "";
            end
            % use the transmitted DMT and the received DMT to get the transfer function

            obj.Nfft = 2*obj.Nacq + 1 ;

            %resample from fadc to fdac (fdac => fs of reference signal)
            data_in = data_in.resample("fs_out",obj.f_ref);

            %est from dmt
            obj.H_all = obj.estHfromDMT(data_in.signal,obj.refsig.signal);

            obj.H     = mean(obj.H_all,1);

            output = obj.H;

            if options.save
                obj.save('fileName',options.fileName,'savePath',options.savePath);
            end

        end

        function Target = precomp(obj,Target,options)
            % apply the freq response
            arguments
                obj
                Target
                options.fileName = ''
                options.loadPath = ''
                options.maxampdb
            end

            if isempty(obj.H)
                obj.load("fileName",options.fileName,"loadPath",options.loadPath);
            end

            H_inv = 1./obj.H;

            nH    = find(~isnan(H_inv),1,'last'); %last value that is not nan
            H_inv(nH:end)=H_inv(nH); %replace everything after first NaN with the last non-NaN value

            fstarget = Target.fs;

            % Build new frequency axis (with current fs)
            fnew = linspace(0,fstarget/2,length(Target)/2+1);
            fnew = fnew(2:end-1);

            % Old frequency axis (should be much coarser)
            idx_old  = find((obj.faxis > 0) .* (obj.faxis < fstarget/2)); %positions of all Frequencies smaller than fs/2
            int_fold = obj.faxis(idx_old); %old frequencies from 0 to fs/2

            idx_new = find(fnew <= int_fold(end)); %positions of all Frequencies smaller than fs/2
            int_fnew = fnew(idx_new);%new frequencies from 0 to fs/2

            % interpolate the frequency response that had a coarse frequency resolution (e.g. 256 bins) to the current frequency resolution (e.g. 21843 bins)
            iH  = interp1(int_fold, real(H_inv(idx_old)) ,fnew, 'linear') + 1i*interp1(int_fold, imag(H_inv(idx_old)) ,fnew, 'linear');

            % set all NaN values to the fist/ last non-NaN value
            nH    = find(~isnan(iH),1,'first');
            iH(1:nH)=iH(nH);
            nH    = find(~isnan(iH),1,'last');
            iH(nH:end) =iH(nH);

            %smoothing takes time and sometimes the result looks odd,
            %however the performance is most of the time better
            smoothing = 0;
            if smoothing
                iH = smooth(fnew,iH,0.1,'loess')';
            end

            % multiply the complex frequency responce with the phase
            % angle (-pi,pi) at lowest frequency
            iH  = iH.*exp(-1j*angle(iH(1))); % to be checked (<- not from silas, so what needs to be checked?)


            % normalize complex freq. resp. by magnitude at the first
            % five frequencies -> should be the vaue at f=0=DC component?
            iH  = iH./mean(abs(iH)); %why 1:5??

            % set maximum amplification
            % set als values higher than hmax to hmax and keep the
            % phase information by multiplication with respective
            % corresponding phase angles
            maxamp_lin = 10^(options.maxampdb/20);

            iH(abs(iH)>maxamp_lin) = maxamp_lin.*exp(1j*angle(iH(abs(iH)>maxamp_lin)));

            % it could be helpful to clip at linear 1 (=to keep comp from attenuating the signal)
            % iH(abs(iH)<1) = 1.*exp(1j*angle(iH(abs(iH)<1)));

            if 0
                figure(7);hold on;plot(fnew,20*log10(abs(iH)))
            end

            % iH(1) is DC ---> iH(end) is High Freq.
            
            if mod(length(Target.signal),2) %ungerade
                H_inv = [iH(1) iH iH(end) fliplr(conj(iH)) conj(iH(1))];
            else
                H_inv = [iH(1) iH 0 fliplr(conj(iH))];
            end
            % H_inv = [iH(1) iH iH(end) fliplr(conj(iH)) conj(iH(1))];

            obj.H_apply = H_inv;

            Target.signal = real((ifft( ( fft(real( Target.signal' )) .* H_inv ) )));
            Target.signal = Target.signal';

        end

        function plot(obj)

            figure(55);
            %clf;

            Havg = obj.H;

            %1)
            subplot(2,1,1);hold on;box on;title('Magnitude Freq. Response');
            plot(obj.faxis/1e9, 20*log10(abs(obj.H_all)),'linewidth',0.1,'LineStyle','-','Color','#808080') ;
            xlim([0.2 .5*max(obj.faxis)*1e-9]);
            plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2);
            grid on;

            %2)
            subplot(2,1,2); hold on; box on; title('Phase Freq. Response');
            plot(obj.faxis/1e9, angle(obj.H_all),'linewidth',0.1,'LineStyle','-','Color','#808080') ;
            plot(obj.faxis/1e9, unwrap(angle(Havg)),'LineWidth',2) ;
            xlim([0.2 .5*max(obj.faxis)*1e-9]);
            grid on;

            figure(56);
            clf;


            %normalize / remove attenuation
            Havg    = Havg./mean(Havg(2:10));

            %3)
            subplot(2,1,1); hold on; box on; title('Inverse Magnitude Freq. Response');
            plot(obj.faxis/1e9, 20*log10(abs(1./Havg)),"LineWidth",2,"Color",[0.3467    0.5360    0.6907]) ;
            xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;
            ylim([-1 15]);
            hold on;
            yline(3,'LineWidth',2,'LineStyle','--');

            %4)
            subplot(2,1,2); hold on; box on; title('Inverse Phase Freq. Response');
            plot(obj.faxis/1e9, unwrap(angle(1./Havg)),"LineWidth",2,"Color",[0.3467    0.5360    0.6907]) ;
            xlim([0.2 .5*max(obj.faxis)*1e-9]); grid on;

            %%% plot for publication
            figure(98989);hold on;box on;title('Magnitude Freq. Response');
            xlim([0.2 .5*max(obj.faxis)*1e-9]);
            ylim([-20, 2]);
            plot(obj.faxis/1e9, 20*log10(abs(Havg)),'LineWidth',2);
            grid on;


        end

        function save(obj,options)
            arguments
                obj
                options.fileName = ''
                options.savePath = ''
            end

            % Check if the fileName was provided
            if isempty(char(options.fileName))
                % If no file name is provided, prompt the user to enter a file name
                [options.fileName, filePath] = uiputfile('*.mat', 'Save As');

                % If the user cancels the dialog, fileName and filePath will be 0
                if isequal(options.fileName, 0) || isequal(filePath, 0)
                    disp('Save operation cancelled.');
                    return;
                end

                % Update the savePath with the directory chosen by the user
                options.savePath = filePath;
            end

            % Check if the savePath was provided
            if isempty(char(options.savePath))
                % If no path is provided, open a UI window to select the path
                options.savePath = uigetdir('', 'Select a folder to save the file');

                % If the user cancels the dialog, savePath will be 0
                if options.savePath == 0
                    disp('Save operation cancelled.');
                    return;
                end
            else
                % If a path is provided, validate it
                if ~isfolder(options.savePath)
                    error('The specified path does not exist.');
                end
            end

            % Construct the full file path
            fullFileName = fullfile(options.savePath, options.fileName);

            % Save the data to the specified file
            save(fullFileName, 'obj');
            fprintf('Frequency response information successfully saved to %s\n', fullFileName);

        end

        function data = load(obj, options)

            % Function to load data from a specified file and path.
            arguments
                obj
                options.fileName = ''
                options.loadPath = ''
            end

            % Check if the fileName was provided
            if isempty(char(options.fileName))
                % If no file name is provided, open a UI window to select the file
                [options.fileName, options.loadPath] = uigetfile('*.mat', 'Select a file to load');

                % If the user cancels the dialog, fileName and loadPath will be 0
                if isequal(options.fileName, 0) || isequal(options.loadPath, 0)
                    disp('Load operation cancelled.');
                    data = [];
                    return;
                end
            end

            % If loadPath is not provided or empty, use the current folder
            if isempty(char(options.loadPath))
                options.loadPath = pwd;
            else
                % Validate the path if provided
                if ~isfolder(options.loadPath)
                    error('The specified path does not exist.');
                end
            end

            % Construct the full file path
            fullFileName = fullfile(options.loadPath, options.fileName);

            % Check if the file exists
            if ~isfile(fullFileName)
                fullFileName = fullfile(options.loadPath, [char(options.fileName),'.mat']);
                if ~isfile(fullFileName)
                    error('The specified file does not exist.');
                end
            end

            % Load the data from the specified file
            loadedData = load(fullFileName);

            if ~strcmp(fieldnames(loadedData),'obj')
                    % user want to load a moveit precomp file
                    obj.H = 1./loadedData.uFF;
                    obj.faxis = loadedData.f;
            else

                % Replace whole obj here.. is this save or unsave?!
                fn = fieldnames(loadedData.obj);
                for n = 1:numel(fn)
                    try
                        obj.(fn{n}) = loadedData.obj.(fn{n});
                    end
                end

            end

            fprintf('Frequency response information successfully loaded from %s\n', fullFileName);
        end

    end

    methods (Access=private)
        % Cant be seen from outside! So put all your functions here that can/
        % shall not be called from outside

        function rOFDM = genRNDDMT(obj,randkey)

            rOFDM = NaN(1,(2*obj.Nacq+obj.Ncp+1)*obj.Navg);

            s = RandStream('mt19937ar','Seed',randkey,'NormalTransform','Polar');

            ref = 2 * round(rand(s,obj.Nacq, obj.Navg)) - 1 ;              % Navg Random BPSK sequences

            ofdm = [ones(1, obj.Navg) ; ref ; conj(ref(end:-1:1,:)) ] ;    % DMT

            ofdm = ifft(ofdm) ;                                            % Navg real OFDM sequences

            ofdm = [ofdm(end-obj.Ncp+1 : end, :) ; ofdm] ;                 % Add cyclic prefix

            rOFDM = reshape(ofdm, 1, size(ofdm,1)*size(ofdm,2)) ;

            rOFDM = rOFDM./max(abs(rOFDM(:)));

        end

        function [rH] = estHfromDMT(obj, data_in, ref_in)


            Nfft = 2*obj.Nacq + 1 ;
            data_in = reshape(data_in,1,length(data_in));
            ref_in = reshape(ref_in,1,length(ref_in));

            if length(data_in) ~= length(ref_in)
                % 0. cross-correlate the received signal with its reference to extract the periods
                corr = abs(ifft( fft(data_in(1:length(ref_in))) .* conj(fft(ref_in)) )) ;

                % find max
                [~, peak] = max(corr) ;
                peak=max(1,peak-1);

                %! Dont (circ)shift the signal here (if signals are equally long) as this would remove the phase information!
                %tested with a butterworth filter this exactly reconstructs the
                %phase and the magnitude. However, the option is here
                
                % data_in = circshift(data_in,-peak) ;
                data_in = data_in(peak:peak+(Nfft+obj.Ncp)*obj.Navg-1) ;
            end

            % 1. Reshape signal to a matrix to support noise averaging
            
            Y = reshape(data_in, Nfft+obj.Ncp, obj.Navg).' ;

            X = reshape(ref_in, Nfft + obj.Ncp, obj.Navg).' ;

            % 2. Remove cyclic prefix and apply FFT transformation
            Y = fft(Y(:, obj.Ncp+1 : obj.Ncp + Nfft), [], 2) ;
            X = fft(X(:, obj.Ncp+1 : obj.Ncp + Nfft), [], 2) ;

            % 4. Caclulate the frequency response using H = Y/X
            rH = Y./X ;

        end

    end
end