radioFound = false;
radiolist = findsdru;
for i = 1:length(radiolist)
  if strcmp(radiolist(i).Status, 'Success')
    if strcmp(radiolist(i).Platform, 'B210')
        radio = comm.SDRuReceiver('Platform','B210', ...
                 'SerialNum', radiolist(i).SerialNum);
        radio.MasterClockRate = 1.92e6 * 4; % Need to exceed 5 MHz minimum
        radio.DecimationFactor = 4;         % Sampling rate is 1.92e6
        radioFound = true;
        break;
    end
    if (strcmp(radiolist(i).Platform, 'X300') || ...
        strcmp(radiolist(i).Platform, 'X310'))
        radio = comm.SDRuReceiver('Platform',radiolist(i).Platform, ...
                 'IPAddress', radiolist(i).IPAddress);
        radio.MasterClockRate = 184.32e6;
        radio.DecimationFactor = 96;        % Sampling rate is 1.92e6
        radioFound = true;
    end
  end
end

if ~radioFound
    error(message('sdru:examples:NeedMIMORadio'));
end

radio.ChannelMapping = [1 2];     % Receive signals from both channels
radio.CenterFrequency = 900e6;
radio.Gain = 30;
radio.SamplesPerFrame = 19200; % Sampling rate is 1.92 MHz. LTE frames are 10 ms long
radio.OutputDataType = 'double';
radio.EnableBurstMode = true;
radio.NumFramesInBurst = 4;
radio.OverrunOutputPort = true;

radio

burstCaptures = zeros(19200,4,2);

len = 0;
for frame = 1:4
    while len == 0
        [data,len,lostSamples] = step(radio);
        burstCaptures(:,frame,:) = data;
    end
    len = 0;
end
release(radio);

eNodeBOutput = reshape(burstCaptures,[],2);
sr = 1.92e6 ; % LTE sampling rate

% Check for presence of LTE System Toolbox
if isempty(ver('lte'))
    error(message('sdru:examples:NeedLST'));
end

% Prior to decoding the MIB, the UE does not know the full system
% bandwidth. The primary and secondary synchronization signals (PSS and
% SSS) and the PBCH (containing the MIB) all lie in the central 72
% subcarriers (6 resource blocks) of the system bandwidth, allowing the UE
% to initially demodulate just this central region. Therefore the bandwidth
% is initially set to 6 resource blocks. The I/Q waveform needs to be
% resampled accordingly. At this stage we also display the spectrum of the
% input signal |eNodeBOutput|.

% set up some housekeeping variables:
% separator for command window logging
separator = repmat('-',1,50);
% plots
if (~exist('channelFigure','var') || ~isvalid(channelFigure))
    channelFigure = figure('Visible','off');
end
[spectrumAnalyzer,synchCorrPlot,pdcchConstDiagram] = ...
    hSIB1RecoveryExamplePlots(channelFigure,sr);
% PDSCH EVM
pdschEVM = comm.EVM();
pdschEVM.MaximumEVMOutputPort = true;

% Set eNodeB basic parameters
enb = struct;                   % eNodeB config structure
enb.DuplexMode = 'FDD';         % assume FDD duxplexing mode
enb.CyclicPrefix = 'Normal';    % assume normal cyclic prefix
enb.NDLRB = 6;                  % Number of resource blocks
ofdmInfo = lteOFDMInfo(enb);    % Needed to get the sampling rate

if (isempty(eNodeBOutput))
    fprintf('\nReceived signal must not be empty.\n');
    return;
end

% Display received signal spectrum
fprintf('\nPlotting received signal spectrum...\n');
step(spectrumAnalyzer, awgn(eNodeBOutput, 100.0));

if (sr~=ofdmInfo.SamplingRate)
    fprintf('\nResampling from %0.3fMs/s to %0.3fMs/s for cell search / MIB decoding...\n',sr/1e6,ofdmInfo.SamplingRate/1e6);
else
    fprintf('\nResampling not required; received signal is at desired sampling rate for cell search / MIB decoding (%0.3fMs/s).\n',sr/1e6);
end
% Downsample received signal
nSamples = ceil(ofdmInfo.SamplingRate/round(sr)*size(eNodeBOutput,1));
nRxAnts = size(eNodeBOutput, 2);
downsampled = zeros(nSamples, nRxAnts);
for i=1:nRxAnts
    downsampled(:,i) = resample(eNodeBOutput(:,i), ofdmInfo.SamplingRate, round(sr));
end

fprintf('\nPerforming frequency offset estimation...\n');
% Note that the duplexing mode is set to FDD here because timing synch has
% not yet been performed - for TDD we cannot use the duplexing arrangement
% to indicate which time periods to use for frequency offset estimation
% prior to doing timing synch.
delta_f = lteFrequencyOffset(setfield(enb,'DuplexMode','FDD'), downsampled); %#ok<SFLD>
fprintf('Frequency offset: %0.3fHz\n',delta_f);
downsampled = lteFrequencyCorrect(enb, downsampled, delta_f);

