%cellQualitySearch LTE cell search + RSRQ reporting
%   REPORT = cellQualitySearch(WAVEFORM,SR) detects LTE cells in the input
%   waveform WAVEFORM with sampling rate SR and returns a structure array
%   REPORT (one element per LTE cell) containing the following fields:
%   FrequencyOffset - The frequency offset obtained by lteFrequencyOffset
%                     after resampling the waveform to 1.92Ms/s
%   NCellID         - Physical layer cell identity
%   TimingOffset    - The timing offset of the first frame head in the 
%                     input waveform after resampling to 1.92Ms/s
%   RSRQdB          - The Reference Signal Received Quality (RSRQ) value in 
%                     dB according to TS 36.214 Section 5.1.3.
%   ReportedRSRQ    - The RSRQ measurement report (integer between 0 and 34
%                     ) according to TS 36.133 Section 9.1.7.
%
%   A structure is present in the array for each detected LTE cell, where a
%   cell is considered to be detected if ReportedRSRQ>0 (therefore a value
%   of ReportedRSRQ=0 will never be returned). The RSRQ is measured over
%   each subframe of the first complete frame of the waveform (or the
%   remainder of the waveform after the first frame head, if shorter than a
%   frame).
%
%   If SR is greater than 1.92Ms/s, WAVEFORM will be resampled to 1.92Ms/s.
%   If SR is equal to 1.92M/s, no resampling will occur. If SR is less than
%   1.92Ms/s, the function will error. Note that the returned timing offset
%   is w.r.t. the resampled waveform.
%
%   REPORT = cellQualitySearch(WAVEFORM,SR,ENB) provides control over the
%   assumed duplexing mode as specified in the structure ENB, which
%   must include the following field:
%   DuplexMode      - Duplex mode ('FDD','TDD')
%
%   This function currently assumes normal cyclic prefix length.
%
%   Example:
%   enb = lteRMCDL('R.4');
%   enb.TotSubframes = 20;
%   enb.NCellID = 42;
%   [waveform,~,info] = lteRMCDLTool(enb,[]);
%   waveform = circshift(waveform,999);
%   report = cellQualitySearch(waveform,info.SamplingRate)
%   
%   report = 
% 
%     struct with fields:
% 
%    FrequencyOffset: 0
%            NCellID: 42
%       TimingOffset: 999
%             RSRQdB: -5.177477803272027
%       ReportedRSRQ: 29
%
%   See also lteOFDMInfo, lteFrequencyOffset, lteFrequencyCorrect,
%   lteCellSearch, lteDLFrameOffset, lteOFDMDemodulate.

%   Copyright 2017-2018 The MathWorks, Inc.

