function [CNo,st] = estCNo( st, I, Q )

%==============================================================================%
%                                                                              %
%   CNo = estCNo( msPerDataBit, I, Q )                                         %
%                                                                              %
%   Estimate C/No for a sequence of I & Q correlator output values             %
%   David S. De Lorenzo                                                        %
%   09 November 2005                                                           %
%                                                                              %
%------------------------------------------------------------------------------%
%                                                                              %
%   Based on the method described by A.J. Van Dierendonck in "GPS Receivers",  %
%   Ch. 8 of Global Positioning System: Theory and Applications, vol. 1, 1996. %
%   The equations are 104, 105, 106, 116, and 117 on pages 391-392.            %
%                                                                              %
%   Note that the only difficult part of applying this method is determining   %
%   the navigation data bit transition boundary so that 20 ms of I&Q samples   %
%   can be grouped together without a navigation data bit transition.          %
%   To make this function most general, it is assumed that the data passed     %
%   in starts on the first sample (of 20) in a 20 ms navigation data bit.      %
%                                                                              %
%==============================================================================%
%                       GPS Software Receiver in Matlab                        %
%     Copyright (c) 2003-2017 David S. De Lorenzo / 2018-2023 Trimble Inc.     %
%==============================================================================%

M = st.msPerDataBit;   % Number of I & Q samples per symbol

if st.method==1
    if M~=20
        error("Only handle 20ms so far");
    end
    
    I5ms = sum(reshape(I,[4,5]),2);
    Q5ms = sum(reshape(Q,[4,5]),2);
    if st.cno >= 36
        % Stinger rotates I/Q out for dynamics
        cIQ = I5ms + 1i*Q5ms;
        cIQ = cIQ .* exp(-1i*angle(cIQ));
        I5ms = real(cIQ);
        Q5ms = imag(cIQ);
    end
    I10ms = sum(reshape(I5ms,[2,2]),2);
    Q10ms = sum(reshape(Q5ms,[2,2]),2);
    
    wpow = sum(I5ms.^2 + Q5ms.^2)/4;
    npow = sum(I10ms.^2 + Q10ms.^2)/8;
    
    noise_pow = wpow - npow;
    signal_pow = 2*npow - wpow;
    
    if isnan(st.signal_power)
        st.signal_power = signal_pow;
        st.noise_power = noise_pow;
        st.count = 1;
    else
        st.count = st.count+1;
        ks = 1/8;
        st.signal_power = (1-ks)*st.signal_power + ks*signal_pow;
        if st.count < 20
            kn = 1/8;
        else
            % Stinger uses 1/127, but.. we accept a bit noisier
            % SNR for faster response times here.
            kn = 1/16;
        end
        st.noise_power = (1-kn)*st.noise_power + kn*noise_pow;
    end
    CNo = 10*log10(st.signal_power/st.noise_power) + 10*log10(100);
    st.cno = CNo;
else
    T = 0.001;          % Correlation time per sample, [ms]
    
    % The number of symbol periods for C/No computation
    kMax = floor(length(I) / M);

    % Dimension arrays for wide-band power and narrow-band power
    WBP = zeros(1, kMax);
    NBP = zeros(1, kMax);
    
    % Loop over symbols in the tracked data record
    for kk = 1:kMax
        
        % Determine the array index range for each symbol
        idx = (1:M) + (kk - 1) * M;
        
        % Determine the wide-band power
        WBP(kk) = sum(I(idx).^2 + Q(idx).^2);
        
        % Determine the narrow-band power
        NBP(kk) = sum(I(idx))^2 + sum(Q(idx))^2;
        
    end
    
    % Determine the normalized power
    NP = NBP ./ WBP;
    
    % Determine the mean normalized power
    muNP = mean(NP);
    
    % Estimate the carrier-to-noise-density ratio
    CNo = 10 * log10(1/T * ((muNP - 1) / (M - muNP)));
end

%===============================================================================
