cfgVHT = wlanVHTConfig;            % Create packet configuration
cfgVHT.ChannelBandwidth = 'CBW80'; % 80 MHz
cfgVHT.NumTransmitAntennas = 1;    % One transmit antenna
cfgVHT.NumSpaceTimeStreams = 1;    % One space-time stream
cfgVHT.STBC = false;               % No STBC so one spatial stream
cfgVHT.MCS = 8;                    % Modulation: 256-QAM
cfgVHT.APEPLength = 3000;          % A-MPDU length pre-EOF padding in bytes

osf = 3;          % 3x oversampling factor
numPackets = 20;  % Generate 20 packets
idleTime = 10e-6; % 10 microseconds idle time between packets

% Create random data; PSDULength is in bytes
savedState = rng(0); % Set random state
data = randi([0 1],cfgVHT.PSDULength*8*numPackets,1);

% Generate a multi-packet waveform
txWaveform = wlanWaveformGenerator(data,cfgVHT, ...
    'OversamplingFactor',osf,'NumPackets',numPackets,'IdleTime',idleTime);

fs = wlanSampleRate(cfgVHT); % Baseband sampling rate of the waveform

modelIQImbalance = true; % Set to true to add IQ gain and phase imbalance

if modelIQImbalance
    iqGaindB = 1; % IQ gain imbalance in dB, specify from the range [-1 1]
    iqPhaseDeg = 1; % IQ phase imbalance in degrees, specify from the range [-2 2]
    iqGainLin = db2mag(iqGaindB); % Convert gain from dB to linear value
    txWaveform = real(txWaveform) + 1i*imag(txWaveform)*iqGainLin*exp(1j*iqPhaseDeg*pi/180); % As specified in Equation-1 of [5]
end

pSaturation = 25; % Saturation power of a power amplifier in dBm
hpaBackoff = 13; % Power amplifier backoff in dB

% Create and configure a memoryless nonlinearity to model the amplifier
nonLinearity = comm.MemorylessNonlinearity;
nonLinearity.Method = 'Rapp model';
nonLinearity.Smoothness = 3; % p parameter
nonLinearity.LinearGain = -hpaBackoff;
nonLinearity.OutputSaturationLevel = db2mag(pSaturation-30);

% Apply the model to each transmit antenna
txWaveform = nonLinearity(txWaveform);

thNoise = comm.ThermalNoise('NoiseMethod','Noise Figure','SampleRate',fs*osf,'NoiseFigure',6);
for i = 1:cfgVHT.NumTransmitAntennas
    txWaveform(:,i) = thNoise(txWaveform(:,i));
end

aStop = 40; % Stopband attenuation
ofdmInfo = wlanVHTOFDMInfo('VHT-Data',cfgVHT); % OFDM parameters
SCS = fs/ofdmInfo.FFTLength; % Subcarrier spacing
txbw = max(abs(ofdmInfo.ActiveFrequencyIndices))*2*SCS; % Occupied bandwidth
[L,M] = rat(osf);
maxLM = max([L M]);
R = (fs-txbw)/fs;
TW = 2*R/maxLM; % Transition width
b = designMultirateFIR(L,M,TW,aStop);

firinterp = dsp.FIRRateConverter(M,L,b);
rxWaveform = firinterp(txWaveform);

% Setup EVM measurements
[EVMPerPkt,EVMPerSC] = vhtEVMSetup(cfgVHT);

% Get indices for accessing each field within the time-domain packet
ind = wlanFieldIndices(cfgVHT);

rxWaveformLength = size(rxWaveform,1);
pktLength = double(ind.VHTData(2));

% Define the minimum length of data we can detect; length of the L-STF in
% samples
minPktLen = double(ind.LSTF(2)-ind.LSTF(1))+1;

% Setup the measurement plots
[hSF,hCon,hEVM] = vhtTxSetupPlots(cfgVHT);

rmsEVM = zeros(numPackets,1);
pktOffsetStore = zeros(numPackets,1);

rng(savedState); % Restore random state

