numPackets = 20;
idleTime   = 10; % In microseconds

chanBW     = 'CBW320'; % Channel bandwidth
mcs        = 12; % Modulation and coding scheme
numTx      = 1; % Number of transmit antennas
apepLength = 8000; % A-MPDU length pre-EOF padding in bytes

savedState = rng(0); % Set random state

cfgEHT = ehtMUConfig(chanBW); % EHT MU configuration object 
cfgEHT.NumTransmitAntennas = numTx;
cfgEHT.User{1}.MCS = mcs;
cfgEHT.User{1}.NumSpaceTimeStreams = numTx;
cfgEHT.User{1}.APEPLength = apepLength;
cfgEHT.RU{1}.SpatialMapping = 'Direct';

osf = 4; % Oversampling factor

psduLength = cfgEHT.getPSDULength*8; % PSDU length in bits per packet
data = randi([0 1],psduLength*numPackets,1);

txWaveform = ehtWaveformGenerator(data,cfgEHT,'NumPackets',numPackets,'IdleTime',idleTime*1e-6,'OversamplingFactor',osf);

fs = ehtSampleRate(cfgEHT);
disp(['Baseband sampling rate: ' num2str(fs/1e6) ' Msps']);

txWaveform = [zeros(round(idleTime*1e-6*fs),numTx); txWaveform];

pSaturation = 25; % Saturation power of a power amplifier in dBm
hpaBackoff = 16;  % Power amplifier backoff in dB
nonLinearity = comm.MemorylessNonlinearity;
nonLinearity.Method = 'Rapp model';
nonLinearity.Smoothness = 3; % p parameter
nonLinearity.LinearGain = -hpaBackoff;
nonLinearity.OutputSaturationLevel = db2mag(pSaturation-30); 

for i=1:cfgEHT.NumTransmitAntennas
    txWaveform(:,i) = nonLinearity(txWaveform(:,i));
end

NF = 6;         % Noise figure (dB)
BW = fs*osf;    % Bandwidth (Hz)
k = 1.3806e-23; % Boltzman constant (J/K)
T = 290;        % Ambient temperature (K)
noisePower = 10*log10(k*T*BW)+NF;

awgnChannel = comm.AWGNChannel('NoiseMethod','Variance','Variance',10^(noisePower/10));
txWaveform = awgnChannel(txWaveform);

aStop = 40; % Stopband attenuation
ofdmInfo = ehtOFDMInfo('EHT-Data',cfgEHT); % 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);

[EVMPerPkt,EVMPerSC] = evmSetup(cfgEHT);

ind = ehtFieldIndices(cfgEHT);

minPktLen = double(ind.LSTF(2)-ind.LSTF(1))+1;

rxWaveformLength = size(rxWaveform,1);
pktLength = double(ind.EHTData(2));
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
    
    % Detect packet and determine coarse packet offset
    pktOffset = ehtPacketDetect(rxWaveform,cfgEHT.ChannelBandwidth,searchOffset);
    % Packet offset from start of the waveform
    pktOffset = searchOffset+pktOffset; 
    % Skip packet if L-STF is empty
    if isempty(pktOffset) || (pktOffset<0) || ...
            ((pktOffset+ind.LSIG(2))>rxWaveformLength)
        break;
    end
    
    % Extract L-STF and perform coarse frequency offset correction
    nonht = rxWaveform(pktOffset+(ind.LSTF(1):ind.LSIG(2)),:);  
    coarsefreqOff = ehtCoarseCFOEstimate(nonht,cfgEHT.ChannelBandwidth);
    nonht = helperFrequencyOffset(nonht,fs,-coarsefreqOff);
    
    % Extract the legacy fields and determine fine packet offset
    lltfOffset = ehtSymbolTimingEstimate(nonht,cfgEHT.ChannelBandwidth);
    pktOffset = pktOffset+lltfOffset; % Determine packet offset

    % If offset is outwith bounds of the waveform, then 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 = ehtFineCFOEstimate(lltf,cfgEHT.ChannelBandwidth);
    rxPacket = helperFrequencyOffset(rxPacket,fs,-fineFreqOff);
    
    % Extract EHT-LTF samples, demodulate, and perform channel estimation
    ehtLTF = rxPacket(ind.EHTLTF(1):ind.EHTLTF(2),:);
    ehtLTFDemod = ehtDemodulate(ehtLTF,'EHT-LTF',cfgEHT);

    % Channel estimate
    [chanEst,pilotEst] = ehtLTFChannelEstimate(ehtLTFDemod,cfgEHT);
    
    % Spectral flatness measurement
    ehtTxSpectralFlatnessMeasurement(chanEst,cfgEHT,pktNum);
    
    % Data demodulate
    rxData = rxPacket(ind.EHTData(1):ind.EHTData(2),:);
    demodSym = ehtDemodulate(rxData,'EHT-Data',cfgEHT);
    
    % Pilot phase tracking
    demodSym = ehtCommonPhaseErrorTracking(demodSym,chanEst,cfgEHT);
    
    % Estimate noise power in EHT fields
    nVarEst = ehtNoiseEstimate(demodSym(ofdmInfo.PilotIndices,:,:),pilotEst,cfgEHT);

    % Extract data subcarriers from demodulated symbols and channel
    % estimate
    demodDataSym = demodSym(ofdmInfo.DataIndices,:,:);
    chanEstData = chanEst(ofdmInfo.DataIndices,:,:);

    % Equalization and STBC combining
    [eqSym,csi] = ehtEqualizeCombine(demodDataSym,chanEstData,nVarEst,cfgEHT);
 
    % Compute RMS EVM over all spatial streams for a 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
    ehtTxEVMConstellationPlots(eqSym,evmPerSC,cfgEHT,pktNum);

    % Recover data field bits
    rxPSDU = ehtDataBitRecover(eqSym,nVarEst,csi,cfgEHT);

    if isequal(rxPSDU,data((1:psduLength)+(pktNum-1)*psduLength))
        fprintf('    Decode success\n');
    else
        fprintf('    Decode failure\n');
    end

    % 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

startIdx = osf*(ind.EHTData(1)-1)+1; % Upsampled start of EHT Data
endIdx = osf*ind.EHTData(2);         % Upsampled end of EHT Data
delay = grpdelay(firinterp,1); % Group delay of downsampling filter
idx = zeros(endIdx-startIdx+1,pktNum);
for i = 1:pktNum
    % Start of packet in txWaveform
    pktOffset = round(osf*pktOffsetStore(i))-delay;
    % Indices of EHT-Data in txWaveform
    idx(:,i) = (pktOffset+(startIdx:endIdx));
end
gatedEHTTData = txWaveform(idx(:),:);

if pktNum>0
    ehtSpectralMaskTest(gatedEHTTData,fs,osf);
end