function report = cellQualitySearch(waveform,sr,enbIn)
    
    % initialise empty report
    report = [];
    
    % assumed eNodeB configuration
    enb.NDLRB = 6;
    enb.CyclicPrefix = 'Normal';
    if nargin <3
        enb.DuplexMode = 'FDD';
    else
        enb.DuplexMode = enbIn.DuplexMode;
    end
    
    % resampling
    ofdmInfo = lteOFDMInfo(enb);
    frameLength = ofdmInfo.SamplingRate * 0.01;
    if (sr<ofdmInfo.SamplingRate)
        error('Input sampling rate (%0.3fMs/s) must be at least %0.3fMs/s.',sr/1e6,ofdmInfo.SamplingRate/1e6);
    end
    waveform = resample(waveform,ofdmInfo.SamplingRate,round(sr));
    
    % frequency synch
    % Note: Frequency synch is performed assuming FDD (i.e. assuming that
    % all OFDM symbols in a frame are active) because timing synch has not
    % yet been achieved, so for TDD the pattern of active and inactive OFDM
    % symbols w.r.t. the start of the waveform is not known. The
    % time-domain PSS detection used in lteCellSearch below is robust to
    % small frequency offsets but larger offsets may degrade the PSS
    % correlation, so it is advantageous to perform frequency synch first
    % in case a large frequency offset is present. The assumption of FDD
    % during frequency synch will mean that correlations for inactive OFDM
    % symbols for TDD will add noise, but this should be outweighed by the
    % benefit of applying frequency synch prior to timing synch, at least
    % for larger frequency offsets.
    delta_f = lteFrequencyOffset(setfield(enb,'DuplexMode','FDD'),waveform); %#ok<SFLD>
    waveform = lteFrequencyCorrect(enb,waveform,delta_f);
    
    % cell search
    searchalg.SSSDetection = 'PostFFT';
    searchalg.MaxCellCount = 5;
    cellIDs = lteCellSearch(enb,waveform,searchalg);
    
    % for each cell:
    for c = 1:numel(cellIDs)
    
        enb.NCellID = cellIDs(c);
        
        % timing synch
        offset = lteDLFrameOffset(enb,waveform);
        synched = waveform(1+offset:end,:);
        enb.NSubframe = 0;
        synched = synched(1:min(frameLength,size(synched,1)),:);
        
        % OFDM demodulation
        rxgrid = lteOFDMDemodulate(enb,synched);
        dims = lteDLResourceGridSize(enb,1);
        L = dims(2);
        nsf = floor(size(rxgrid,2)/L);
        
        % RS measurements, taking the maximum RSRQ over hypotheses of
        % CellRefP = 1 and 2 and hypotheses of each uplink/downlink
        % configuration (TDDConfig) for TDD. Any cases where more than half
        % of the subframes have zero reported RSRQ are considered to have
        % zero reported RSRQ overall.
        maxRSRQdB = 0;
        maxReportedRSRQ = 0;
        for cellrefp = [1 2]
            sfRSRQ = zeros(1,nsf);
            sfReportedRSRQ = zeros(1,nsf);
            enb.CellRefP = cellrefp;
            for sf = 1:nsf
                enb.NSubframe = mod(sf-1,10);
                sfgrid = rxgrid(:,(sf-1)*L + (1:L),:);
                % Note: For TDD duplex mode, RS measurement is performed
                % assuming FDD in order that an RS measurement is obtained
                % for every subframe (later filtered on TDDConfig)
                rsmeas = hRSMeasurements(setfield(enb,'DuplexMode','FDD'),sfgrid); %#ok<SFLD>
                sfRSRQ(sf) = rsmeas.RSRQ;
                sfReportedRSRQ(sf) = rsrqReport(rsmeas.RSRQdB);
            end
            if (strcmpi(enb.DuplexMode,'TDD'))
                tddConfigs = 0:6;
            else
                tddConfigs = 0; % not used for FDD, only used to control for loop
            end
            for tddConfig = tddConfigs
                duplexingInfo = arrayfun(@(x)lteDuplexingInfo(setfield(setfield(enb,'TDDConfig',tddConfig),'NSubframe',mod(x-1,10))),1:nsf);
                active = [duplexingInfo.NSymbolsDL]~=0;
                % Note: TDDConfig detection may be possible by taking into
                % consideration how weak the RSRQ is for inactive subframes
                % versus how strong it is for active subframes
                if (sum(sfReportedRSRQ(active)==0)>(length(sfReportedRSRQ(active))/2))
                    reportedRSRQ = 0;
                else
                    RSRQdB = 10*log10(mean(sfRSRQ(active)));
                    reportedRSRQ = rsrqReport(RSRQdB);
                end
                if (reportedRSRQ > maxReportedRSRQ)
                    maxRSRQdB = RSRQdB;
                    maxReportedRSRQ = reportedRSRQ;
                end
            end
        end
        
        % report cell if detected
        if (maxReportedRSRQ~=0)
            report = [report struct('FrequencyOffset',delta_f,'NCellID',enb.NCellID,'TimingOffset',offset,'RSRQdB',maxRSRQdB,'ReportedRSRQ',maxReportedRSRQ)]; %#ok<AGROW>
        end
        
    end
    
end

