classdef cipCube < handle
    %cipCube is a class to do coherent processing for the radar
    %   the dataCube consists of a set of data fully describes an coherent
    %   processing interval for the radar.  Objective is lightweight and
    %   fast.

    properties(SetAccess=protected,GetAccess=public)
        myRadar;
        myRadarIO;
        dataCube;
        magSqCube;
        magSqValid;
        magSqCubeSum; % nci for Wx
        magSqCubeNum; % nci for Wx
        cipMetaData;
        dfBank;
        calCorr;
        LrChans;
        SrChans;
        sumLr;
        sumSr;
        LrPc;
        SrPc;
        eof;
        stcCorrValues;
        dopplerCells;
        bgEst;
        hannWgts;    
    end

    methods
        % constructor...
        function obj = cipCube(myRadar)
            %dataCube Construct an instance of this class
            %   Detailed explanation goes here
            obj.myRadar=myRadar;

            obj.cipMetaData=myRadar.getCipMetaData();

            singleEx=zeros(1,1,'single');
            Cplx=complex(singleEx,singleEx);
            obj.dataCube=zeros(myRadar.getNumChan(),...
                myRadar.getPriSize(),...
                myRadar.getNumPris(),...
                'like',Cplx);
            obj.magSqCube=zeros(myRadar.getNumChan(),...
                myRadar.getPriSize(),...
                myRadar.getNumPris(),'single');
            obj.magSqValid=false;

            obj.sumLr=uint32(radChan.sumLr);
            obj.sumSr=uint32(radChan.sumSr);

            load("CalCorr.mat","calCorr");
            obj.calCorr=cast(calCorr,'single');
            if (obj.cipMetaData.wxMode)
                load("dfbwts_wx.mat","dfbank");
            else
                load("dfbwts_uas.mat","dfbank");
            end
            obj.dfBank=cast(fliplr(fftshift(dfbank,2)),'single');
            obj.eof = false;
            [pcSll,pcNbar]=myRadar.getPcWeightConfig();


            obj.LrPc = PulseCompressor(waveformTypes.LfmUpChirp,...
                obj.cipMetaData.bandWidth(uint32(radChan.sumLr)),...
                obj.cipMetaData.pulseLength_s(uint32(radChan.sumLr)),...
                obj.cipMetaData.sampleRate_sps(uint32(radChan.sumLr)),...
                obj.cipMetaData.priSize,...
                pcSll,pcNbar);

            obj.SrPc = PulseCompressor(waveformTypes.LfmUpChirp,...
                obj.cipMetaData.bandWidth(uint32(radChan.sumSr)),...
                obj.cipMetaData.pulseLength_s(uint32(radChan.sumSr)),...
                obj.cipMetaData.sampleRate_sps(uint32(radChan.sumSr)),...
                obj.cipMetaData.priSize,...
                pcSll,pcNbar);
            obj.stcCorrValues=obj.myRadar.getStcCorr();
            obj.dopplerCells=obj.myRadar.getDopCells();
            obj.hannWgts=cast(hanning(obj.cipMetaData.numPris),'single');
        end

        %%
        function [cipMetaData,eof] = loadCip(obj,cipNum,cipTime)
            %loadCip -- load a CIP from the system or
            %record file.
            % arguments:
            %   optional: cipNum -- record files only jump to the cip
            %   number identified.
            %   optional: cipTime -- jump to the nearest PRI to the time of
            %   the cip (used to make real-time videos to match 30
            %   frames/second).
            % optional cipNum argument for selecting a cip from a record
            % file.
            if exist('cipNum','var') && ~isempty(cipNum)
                [obj.dataCube,obj.cipMetaData,obj.eof]=obj.myRadarIO.getCIP(obj.dataCube,cipNum);
            elseif exist('cipTime','var')
                [obj.dataCube,obj.cipMetaData,obj.eof]=obj.myRadarIO.getCIP(obj.dataCube,[],cipTime);
            else
                try
                    [obj.dataCube,obj.cipMetaData,obj.eof]=obj.myRadarIO.getCIP(obj.dataCube);
                catch
                    obj.eof = true; 
                    eof = obj.eof;
                end
            end
            cipMetaData = obj.cipMetaData;
            obj.magSqValid=false;
            eof=obj.eof;
        end
        %%
        function calibrate(obj)

            if ~obj.cipMetaData.wxMode % no need to calibrate in wx mode...
                obj.dataCube = obj.dataCube.*obj.calCorr;
                obj.magSqValid=false;
            end
        end
        % not sure this helps... due to the uncertian delays between
        % tx/rx... it is going to be hard to align the corrections to the
        % data.  without alignmnent I don't think it is a good idea to
        % apply PTM.
        %%
        function stcCorr(obj)
            obj.dataCube = obj.dataCube.*obj.stcCorrValues;
            obj.magSqValid=false;
        end
        %%
        function compress(obj)
            % need to adjust for rx/tx bias BEFORE pulse compression MPS.
            %obj.dataCube(obj.LrChans,:,:) = circshift(obj.dataCube(obj.LrChans,:,:),-4,2);
            %obj.dataCube(obj.LrChans,1:65,:) = 0; 
            obj.dataCube(obj.LrChans,:,:)=obj.LrPc.compress(obj.dataCube(obj.LrChans,:,:));
            obj.dataCube(obj.SrChans,:,:)=obj.SrPc.compress(obj.dataCube(obj.SrChans,:,:));
            obj.magSqValid=false;
        end
        %% perform doppler filtering
        function dopplerFilter(obj)
            [aa, bb , ~]=size(obj.dataCube);
            obj.dataCube=reshape(obj.dataCube,(aa*bb),[]);
            switch obj.myRadar.configVals.dopplerType
                % Byron and Peter decided that they prefered a hanning
                % weighted FFT over the descrite dft filter bank.
                case 'hanningFft'
                    obj.dataCube=obj.dataCube.*obj.hannWgts';
                    obj.dataCube=fftshift(fft(obj.dataCube,obj.cipMetaData.numPris,2),2)/sqrt(obj.cipMetaData.numPris);
                case 'dft'
                    obj.dataCube = obj.dataCube * obj.dfBank;
                otherwise
                    error("Unknown Doppler");
            end
            obj.dataCube=reshape(obj.dataCube,aa,bb,[]);
            obj.magSqValid=false;
        end
        %% calculate the magnitude squared....
        function magSquared(obj)
            obj.magSqCube=obj.dataCube.*conj(obj.dataCube);
            obj.magSqValid=true;
        end
        %% sum the magnitude squared values (non-coherent integration)
        function [finePsd,nBurst] = magSquaredSum(obj,initFlag)
            if (initFlag == 0)
                obj.magSqCubeSum = obj.dataCube.*conj(obj.dataCube);
                obj.magSqCubeNum = 1;
            else
                obj.magSqCubeSum = obj.magSqCubeSum + obj.dataCube.*conj(obj.dataCube);
                obj.magSqCubeNum = obj.magSqCubeNum + 1;
            end
            nBurst = obj.magSqCubeNum;
            finePsd = obj.magSqCubeSum ./ nBurst;
        end
        %% sum the magnitude squared values (non-coherent integration)
        function [finePsd,nBurst] = magSquaredSum_alt(obj,initFlag)
            if (~obj.magSqValid)
                obj.magSquared()
            end
            if (initFlag == 0)
                obj.magSqCubeSum = obj.magSqCube;
                obj.magSqCubeNum = 1;
            else
                obj.magSqCubeSum = obj.magSqCubeSum + obj.magSqCube;
                obj.magSqCubeNum = obj.magSqCubeNum + 1;
            end
            nBurst = obj.magSqCubeNum;
            finePsd = obj.magSqCubeSum ./ nBurst;
        end

        %% normalize the magnitude squared (to the background)
        % normalize to the mean of the doppler cells... this assumes the
        % peaks are small compared to the number of cells and will not skew
        % the value of the backgrond
        function normalize(obj)
            procChans=[uint32(radChan.sumLr) uint32(radChan.sumSr)];

            if (~obj.magSqValid)
                obj.magSquared()
            end      
            switch obj.myRadar.configVals.normalizerType
                case 'meanBg'
                    norm_range_int= [obj.cipMetaData.cells.rangeZero(radChan.sumSr) + ...
                        fix(obj.myRadar.configVals.normalizerMinRange_m / ...
                        obj.myRadar.cipMetaData.rangeCellSize_m(radChan.sumSr)): ...
                        obj.myRadar.cipMetaData.cells.procEnd(radChan.sumSr)];
                    lBgEst =mean(obj.magSqCube(procChans(2),norm_range_int,obj.dopplerCells),3);
                    obj.bgEst(2) = mean(lBgEst);
                    norm_range_int= [obj.cipMetaData.cells.rangeZero(radChan.sumLr) + ...
                        fix(obj.myRadar.configVals.normalizerMinRange_m / ...
                        obj.myRadar.cipMetaData.rangeCellSize_m(radChan.sumLr)): ...
                        obj.myRadar.cipMetaData.cells.procEnd(radChan.sumLr)];
                    lBgEst =mean(obj.magSqCube(procChans(1),norm_range_int,obj.dopplerCells),3);
                    obj.bgEst(1) = mean(lBgEst);

                    for ii=obj.LrChans
                        obj.magSqCube(ii,:,:)=obj.magSqCube(ii,:,:)./obj.bgEst(1);
                    end
                    for ii=obj.SrChans
                        obj.magSqCube(ii,:,:)=obj.magSqCube(ii,:,:)./obj.bgEst(2);
                    end
                case 'doppler'
                    obj.bgEst =mean(obj.magSqCube(procChans,:,obj.dopplerCells),3);
                    for ii=obj.LrChans
                        obj.magSqCube(ii,:,:)=obj.magSqCube(ii,:,:)./obj.bgEst(1,:);
                    end
                    for ii=obj.SrChans
                        obj.magSqCube(ii,:,:)=obj.magSqCube(ii,:,:)./obj.bgEst(2,:);
                    end
                case 'dopplerShear'
                    shearCount = obj.myRadar.configVals.normalizerShearCount;
                    [maxVal,~]=maxk(obj.magSqCube(procChans,:,obj.dopplerCells),...
                        shearCount,3);
                    obj.bgEst=(sum(obj.magSqCube(procChans,:,obj.dopplerCells),3) -...
                        sum(maxVal,3))./ (obj.cipMetaData.numPris-shearCount);
                    for ii=obj.LrChans
                        obj.magSqCube(ii,:,:)=obj.magSqCube(ii,:,:)./obj.bgEst(1,:);
                    end
                    for ii=obj.SrChans
                        obj.magSqCube(ii,:,:)=obj.magSqCube(ii,:,:)./obj.bgEst(2,:);
                    end
                case 'range'
                    error('rangeNorm not implemented'); %TODO
                case 'rangeShear'
                    error('rangeShear not implemented'); % TODO
                case 'none'
            end
        end
        %% get the data from the cipCube...
        function [aCipCube,aCipMetaData,aEof] = getCipCube(obj)
            %obj.dataCube = circshift(obj.dataCube,-7,2);
            aCipCube=obj.dataCube;
            
            aCipMetaData=obj.cipMetaData;
            aEof = obj.eof;
        end
        %% get the data from the cipMagSqCube...
        function [aMagSqCube,aCipMetaData,aEof] = getMagSqCube(obj)
            if (~obj.magSqValid)
                obj.magSquared()
            end
            aMagSqCube=obj.magSqCube;
            aCipMetaData=obj.cipMetaData;
            aEof = obj.eof;
        end

        %% perform DC bias correction (necessary because of hardware truncation issue?)
        function correctDcBias(obj)
            obj.dataCube=obj.dataCube + 0.25 + 1i*0.25;
        end
        function returnVal = threshold(~,~)
            %null
            returnVal.range=[];
            returnVal.doppler=[];
        end
        %% get the calculated thresholds for the current CIP... must be called after executing threshold function for the CIP.
        function thresholds=getThresholds(~,~)
            thresholds.peakThreshold=0;
            thresholds.cfar=0;
            thresholds.slb=0;
            thresholds.ranges=0;
        end

        function setCipNum(obj, aCipNum)
            obj.cipMetaData.cipNum=aCipNum;
        end
    end

    methods(Static)
        function obj=factory(radarObj, UDP, varargin)
            thisRadarIO = radarIO.factory(radarObj, UDP, varargin{:});
            radarObj.setClutterCells(thisRadarIO.clutterCells);
            radarObj.setDopplerCells(thisRadarIO.dopplerCells);
            radarObj.setCipMetaData(thisRadarIO.getCipMetaData());
            if thisRadarIO.cipMetaData.wxMode
                obj=wxCipCube(radarObj);
            else
                obj=uasCipCube(radarObj);
            end
            obj.myRadarIO=thisRadarIO;
            radarObj.setCipMetaData(obj.cipMetaData)
            obj.myRadar=radarObj;
            obj.cipMetaData=obj.myRadarIO.getCipMetaData();
        end
    end
end

