function [gainEst, phaseEst, alphaEst, betaEst, gamma, dataRot] = ...
                                helperIQImbalanceEstimate(chanEst, cfgVHT)
% helperIQImbalanceEstimate IQ gain and phase imbalance estimation using
% the channel estimates
%
%   [GAINEST, PHASEEST, ALPHAEST, BETAEST, GAMMA, DATAROT] =
%   helperIQImbalanceEstimate(CHANEST, CFGVHT) estimates the IQ gain
%   imbalance, GAINEST and phase imbalance, PHASEEST using the channel
%   estimates, CHANEST as specified in [1].
%
%   GAINEST is the estimated IQ gain imbalance in dB.
%
%   PHASEEST is the estimated IQ phase imbalance in degrees.
%
%   ALPHAEST and BETAEST are the IQ imbalance coefficients derived from IQ
%   gain and phase imbalance, GAINEST and PHASEEST respectively.
%
%   GAMMA is the ratio of conjugate mirrored and actual Very High Throughput Long
%   Training Field (VHT-LTF) sequence (e.g. conj(wrev(vhtLTF))./vhtLTF) as
%   specified in equation 7 of [1].
%
%   DATAROT is the phase rotation on data subcarriers: IEEE Std 802.11ac-2013 Section 22.3.7.5.
%
%   CHANEST is the estimated channel using VHT-LTF sequence.
%
%   CFGVHT is a packet format configuration object of type <a
%   href="matlab:help('wlanVHTConfig')">wlanVHTConfig</a>.
%
%   See also helperIQImbalance.

%   Reference:
%   [1] M. Janaswamy, N. K. Chavali and S. Batabyal, "Measurement of
%   transmitter IQ parameters in HT and VHT wireless LAN systems," 2016
%   International Conference on Signal Processing and Communications
%   (SPCOM), Bangalore.

%   Copyright 2021 The MathWorks, Inc.

% Get the known VHT-LTF sequence
chanBW = cfgVHT.ChannelBandwidth;
numSTSVec = cfgVHT.NumSpaceTimeStreams;
ltf = wlan.internal.vhtltfSequence(chanBW,numSTSVec);
cfgOFDM = wlan.internal.wlanGetOFDMConfig(chanBW,'Long','VHT',numSTSVec);
ind = sort([cfgOFDM.DataIndices; cfgOFDM.PilotIndices]);
vhtLTF = ltf(ind).*cfgOFDM.CarrierRotations(ind);
gamma = conj(vhtLTF(end:-1:1))./vhtLTF; % As specified in Equation-7 of [1]
dataRot = cfgOFDM.CarrierRotations(cfgOFDM.DataIndices); % Carrier rotation: IEEE Std 802.11ac-2013 Section 22.3.7.5.

% Estimate the IQ gain and phase imbalances as specified in Section-III of [1]
if strcmp(chanBW,'CBW40')
    indices = find(diff(imag(gamma))~=0);
    chanRatio = chanEst(indices,1,1)./chanEst(indices+1,1,1); % Equation-16 of [1]
    % Estimate the IQ gain and phase imbalances using the indices
    % corresponding to the gamma transitions from -j to +j.
    % IQ gain and phase imbalance is averaged over different values of
    % frequency indices as specified in Equation-26 & 27 of [1].
    chanRatio1 = chanRatio(1:2:end);
    m4 = tan(angle(chanRatio1));
    gain = m4+sqrt(m4.^2+1);
    b = (gain.^2+1).*(1-abs(chanRatio1).^2)./(2*gain.*(1+abs(chanRatio1).^2));
    phaseEst1 = mean(asind(b));
    approxError = cosd(phaseEst1); % Approximation error as per Equation-24 of [1]
    m4 = m4*approxError;
    gain = (m4+sqrt(m4.^2+1));
    gainEst1 = mean(gain);

    % Estimate the IQ gain and phase imbalances using the indices
    % corresponding to the gamma transitions from +j to -j.
    % IQ gain and phase imbalance is averaged over different values of
    % frequency indices as specified in Equation-23 & 25 of [1].
    chanRatio2 = chanRatio(2:2:end);
    m3 = tan(angle(chanRatio2));
    gain = -m3+sqrt(m3.^2+1);
    b = (gain.^2+1).*(-1+abs(chanRatio2).^2)./(2*gain.*(1+abs(chanRatio2).^2));
    phaseEst2 = mean(asind(b));

    approxError = cosd(phaseEst2); % Approximation error as per Equation-22 of [1]
    m3 = m3*approxError;
    gain = (-m3+sqrt(m3.^2+1));
    gainEst2 = mean(gain);
else
    % Estimate the IQ gain and phase imbalances using the indices
    % corresponding to the gamma transitions from -1 to +1
    indices = find(diff((gamma))>0);
    chanRatio = chanEst(indices,1,1)./chanEst(indices+1,1,1); % Equation-8 of [1]
    % IQ gain and phase imbalance is averaged over different values of
    % frequency indices as specified in Equation-14 & 15 of [1]
    gainEst1 = mean(abs(chanRatio));
    phaseEst1 = mean(angle(chanRatio))*180/pi;

    % Estimate the IQ gain and phase imbalances using the indices
    % corresponding to the gamma transitions from +1 to -1
    indices = find(diff((gamma))<0);
    chanRatio = chanEst(indices+1,1,1)./chanEst(indices,1,1);
    % IQ gain and phase imbalance is averaged over different values of
    % frequency indices as specified in Equation-11 & 12 of [1]
    gainEst2 = mean(abs(chanRatio));
    phaseEst2 = mean(angle(chanRatio))*180/pi;
end

% Over all IQ gain and phase imbalance estimates
gainEstLin = (gainEst1 + gainEst2)/2; % Gain imbalance in linear scale
gainEst = 20*log10(gainEstLin); % Gain imbalance in dB
phaseEst = (phaseEst1 + phaseEst2)/2; % Phase imbalance in degrees

alphaEst = (1+gainEstLin*exp(1i*phaseEst*pi/180))/2; % As specified in Equation-4 of [1]
betaEst = (1-gainEstLin*exp(1i*phaseEst*pi/180))/2;
end