% RSRQ measurement report according to TS 36.133 Section 9.1.7
function value = rsrqReport(RSRQdB)
    
    if (RSRQdB<-19.5)
        value = 0;
    elseif (RSRQdB>=-3)
        value = 34;
    else
        value = floor(2*(RSRQdB+19.5)+1);
    end
    
end

%hRSMeasurements TS36.214 Reference Signal measurements.
%   MEAS = hRSMeasurements(ENB,GRID) performs the measurements defined in
%   TS36.214 Sections 5.1.1 and 5.1.3 (Reference Signal Received Power
%   (RSRP), Received Signal Strength Indicator (RSSI) and Reference Signal
%   Received Quality (RSRQ)) on the input resource array GRID given
%   cell-wide settings structure ENB which must include the fields:
%   NDLRB           - Number of downlink resource blocks 
%   NCellID         - Physical layer cell identity
%   CellRefP        - Number of cell-specific reference signal antenna 
%                     ports (1,2,4)
%   CyclicPrefix    - Optional. Cyclic prefix length 
%                     ('Normal'(default),'Extended')
%   DuplexMode      - Optional. Duplex mode ('FDD'(default),'TDD')
%   Only required for 'TDD' duplex mode:
%       TDDConfig       - Optional. Uplink/Downlink Configuration (0...6) 
%                         (default 0)
%       SSC             - Optional. Special Subframe Configuration (0...8) 
%                         (default 0)
%       NSubframe       - Subframe number   
%
%   GRID is a 3-dimensional array of the resource elements for one or more
%   subframes across all receive antennas. GRID is an M-by-N-by-R array
%   where M is the number of subcarriers, N is the number of OFDM symbols
%   and R is the number of receive antennas. Dimension M must be 12*NDLRB
%   where NDLRB must be (6...110). Dimension N must be a multiple of number
%   of symbols in a subframe L, where L=14 for normal cyclic prefix and
%   L=12 for extended cyclic prefix.
%
%   MEAS is a structure including the fields:
%   AntennaRSRP         - The linear Reference Signal Received Power (RSRP)
%                         for each receive antenna.
%   AntennaRSSI         - The linear Received Signal Strength Indicator
%                         (RSSI) for each receive antenna.
%   AntennaRSRQ         - The linear Reference Signal Received Quality 
%                         (RSRQ) for each receive antenna.
%   RSRP                - The reported RSRP value, the maximum of the 
%                         linear RSRPs for each receive antenna, i.e. the 
%                         maximum element of AntennaRSRP.
%   RSSI                - The reported RSSI value, the RSSI of the receive 
%                         antenna with the greatest RSRQ.
%   RSRQ                - The reported RSRQ value, the maximum of the 
%                         linear RSRQs for each receive antenna, i.e. the 
%                         maximum element of AntennaRSRQ. 
%   RSRPdBm             - The RSRP expressed in deciBels relative to 1mW in
%                         1ohm. 
%   RSSIdBm             - The RSSI expressed in deciBels relative to 1mW in
%                         1ohm. 
%   RSRQdB              - The RSRQ expressed in deciBels.
%
function meas = hRSMeasurements(enb,grid)

    if (enb.CellRefP>=2)
        ports = [0 1];
    else
        ports = 0;
    end

    dims = lteResourceGridSize(enb);
    K = dims(1);
    L = dims(2);
    R = size(grid,3);
    nsf = floor(size(grid,2)/L);
    cellRSIndices = [];
    for p = 1:length(ports)
        cellRSIndices = [cellRSIndices; lteCellRSIndices(enb,ports(p))-(p-1)*K*L]; %#ok<AGROW>
    end                
    [ksubs,lsubs] = ind2sub([K L],reshape(double(cellRSIndices),[],length(ports)));
    ksubs(ksubs>K/2) = ksubs(ksubs>K/2) + 1;
    cellRSSyms = unique(lsubs(:,1));
    rssiIndices = repmat((cellRSSyms-1)*K,1,K)+repmat(1:K,length(cellRSSyms),1);
    rssiIndices = reshape(rssiIndices,numel(rssiIndices),1);
    meas.AntennaRSRP = NaN(nsf,R);    
    meas.AntennaRSSI = NaN(nsf,R);    
    meas.AntennaRSRQ = NaN(nsf,R);  
    startSubframe = enb.NSubframe;
    for n = 1:nsf
        enb.NSubframe = mod(startSubframe+n-1,10);
        duplexingInfo = lteDuplexingInfo(enb);
        if (strcmpi(duplexingInfo.SubframeType,'Downlink'))
            subframe = grid(:,((n-1)*L)+(1:L),:);
            cellRS = [];
            for p = 1:length(ports)
                cellRS = [cellRS; lteCellRS(enb,ports(p))]; %#ok<AGROW>
            end
            for r=1:R        
                subframe(:,:,r) = linearPhaseEqualize(subframe(:,:,r),cellRSIndices,cellRS,K,ksubs,lsubs);
                meas.AntennaRSRP(n,r) = abs(mean(subframe(cellRSIndices+(r-1)*K*L).*conj(cellRS))*length(ports)).^2;
                meas.AntennaRSSI(n,r) = sum(abs(subframe(rssiIndices+(r-1)*K*L).^2))/(length(cellRSSyms));   
                meas.AntennaRSRQ(n,r) = meas.AntennaRSRP(n,r)*double(enb.NDLRB)./meas.AntennaRSSI(n,r);
            end
        end
    end
    meas.AntennaRSRP = perAntennaMeanFn(meas.AntennaRSRP);
    meas.AntennaRSSI = perAntennaMeanFn(meas.AntennaRSSI);
    meas.AntennaRSRQ = perAntennaMeanFn(meas.AntennaRSRQ);
    idxRSRP = find(meas.AntennaRSRP==max(meas.AntennaRSRP),1);
    idxRSRQ = find(meas.AntennaRSRQ==max(meas.AntennaRSRQ),1);
    meas.RSRP = mean(meas.AntennaRSRP(idxRSRP));
    meas.RSSI = mean(meas.AntennaRSSI(idxRSRQ));    
    meas.RSRQ = mean(meas.AntennaRSRQ(idxRSRQ));    
    meas.RSRPdBm = 10*log10(meas.RSRP)+30;
    meas.RSSIdBm = 10*log10(meas.RSSI)+30; 
    meas.RSRQdB = 10*log10(meas.RSRQ);
    