% Cell search to find cell identity and timing offset
fprintf('\nPerforming cell search...\n');
[cellID, offset] = lteCellSearch(enb, downsampled);

% Compute the correlation for each of the three possible primary cell
% identities; the peak of the correlation for the cell identity established
% above is compared with the peak of the correlation for the other two
% primary cell identities in order to establish the quality of the
% correlation.
corr = cell(1,3);
for i = 0:2
    enb.NCellID = mod(cellID + i,504);
    [~,corr{i+1}] = lteDLFrameOffset(enb, downsampled);
    corr{i+1} = sum(corr{i+1},2);
end
threshold = 1.3 * max([corr{2}; corr{3}]); % multiplier of 1.3 empirically obtained
if (max(corr{1})<threshold)
    warning('sdru:examples:WeakSignal','Synchronization signal correlation was weak; detected cell identity may be incorrect.');
end
enb.NCellID = cellID;

% plot PSS/SSS correlation and threshold
synchCorrPlot.YLimits = [0 max([corr{1}; threshold])*1.1];
step(synchCorrPlot, [corr{1} threshold*ones(size(corr{1}))]);

% perform timing synchronisation
fprintf('Timing offset to frame start: %d samples\n',offset);
downsampled = downsampled(1+offset:end,:);
enb.NSubframe = 0;

% show cell-wide settings
fprintf('Cell-wide settings after cell search:\n');
disp(enb);

% Channel estimator configuration
cec.PilotAverage = 'UserDefined';     % Type of pilot averaging
cec.FreqWindow = 9;                   % Frequency window size
cec.TimeWindow = 9;                   % Time window size
cec.InterpType = 'cubic';             % 2D interpolation type
cec.InterpWindow = 'Centered';        % Interpolation window type
cec.InterpWinSize = 1;                % Interpolation window size

% Assume 4 cell-specific reference signals for initial decoding attempt;
% ensures channel estimates are available for all cell-specific reference
% signals
enb.CellRefP = 4;

fprintf('Performing OFDM demodulation...\n\n');

griddims = lteResourceGridSize(enb); % Resource grid dimensions
L = griddims(2);                     % Number of OFDM symbols in a subframe
% OFDM demodulate signal
rxgrid = lteOFDMDemodulate(enb, downsampled);
if (isempty(rxgrid))
    fprintf('After timing synchronization, signal is shorter than one subframe so no further demodulation will be performed.\n');
    return;
end
% Perform channel estimation
if (strcmpi(enb.DuplexMode,'TDD'))
    enb.TDDConfig = 0;
    enb.SSC = 0;
end
[hest, nest] = lteDLChannelEstimate(enb, cec, rxgrid(:,1:L,:));

% Decode the MIB
% Extract resource elements (REs) corresponding to the PBCH from the first
% subframe across all receive antennas and channel estimates
fprintf('Performing MIB decoding...\n');
pbchIndices = ltePBCHIndices(enb);
[pbchRx, pbchHest] = lteExtractResources( ...
    pbchIndices, rxgrid(:,1:L,:), hest(:,1:L,:,:));

% Decode PBCH
[bchBits, pbchSymbols, nfmod4, mib, enb.CellRefP] = ltePBCHDecode( ...
    enb, pbchRx, pbchHest, nest);

% Parse MIB bits
enb = lteMIB(mib, enb);

% Incorporate the nfmod4 value output from the function ltePBCHDecode, as
% the NFrame value established from the MIB is the System Frame Number
% (SFN) modulo 4 (it is stored in the MIB as floor(SFN/4))
enb.NFrame = enb.NFrame+nfmod4;

% Display cell wide settings after MIB decoding
fprintf('Cell-wide settings after MIB decoding:\n');
disp(enb);

if (enb.CellRefP==0)
    fprintf('MIB decoding failed (enb.CellRefP=0).\n\n');
    return;
end
if (enb.NDLRB==0)
    fprintf('MIB decoding failed (enb.NDLRB=0).\n\n');
    return;
end

fprintf('Restarting reception now that bandwidth (NDLRB=%d) is known...\n',enb.NDLRB);

% Resample now we know the true bandwidth
ofdmInfo = lteOFDMInfo(enb);
if (sr~=ofdmInfo.SamplingRate)
    fprintf('\nResampling from %0.3fMs/s to %0.3fMs/s...\n',sr/1e6,ofdmInfo.SamplingRate/1e6);
else
    fprintf('\nResampling not required; received signal is at desired sampling rate for NDLRB=%d (%0.3fMs/s).\n',enb.NDLRB,sr/1e6);