pktNum = 0;
searchOffset = 0; % Start at first sample (no offset)
while (searchOffset+minPktLen)<=rxWaveformLength
    % Packet detect
    pktOffset = wlanPacketDetect(rxWaveform,cfgVHT.ChannelBandwidth);
    % Packet offset from start of waveform
    pktOffset = searchOffset+pktOffset;
    % If no packet detected or offset outwith bounds of waveform then stop
    if isempty(pktOffset) || (pktOffset<0) || ...
            ((pktOffset+ind.LSIG(2))>rxWaveformLength)
        break;
    end

    % Extract non-HT fields and perform coarse frequency offset correction
    % to allow for reliable symbol timing
    nonht = rxWaveform(pktOffset+(ind.LSTF(1):ind.LSIG(2)),:);
    coarsefreqOff = wlanCoarseCFOEstimate(nonht,cfgVHT.ChannelBandwidth);
    nonht = helperFrequencyOffset(nonht,fs,-coarsefreqOff);

    % Determine offset between the expected start of L-LTF and actual start
    % of L-LTF
    lltfOffset = wlanSymbolTimingEstimate(nonht,cfgVHT.ChannelBandwidth);
    % Determine packet offset
    pktOffset = pktOffset+lltfOffset;
    % If offset is without bounds of waveform skip samples and continue
    % searching within remainder of the waveform
    if (pktOffset<0) || ((pktOffset+pktLength)>rxWaveformLength)
        searchOffset = pktOffset+double(ind.LSTF(2))+1;
        continue;
    end

    % Timing synchronization complete; extract the detected packet
    rxPacket = rxWaveform(pktOffset+(1:pktLength),:);
    pktNum = pktNum+1;
    disp(['  Packet ' num2str(pktNum) ' at index: ' num2str(pktOffset+1)]);

    % Apply coarse frequency correction to the extracted packet
    rxPacket = helperFrequencyOffset(rxPacket,fs,-coarsefreqOff);

    % Perform fine frequency offset correction on the extracted packet
    lltf = rxPacket(ind.LLTF(1):ind.LLTF(2),:); % Extract L-LTF
    fineFreqOff = wlanFineCFOEstimate(lltf,cfgVHT.ChannelBandwidth);
    rxPacket = helperFrequencyOffset(rxPacket,fs,-fineFreqOff);

    % Extract VHT-LTF samples, demodulate and perform channel estimation
    vhtltf = rxPacket(ind.VHTLTF(1):ind.VHTLTF(2),:);
    vhtltfDemod = wlanVHTLTFDemodulate(vhtltf,cfgVHT);

    % Get single stream channel estimate
    chanEstSSPilots = vhtSingleStreamChannelEstimate(vhtltfDemod,cfgVHT);

    % Channel estimate
    chanEst = wlanVHTLTFChannelEstimate(vhtltfDemod,cfgVHT);

    % Perform IQ gain and phase imbalance estimation
    [gainEst,phaseEst,alphaEst,betaEst,gamma,dataRot] = ...
                                    helperIQImbalanceEstimate(chanEst,cfgVHT);
    fprintf('    Measured IQ gain & phase imbalance: %2.2f dB, %2.2f deg\n',gainEst,phaseEst);

    % Perform IQ gain and phase imbalance correction on channel
    % estimates
    chanEst = chanEst./(alphaEst + betaEst.*gamma); % As specified in Equation-29 of [5]

    % Spectral flatness measurement
    vhtTxSpectralFlatnessMeasurement(chanEst,cfgVHT,pktNum,hSF);

    % Extract VHT Data samples from the waveform
    vhtdata = rxPacket(ind.VHTData(1):ind.VHTData(2),:);

    % Estimate the noise power in VHT data field
    noiseVarVHT = vhtNoiseEstimate(vhtdata,chanEstSSPilots,cfgVHT);

    % Extract VHT Data samples and perform OFDM demodulation, equalization
    % and phase tracking
    [~,~,eqSym] = wlanVHTDataRecover(vhtdata,chanEst,noiseVarVHT,cfgVHT,...
        'EqualizationMethod','ZF','PilotPhaseTracking','PreEQ'); % Use zero forcing algorithm for equalization

    % Perform IQ gain and phase imbalance correction on VHT data
    eqSym = eqSym.*dataRot; % Carrier rotation on data subcarriers
    eqSym = ((conj(alphaEst)*eqSym)-(betaEst*conj(eqSym(end:-1:1,:,:))))/((abs(alphaEst)^2)-(abs(betaEst)^2));  % As specified in Equation-30 of [5]

    % Compute RMS EVM over all spatial streams for packet
    rmsEVM(pktNum) = EVMPerPkt(eqSym);
    fprintf('    RMS EVM: %2.2f%%, %2.2fdB\n',rmsEVM(pktNum),20*log10(rmsEVM(pktNum)/100));

    % Compute RMS EVM per subcarrier and spatial stream for the packet
    evmPerSC = EVMPerSC(eqSym); % Nst-by-1-by-Nss

    % Plot RMS EVM per subcarrier and equalized constellation
    vhtTxEVMConstellationPlots(eqSym,evmPerSC,cfgVHT,pktNum,hCon,hEVM);

    % Store the offset of each packet within the waveform
    pktOffsetStore(pktNum) = pktOffset;

    % Increment waveform offset and search remaining waveform for a packet
    searchOffset = pktOffset+pktLength+minPktLen;
end

if pktNum>0
    fprintf('Average EVM for %d packets: %2.2f%%, %2.2fdB\n', ...
        pktNum,mean(rmsEVM(1:pktNum)),20*log10(mean(rmsEVM(1:pktNum))/100));
else
    disp('No complete packet detected');
end

% Indices for accessing each field within the time-domain packet
ind = wlanFieldIndices(cfgVHT,'OversamplingFactor', osf);
startIdx = ind.VHTData(1); % Start of Data
endIdx = ind.VHTData(2);   % End of Data
idleNSamps = idleTime*fs*osf; % Idle time samples
perPktLength = endIdx+idleNSamps;

idx = zeros(endIdx-startIdx+1,numPackets);
for i = 1:numPackets
    % Start of packet in txWaveform, accounting for the filter delay
    pktOffset = (i-1)*perPktLength;
    % Indices of non-HT Data in txWaveform
    idx(:,i) = pktOffset+(startIdx:endIdx);
end
gatedVHTData = txWaveform(idx(:),:);

if pktNum>0
    helperSpectralMaskTest(gatedVHTData,fs,osf);
end