end

function out = linearPhaseEqualize(in,cellRSIndices,cellRS,K,ksubs,lsubs)
    
    nports = size(ksubs,2);
    cellRSIndices = reshape(cellRSIndices,[],nports);
    cellRS = reshape(cellRS,[],nports);
    cellRSSyms = unique(lsubs);
    nsyms = length(cellRSSyms);
    thetapoly = zeros(nsyms*nports,2);
    rms = zeros(nsyms*nports,1);
    for p = 1:nports
        lsubsp = lsubs(:,p);
        for li = 1:nsyms
            l = cellRSSyms(li);
            theta = unwrap(angle(in(cellRSIndices(lsubsp==l,p)).*conj(cellRS(lsubsp==l,p))));
            [thetapoly(li + (p-1)*nsyms,:),s] = polyfit(ksubs(lsubsp==l,p),theta,1);
            rms(li + (p-1)*nsyms) = s.normr/sqrt(numel(theta)-1);
        end
    end
    threshold = 1.2;
    if (sum(rms<threshold)>0)
        thetapoly = mean(thetapoly(rms<threshold,:),1);
        k = [1:K/2 ((K/2)+2):(K+1)].';
        theta = polyval(thetapoly,k);
    else
        theta = 0;
    end
    out = in .* exp(-1i*theta);
    
end

function x = perAntennaMeanFn(r)

    x = arrayfun(@(p)mean(r(~isnan(r(:,p)),p)),1:size(r,2));
    x(isnan(x)) = 0;
    
end