end
nSamples = ceil(ofdmInfo.SamplingRate/round(sr)*size(eNodeBOutput,1));
resampled = zeros(nSamples, nRxAnts);
for i = 1:nRxAnts
    resampled(:,i) = resample(eNodeBOutput(:,i), ofdmInfo.SamplingRate, round(sr));
end

% Perform frequency offset estimation and correction
fprintf('\nPerforming frequency offset estimation...\n');
% Note that the duplexing mode is set to FDD here because timing synch has
% not yet been performed - for TDD we cannot use the duplexing arrangement
% to indicate which time periods to use for frequency offset estimation
% prior to doing timing synch.
delta_f = lteFrequencyOffset(setfield(enb,'DuplexMode','FDD'), resampled); %#ok<SFLD>
fprintf('Frequency offset: %0.3fHz\n',delta_f);
resampled = lteFrequencyCorrect(enb, resampled, delta_f);

% Find beginning of frame
fprintf('\nPerforming timing offset estimation...\n');
offset = lteDLFrameOffset(enb, resampled);
fprintf('Timing offset to frame start: %d samples\n',offset);
% aligning signal with the start of the frame
resampled = resampled(1+offset:end,:);

% OFDM demodulation
fprintf('\nPerforming OFDM demodulation...\n\n');
rxgrid = lteOFDMDemodulate(enb, resampled);

% Check this frame contains SIB1, if not advance by 1 frame provided we
% have enough data, terminate otherwise.
if (mod(enb.NFrame,2)~=0)
    if (size(rxgrid,2)>=(L*10))
        rxgrid(:,1:(L*10),:) = [];
        fprintf('Skipping frame %d (odd frame number does not contain SIB1).\n\n',enb.NFrame);
    else
        rxgrid = [];
    end
    enb.NFrame = enb.NFrame + 1;
end

% Advance to subframe 5, or terminate if we have less than 5 subframes
if (size(rxgrid,2)>=(L*5))
    rxgrid(:,1:(L*5),:) = [];   % Remove subframes 0 to 4
else
    rxgrid = [];
end
enb.NSubframe = 5;

if (isempty(rxgrid))
    fprintf('Received signal does not contain a subframe carrying SIB1.\n\n');
end

% Reset the HARQ buffers
decState = [];

