%% Joint Sampling Rate and Carrier Frequency Offset Tracking
% This example demonstrates joint sampling rate and carrier frequency
% offset tracking in a WLAN receiver.

% Copyright 2016-2018 The MathWorks, Inc.

%% Introduction
% In a WLAN radio a single oscillator is typically used to derive clocks
% for sampling and modulation. The oscillators in the transmitter and
% receiver radios do not run at the exact same frequency. Due to this
% mismatch, the sampling instants at the receiver shift relative to the
% transmitter. Therefore, a sample rate offset (SRO) exists between
% receiver and transmitter. Similarly, a carrier frequency offset (CFO)
% exists between the receiver and transmitter due to mismatched carrier
% frequencies. The inclusion of pilot subcarriers in the IEEE(R) 802.11(TM)
% standard allows for tracking and correction of SRO and CFO impairments.
%
% In OFDM systems SRO manifests itself as a subcarrier and symbol dependent
% phase rotation and inter-carrier interference (ICI) [ <#15 1> ]. When the
% SRO is large, and a packet is long, subcarriers far away from DC will
% experience a substantial impairment. CFO manifests itself as ICI and a
% symbol dependent phase rotation common to all subcarriers. The phase
% rotation on subcarriers from one OFDM symbol to the next due to these
% impairments is illustrated below. $\Phi_k$ is the phase error, $k$ is the
% subcarrier index, $K$ is the number of subcarriers. $\zeta$ is the SRO,
% ${\Delta}f$ is the carrier frequency offset, $T_u$ is the period of the
% symbol, $\delta$ is the phase error gradient (PEG), and $\omega$ is the
% common phase error (CPE). The PEG and CPE can be used to estimate SRO and
% residual CFO.
% 
% <<../PilotTracking.png>>
%
% In this example, an IEEE 802.11ac(TM) VHT waveform is generated with
% fixed SRO and CFO impairments [ <#15 2> ]. The impaired waveform is
% synchronized and the data portion demodulated and decoded with and
% without joint timing and phase tracking to correct for SRO and CFO. The
% equalized constellation with and without tracking is shown to demonstrate
% the effectiveness of tracking.

%% Generate a Baseband Waveform
% A VHT configuration object is created to parameterize the transmission.
% Because the impairments are easier to visualize with fewer OFDM symbols,
% a data payload with only 500 bytes is transmitted using 16-QAM
% modulation.

% Create a VHT configuration
cfgVHT = wlanVHTConfig;
cfgVHT.ChannelBandwidth = 'CBW20';
cfgVHT.NumTransmitAntennas = 1;
cfgVHT.NumSpaceTimeStreams = 1;
cfgVHT.MCS = 4;          % 16-QAM and 3/4 coding rate
cfgVHT.APEPLength = 500; % Bytes

% Create a random PSDU
s = rng(10); % Seed the random number generator
psdu = randi([0 1],cfgVHT.PSDULength*8,1,'int8');

% Generate a VHT packet
tx = wlanWaveformGenerator(psdu,cfgVHT);

%% Model Impairments
% Sample rate offset between the transmitter and receiver is modeled by
% resampling the transmitted waveform. The <matlab:doc('resample')
% resample> function can be used to model a limited range of large sample
% rate offsets. In this example a -100 parts per million (PPM) sample rate
% offset is modeled.

% Model sample rate offset
p = 1e4;   % Interpolation factor
q = 1e4-1; % Decimation factor
sro = (1-p/q)*1e6;
disp('Impairments:');
disp(['  Sample rate offset (SRO): ' num2str(sro,'%.1f') ' PPM']);

% Resample the waveform, appending zeros to allow for filter delay
N = 100; % Size of filter used for resampling
rx = resample([tx; zeros(N,cfgVHT.NumTransmitAntennas)],p,q,N);

%%
% Residual carrier frequency offset is added to the waveform using the
% helper function <matlab:edit('helperFrequencyOffset')
% helperFrequencyOffset>. In this example we assume the same oscillator is
% used for sampling and modulation, therefore the CFO will be a function of
% the SRO and carrier frequency.

fc = 5.25e9;         % Carrier frequency, Hertz
cfo = (sro*1e-6)*fc; % Carrier frequency offset, Hertz
disp(['  Carrier frequency offset (CFO): ' num2str(cfo,'%.1f') ' Hz']);

fs = wlanSampleRate(cfgVHT);           % Baseband sample rate
rx = helperFrequencyOffset(rx,fs,cfo); % Add frequency offset

%%
% Noise is added to the waveform with 30 dBW variance.

awgnChannel = comm.AWGNChannel('NoiseMethod','Variance','Variance',10^(-30/10));
rx = awgnChannel(rx);

%% Front-End Synchronization and Receiver Processing
% The following processing steps occur to synchronize the packet, in
% preparation for recovering the data field.
%
% # The packet is detected
% # Coarse carrier frequency offset is estimated and corrected
% # Symbol timing synchronization is established
% # Fine carrier frequency offset is estimated and corrected
% # The L-LTF is OFDM demodulated and noise estimation is performed
% # The VHT-LTF is OFDM demodulated and channel estimation is performed

% Generate field indices
ind = wlanFieldIndices(cfgVHT);

% Packet detection
tOff = wlanPacketDetect(rx,cfgVHT.ChannelBandwidth);

% Coarse frequency offset correction
lstf = rx(tOff+(ind.LSTF(1):ind.LSTF(2)),:);
coarseCFOEst = wlanCoarseCFOEstimate(lstf,cfgVHT.ChannelBandwidth);
rx = helperFrequencyOffset(rx,fs,-coarseCFOEst);

% Symbol timing synchronization
nonhtPreamble = rx(tOff+(ind.LSTF(1):ind.LSIG(2)),:);
symOff = wlanSymbolTimingEstimate(nonhtPreamble,cfgVHT.ChannelBandwidth);
tOff = tOff+symOff;

% Fine frequency offset correction
lltf = rx(tOff+(ind.LLTF(1):ind.LLTF(2)),:);
fineCFOEst = wlanFineCFOEstimate(lltf,cfgVHT.ChannelBandwidth);
rx = helperFrequencyOffset(rx,fs,-fineCFOEst);

% Noise power estimation
lltf = rx(tOff+(ind.LLTF(1):ind.LLTF(2)),:);
lltfDemod = wlanLLTFDemodulate(lltf,cfgVHT);
noiseEst = helperNoiseEstimate(lltfDemod,cfgVHT.ChannelBandwidth,cfgVHT.NumSpaceTimeStreams);

% Channel estimation
vhtltf = rx(tOff+(ind.VHTLTF(1):ind.VHTLTF(2)),:);
vhtltfDemod = wlanVHTLTFDemodulate(vhtltf,cfgVHT);
chanEst = wlanVHTLTFChannelEstimate(vhtltfDemod,cfgVHT);

%% Recovery Without Sample Rate Offset or Residual CFO Tracking
% The coarse and fine frequency offset estimation and correction removes
% the majority of CFO, but residual CFO remains due to the presence of
% impairments in the waveform. This must be tracked and corrected by the
% receiver.

disp('Front-end impairment correction:');
frontEndCFOEst = coarseCFOEst+fineCFOEst;
disp(['  Estimated CFO: ' num2str(frontEndCFOEst,'%.1f') ' Hz']);
residualCFO = cfo-frontEndCFOEst;
disp(['  Residual CFO after initial correction: ' num2str(residualCFO,'%.1f') ' Hz']);

%%
% The helper function <matlab:edit('trackingVHTDataRecover')
% trackingVHTDataRecover> recovers the VHT data field with optional pilot
% tracking to correct for timing and phase errors due to SRO and CFO. Pilot
% tracking is controlled using the helper object
% <matlab:edit('trackingRecoveryConfig') trackingRecoveryConfig>.
%
% The data field is first recovered without pilot tracking. The data field
% is extracted from the waveform using the start and end sample indices of
% the field at the baseband rate. If the receiver sampling rate is higher
% than the transmitter rate, the receiver requires more samples than the
% transmitter produces. To allow for this, |Ne| additional samples are
% extracted from the waveform and passed to the recovery function. The
% maximum number of additional samples required is a function of the
% expected SRO, the baseband sampling rate, and the maximum packet
% duration.

% Recovery configuration with pilot tracking disabled
cfgRec = trackingRecoveryConfig;
cfgRec.PilotTracking = 'None';

% Extract data field with Ne additional samples to allow for negative SRO
maxDuration = 5.484e-3; % Maximum packet duration in seconds
maxSRO = 120;           % PPM
Ne = ceil(fs*maxDuration*maxSRO*1e-6); % Number of extra samples
dataInd = tOff+(ind.VHTData(1):ind.VHTData(2)+Ne);
dataInd = dataInd(dataInd<=length(rx)); % Only use indices within waveform
data = rx(dataInd,:); 

% Perform demodulation and decoding
[rxPSDUNoTrack,~,eqSymNoTrack] = trackingVHTDataRecover(data,chanEst,noiseEst,cfgVHT,cfgRec);

%%
% The equalized constellation is plotted which shows a rotation of all
% constellation points caused by residual CFO, and a spreading of
% constellation points due to SRO. Despite the modest AWGN added to the
% waveform, the impairments cause bit errors within the decoded PSDU.

ConstNoTrack = comm.ConstellationDiagram;
ConstNoTrack.Title = 'Equalized symbols with no pilot tracking';
ConstNoTrack.ReferenceConstellation = wlanReferenceSymbols(cfgVHT);
ConstNoTrack(eqSymNoTrack(:));

[~,berNoTrack] = biterr(rxPSDUNoTrack,psdu);
disp('Bit error rate:');
disp(['  Without tracking: ' num2str(berNoTrack)]);

%% Recovery With Sample Rate Offset Tracking and Residual CFO Tracking
% Now the data field is recovered with joint timing and phase pilot
% tracking to correct for SRO and residual CFO.
%
% The tracking algorithm used in this example estimates absolute values of
% $\delta$ and $\omega$ per OFDM symbol, and applies a per subcarrier and
% symbol phase correction to the demodulated symbols to reverse the phase
% errors caused by SRO and CFO. The phase error between each received pilot
% subcarrier and the expected value is calculated per symbol and averaged
% over |PilotTrackingWindow| OFDM symbols. From this, least-square
% estimates of $\delta$ and $\omega$ are calculated per symbol. These
% estimates are used to apply a phase correction to each symbol and
% subcarrier [ <#15 3>, <#15 4> ].

% Recovery configuration with pilot tracking enabled
cfgRec = trackingRecoveryConfig;
cfgRec.PilotTracking = 'Joint'; % Joint timing and phase tracking
cfgRec.PilotTrackingWindow = 9; % Averaging window in OFDM symbols

% Perform demodulation and decoding
[rxPSDU,~,eqSymTrack,cpe,peg] = trackingVHTDataRecover(data,chanEst,noiseEst,cfgVHT,cfgRec);

%%
% The equalized constellation is plotted which shows a clear 16-QAM
% constellation with no spreading or rotation. There are no bit errors.

ConstTrack = comm.ConstellationDiagram;
ConstTrack.Title = 'Equalized symbols with joint pilot tracking';
ConstTrack.ReferenceConstellation = wlanReferenceSymbols(cfgVHT);
ConstTrack(eqSymTrack(:));

[~,berTrack] = biterr(rxPSDU,psdu);
disp(['  With tracking: ' num2str(berTrack)]);

%%
% The function <matlab:edit('trackingVHTDataRecover')
% trackingVHTDataRecover> returns measurements from which the residual CFO,
% and SRO can be estimated:
%
% * |cpe| - The common phase error (radians) per symbol
% * |peg| - The phase error gradient (radians per subcarrier) per symbol
%
% The SRO and residual CFO are estimated from these measurements using a
% linear least-square fit of the rate of change. The measurements are
% plotted using the helper function
% <matlab:edit('trackingPlotSROCFOEstimates') trackingPlotSROCFOEstimates>.

[residualCFOEst,sroEst] = trackingPlotSROCFOEstimates(cpe,peg,cfgVHT);

% Display estimated SRO, residual CFO and total CFO
fprintf('Tracked impairments:\n');
fprintf('  Estimated residual CFO: %3.1f Hz (%.1f Hz error)\n', ...
    residualCFOEst,residualCFOEst-residualCFO);
fprintf('  Estimated SRO: %3.1f PPM (%.1f PPM error)\n',sroEst,sroEst-sro);
cfoEst = frontEndCFOEst+residualCFOEst; % Initial + tracked CFO estimate
fprintf('Estimated CFO (initial + tracked): %.1f Hz (%.1f Hz error)\n',cfoEst,cfoEst-cfo);

rng(s); % Restore the state of the random number generator

%% Conclusion
% The example shows how sample rate and carrier frequency offsets can be
% tracked and corrected for while recovering the data field of a WLAN
% waveform.
%
% Data field recovery functions with joint pilot tracking for VHT, HT-MF
% and non-HT formats are provided, along with an object to configure the
% recovery algorithms:
%
% * <matlab:edit('trackingVHTDataRecover') trackingVHTDataRecover.m>
% * <matlab:edit('trackingHTDataRecover') trackingHTDataRecover.m>
% * <matlab:edit('trackingNonHTDataRecover') trackingNonHTDataRecover.m>
% * <matlab:edit('trackingRecoveryConfig') trackingRecoveryConfig.m>

%% Appendix
% This example uses the following helper functions:
%
% * <matlab:edit('helperFrequencyOffset') helperFrequencyOffset.m>
% * <matlab:edit('helperNoiseEstimate') helperNoiseEstimate.m>
% * <matlab:edit('trackingPlotSROCFOEstimates') trackingPlotSROCFOEstimates.m>
% * <matlab:edit('trackingRecoveryConfig') trackingRecoveryConfig.m>
% * <matlab:edit('trackingVHTDataRecover') trackingVHTDataRecover.m>

%% Selected Bibliography
% # Speth, Michael, et al. "Optimum receiver design for wireless broad-band
% systems using OFDM." IEEE Transactions on communications 47.11 (1999):
% 1668-1677.
% # IEEE Std 802.11ac(TM)-2013 IEEE Standard for Information technology -
% Telecommunications and information exchange between systems - Local and
% metropolitan area networks - Specific requirements - Part 11: Wireless
% LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications -
% Amendment 4: Enhancements for Very High Throughput for Operation in Bands
% below 6 GHz.
% # Chiueh, Tzi-Dar, Pei-Yun Tsai, and I-Wei Lai. Baseband receiver design
% for wireless MIMO-OFDM communications. John Wiley & Sons, 2012.
% # Horlin, P.F. and Bourdoux, A. Digital Compensation for Analog
% Front-Ends: A New Approach to Wireless Transceiver Design. Wiley, 2008.
