利用matlab和SDR实现LTE信号的采集以及帧同步，MIB解码

        找到这篇文章的同志们，一定正在LTE生死线上苦苦挣扎，大家都知道对LTE信号进行采集的时候，是采用30.72M（15K✖2048）的时钟频率，然后呢，LTE频带的带宽是18M（15K✖1200）为了避免大家少走弯路，我将代码附在文章末尾。我们可以采集到的信号如下图所示。

        大家可以看到这与事实是相符合的，在中心频率的左右各是±9MHz，因此我们可以通过查找移动，联通，电信对应的下行频段信息，然后就可以采集到信息了。

        帧同步如下所示。

频偏估计以及MIB就不一一展示了，大家可以运行程序调参。

 

%% Connect to Radio

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
    if (strcmp(radiolist(i).Platform, 'N300') || ...
        strcmp(radiolist(i).Platform, 'N310'))
        radio = comm.SDRuReceiver('Platform',radiolist(i).Platform, ...
                 'IPAddress', radiolist(i).IPAddress);
        radio.MasterClockRate = 122.88e6;
        radio.DecimationFactor = 64;        % Sampling rate is 1.92e6
        radioFound = true;
    end
    if (strcmp(radiolist(i).Platform, 'N320/N321'))
        radio = comm.SDRuReceiver('Platform',radiolist(i).Platform, ...
                 'IPAddress', radiolist(i).IPAddress);
        radio.MasterClockRate = 245.76e6;
        radio.DecimationFactor = 128;        % 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

%% Capture Signal

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 Toolbox
if isempty(ver('lte')) 
    error(message('sdru:examples:NeedLST'));
end

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

%% Frequency offset estimation and correction
fprintf('\nPerforming frequency offset estimation...\n');
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 and Synchronization
% Cell search to find cell identity and timing offset
fprintf('\nPerforming cell search...\n');
[cellID, offset] = lteCellSearch(enb, downsampled);

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 
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);

%% OFDM Demodulation and Channel Estimation  

% 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,:));

%% PBCH Demodulation, BCH Decoding, MIB parsing

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); 


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

%% OFDM Demodulation on Full Bandwidth

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');

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);
resampled = resampled(1+offset:end,:);   

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

%% SIB1 Decoding

% 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);    
    
    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);   
  
    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) = [];        
      
        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);

       
        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);

      

        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')