% While we have more data left, attempt to decode SIB1
while (size(rxgrid,2) > 0)

    fprintf('%s\n',separator);
    fprintf('SIB1 decoding for frame %d\n',mod(enb.NFrame,1024));
    fprintf('%s\n\n',separator);

    % Reset the HARQ buffer with each new set of 8 frames as the SIB1
    % info may be different
    if (mod(enb.NFrame,8)==0)
        fprintf('Resetting HARQ buffers.\n\n');
        decState = [];
    end

    % Extract current subframe
    rxsubframe = rxgrid(:,1:L,:);

    % Perform channel estimation
    [hest,nest] = lteDLChannelEstimate(enb, cec, rxsubframe);

    % PCFICH demodulation, CFI decoding. The CFI is now demodulated and
    % decoded using similar resource extraction and decode functions to
    % those shown already for BCH reception. lteExtractResources is used to
    % extract REs corresponding to the PCFICH from the received subframe
    % rxsubframe and channel estimate hest.
    fprintf('Decoding CFI...\n\n');
    pcfichIndices = ltePCFICHIndices(enb);  % Get PCFICH indices
    [pcfichRx, pcfichHest] = lteExtractResources(pcfichIndices, rxsubframe, hest);
    % Decode PCFICH
    cfiBits = ltePCFICHDecode(enb, pcfichRx, pcfichHest, nest);
    cfi = lteCFIDecode(cfiBits); % Get CFI
    if (isfield(enb,'CFI') && cfi~=enb.CFI)
        release(pdcchConstDiagram);
    end
    enb.CFI = cfi;
    fprintf('Decoded CFI value: %d\n\n', enb.CFI);

    % For TDD, the PDCCH must be decoded blindly across possible values of
    % the PHICH configuration factor m_i (0,1,2) in TS36.211 Table 6.9-1.
    % Values of m_i = 0, 1 and 2 can be achieved by configuring TDD
    % uplink-downlink configurations 1, 6 and 0 respectively.
    if (strcmpi(enb.DuplexMode,'TDD'))
        tddConfigs = [1 6 0];
    else
        tddConfigs = 0; % not used for FDD, only used to control while loop
    end
    dci = {};
    while (isempty(dci) && ~isempty(tddConfigs))
        % Configure TDD uplink-downlink configuration
        if (strcmpi(enb.DuplexMode,'TDD'))
            enb.TDDConfig = tddConfigs(1);
        end
        tddConfigs(1) = [];
        % PDCCH demodulation. The PDCCH is now demodulated and decoded
        % using similar resource extraction and decode functions to those
        % shown already for BCH and CFI reception
        pdcchIndices = ltePDCCHIndices(enb); % Get PDCCH indices
        [pdcchRx, pdcchHest] = lteExtractResources(pdcchIndices, rxsubframe, hest);
        % Decode PDCCH and plot constellation
        [dciBits, pdcchSymbols] = ltePDCCHDecode(enb, pdcchRx, pdcchHest, nest);
        step(pdcchConstDiagram, pdcchSymbols);

        % PDCCH blind search for System Information (SI) and DCI decoding.
        % The LTE System Toolbox provides full blind search of the PDCCH to
        % find any DCI messages with a specified RNTI, in this case the
        % SI-RNTI.
        fprintf('PDCCH search for SI-RNTI...\n\n');
        pdcch = struct('RNTI', 65535);
        dci = ltePDCCHSearch(enb, pdcch, dciBits); % Search PDCCH for DCI
    end

    % If DCI was decoded, proceed with decoding PDSCH / DL-SCH
    if ~isempty(dci)

        dci = dci{1};
        fprintf('DCI message with SI-RNTI:\n');
        disp(dci);
        % Get the PDSCH configuration from the DCI
        [pdsch, trblklen] = hPDSCHConfiguration(enb, dci, pdcch.RNTI);
        pdsch.NTurboDecIts = 5;
        fprintf('PDSCH settings after DCI decoding:\n');
        disp(pdsch);

        % PDSCH demodulation and DL-SCH decoding to recover SIB bits.
        % The DCI message is now parsed to give the configuration of the
        % corresponding PDSCH carrying SIB1, the PDSCH is demodulated and
        % finally the received bits are DL-SCH decoded to yield the SIB1
        % bits.

        fprintf('Decoding SIB1...\n\n');
        % Get PDSCH indices
        [pdschIndices,pdschIndicesInfo] = ltePDSCHIndices(enb, pdsch, pdsch.PRBSet);
        [pdschRx, pdschHest] = lteExtractResources(pdschIndices, rxsubframe, hest);
        % Decode PDSCH
        [dlschBits,pdschSymbols] = ltePDSCHDecode(enb, pdsch, pdschRx, pdschHest, nest);
        % Decode DL-SCH with soft buffer input/output for HARQ combining
        if ~isempty(decState)
            fprintf('Recombining with previous transmission.\n\n');
        end
        [sib1, crc, decState] = lteDLSCHDecode(enb, pdsch, trblklen, dlschBits, decState);

        % Compute PDSCH EVM
        recoded = lteDLSCH(enb, pdsch, pdschIndicesInfo.G, sib1);
        remod = ltePDSCH(enb, pdsch, recoded);
        [~,refSymbols] = ltePDSCHDecode(enb, pdsch, remod);
        release(pdschEVM);
        [rmsevm,peakevm] = step(pdschEVM, refSymbols{1}, pdschSymbols{1});
        fprintf('PDSCH RMS EVM: %0.3f%%\n',rmsevm);
        fprintf('PDSCH Peak EVM: %0.3f%%\n\n',peakevm);

        fprintf('SIB1 CRC: %d\n',crc);
        if crc == 0
            fprintf('Successful SIB1 recovery.\n\n');
        else
            fprintf('SIB1 decoding failed.\n\n');
        end

    else
        % indicate that DCI decoding failed
        fprintf('DCI decoding failed.\n\n');
    end

    % Update channel estimate plot
    figure(channelFigure);
    surf(abs(hest(:,:,1,1)));
    hSIB1RecoveryExamplePlots(channelFigure);
    channelFigure.CurrentAxes.XLim = [0 size(hest,2)+1];
    channelFigure.CurrentAxes.YLim = [0 size(hest,1)+1];

    % Skip 2 frames and try SIB1 decoding again, or terminate if we
    % have less than 2 frames left.
    if (size(rxgrid,2)>=(L*20))
        rxgrid(:,1:(L*20),:) = [];   % Remove 2 more frames
    else
        rxgrid = []; % Less than 2 frames left
    end
    enb.NFrame = enb.NFrame+2;

end

figure('Position',[0 0 800 600])
plot(abs(hest(:,1,1,1)),'b-')
hold on
plot(abs(hest(:,1,2,1)),'r-.')
plot(abs(hest(:,1,1,2)),'k--')
plot(abs(hest(:,1,2,2)),'m.')
hold off
grid on
m = max(reshape(abs(hest(:,1,:,:)),[],1));
axis([0 size(hest,1)+1 0 m*1.3])
title('Channel Estimate for First OFDM Symbol')
xlabel('Subcarrier Index')
ylabel('Magnitude of Estimated Channel Coefficient')
legend('TX1 to RX1','TX1 to RX2', ...
       'TX2 to RX1','TX2 to RX2', ...
       'Location', 'northeast')