nr uplink channel state information estimation using srs -凯发k8网页登录

numframes = 1; % 10 ms frames
snr = 20; % snr in db

% create ue/carrier configuration
ue = nrcarrierconfig;
ue.nsizegrid = 52;          % bandwidth in number of resource blocks (52rbs at 15khz scs for 10mhz bw)
ue.subcarrierspacing = 15;  % 15, 30, 60, 120, 240 (khz)
ue.cyclicprefix = 'normal'; % 'normal' or 'extended'
ntxants = 2;  % number of transmit antennas (1,2,4)
nrxants = 2;  % number of receive antennas
nlayers = min(ntxants,nrxants);
% configure a periodic multi-port srs and enable frequency hopping
srs = nrsrsconfig;
srs.numsrssymbols = 4;          % number of ofdm symbols allocated per slot (1,2,4)
srs.symbolstart = 8;            % starting ofdm symbol within a slot
srs.numsrsports = ntxants;      % number of srs antenna ports (1,2,4).
srs.frequencystart = 0;         % frequency position of the srs in bwp in rbs
srs.nrrc = 0;                   % additional offset from freqstart specified in blocks of 4 prbs (0...67)
srs.csrs = 14;                  % bandwidth configuration c_srs (0...63). it controls the allocated bandwidth to the srs
srs.bsrs = 0;                   % bandwidth configuration b_srs (0...3). it controls the allocated bandwidth to the srs
srs.bhop = 0;                   % frequency hopping configuration (0...3). set bhop < bsrs to enable frequency hopping
srs.ktc = 2;                    % comb number (2,4). frequency density in subcarriers
srs.repetition = 2;             % repetition (1,2,4). it disables frequency hopping in blocks of |repetition| symbols
srs.srsperiod = [2 0];          % periodicity and offset in slots. srsperiod(2) must be < srsperiod(1)
srs.resourcetype = 'periodic';  % resource type ('periodic', 'semi-persistent','aperiodic'). use 'aperiodic' to disable inter-slot frequency hopping

practicalsynchronization = true;

csisubbandsize = 4; % number of rbs per subband

channel = nrtdlchannel;
channel.delayprofile = 'tdl-c';
channel.delayspread = 40e-9;
channel.maximumdopplershift = 30;
channel.numtransmitantennas = ntxants;
channel.numreceiveantennas = nrxants;
channel.seed = 5;
% set channel sample rate
ofdminfo = nrofdminfo(ue);
channel.samplerate = ofdminfo.samplerate;
% get the maximum delay of the channel
chinfo = info(channel);
maxchdelay = chinfo.maximumchanneldelay;
% reset random generator for reproducibility
rng('default');

% number of slots to simulate
numslots = numframes*ue.slotsperframe;
% total number of subcarriers and symbols per slot
k = ue.nsizegrid * 12;
l = ue.symbolsperslot;
% initialize arrays storing channel estimates
alltxgrid = zeros([k l*numslots ntxants]);
slotgridsize = [k l nrxants ntxants];
hest = zeros(slotgridsize);
hestinterp = zeros(slotgridsize);
hestupdate = zeros(slotgridsize);
totalgridsize = [k l*numslots nrxants ntxants];
allhest = zeros(totalgridsize);
allhestperfect = zeros(totalgridsize);
allhestinterp = zeros(totalgridsize);
% initialize noise power estimate
nvar = 0;
% calculate the number of csi subbands for the carrier
numcsisubbands = ceil(ue.nsizegrid/csisubbandsize);
% initialize sinr per subband, slot, and pmi
maxpmi = hmaxpuschprecodingmatrixindicator(nlayers,ntxants);
sinrsubband = zeros([numcsisubbands numslots maxpmi 1]);
% initialize pmi matrix and sinr loss
pmi = nan(numcsisubbands,numslots);
pmiperfect = pmi;
loss = zeros(size(pmi));
% initialize timing estimation offset
offset = chinfo.channelfilterdelay;
% calculate srs cdm lengths
cdmlengths = hsrscdmlengths(srs);
% ofdm symbols used for csi acquisition
csiselectsymbols = srs.symbolstart   (1:srs.numsrssymbols);
for nslot = 0:numslots-1
    % update slot counter
    ue.nslot = nslot;
    % generate srs and map to slot grid
    [srsindices,srsindinfo] = nrsrsindices(ue,srs);
    srssymbols = nrsrs(ue,srs);
    % create a slot-wise resource grid empty grid and map srs symbols
    txgrid = nrresourcegrid(ue,ntxants);
    txgrid(srsindices) = srssymbols;
    % determine if the slot contains srs
    issrsslot= ~isempty(srssymbols);
    % ofdm modulation
    [txwaveform,waveforminfo] = nrofdmmodulate(ue,txgrid);
    txwaveform = [txwaveform; zeros(maxchdelay, size(txwaveform,2))]; %#ok % required later to flush the channel filter to obtain the received signal
    % transmission through channel
    [rxwaveform,pathgains] = channel(txwaveform);
    % add awgn to the received time domain waveform
    % normalize noise power to take account of sampling rate, which is
    % a function of the ifft size used in ofdm modulation. the snr
    % is defined per re for each receive antenna (ts 38.101-4).
    snr = 10^(snr/10);
    n0 = 1/sqrt(2.0*nrxants*double(waveforminfo.nfft)*snr);
    noise = n0*complex(randn(size(rxwaveform)),randn(size(rxwaveform)));
    rxwaveform = rxwaveform   noise;
    % perform timing offset estimation
    pathfilters = getpathfilters(channel);
    % timing estimation is only performed in the slots where the srs is
    % transmitted to keep the perfect and practical channel estimation
    % synchronized.
    if issrsslot
        if practicalsynchronization
            % practical synchronization. correlate the received waveform
            % with the srs to give timing offset estimate
            offset = nrtimingestimate(ue,rxwaveform,srsindices,srssymbols);
        else
            offset = nrperfecttimingestimate(pathgains,pathfilters); %#ok
        end
    end
    % perform ofdm demodulation
    rxgrid = nrofdmdemodulate(ue,rxwaveform(1 offset:end,:));
    % perform practical channel estimation
    % update channel estimates only in the symbols and rbs containing srs
    % in this slot and hold the estimates from the previous slot in other
    % locations. nvar is not updated when there is no srs transmission. use
    % a time-averaging window that covers all the srs symbols transmitted.
    if issrsslot % this slot contains an srs transmission
        [hest,nvar] = nrchannelestimate(ue,rxgrid,srsindices,srssymbols,'averagingwindow',[0 7],'cdmlengths',cdmlengths);
        % use channel estimate from previous slot for ofdm symbols before the first srs symbol
        hestinterp = repmat(hestinterp(:,end,:,:),1,ue.symbolsperslot);
        % update channel estimate in ofdm symbols and rb where the srs is
        % present and hold all channel estimates until the end of the slot
        firstsymbol = srs.symbolstart 1;
        lastsymbol = srs.symbolstart   srs.numsrssymbols;
        hestupdate(:,firstsymbol:lastsymbol,:,:) = hest(:,firstsymbol:lastsymbol,:,:);
        hestupdate(:,lastsymbol:l,:,:) = repmat(hest(:,lastsymbol,:,:),1,ue.symbolsperslot-lastsymbol 1);
        idxhestupdate = hestupdate ~= 0; % indices of updated channel estimates
        hestinterp(idxhestupdate) = hestupdate(idxhestupdate);
    else % hold previous channel estimates if this slot does not contain srs
        hestinterp = repmat(hestinterp(:,end,:,:),1,ue.symbolsperslot);
    end
    % pmi selection
    % select the precoder matrix indicators for a number of layers using
    % the interpolated channel estimates in the ofdm symbols specified by
    % csiselectsymbols. the pmis are estimated per csi subband
    [pmi(:,nslot 1),sinrsubband(:,nslot 1,:),subbandindices] = hpmiselect(nlayers, hestinterp(:,csiselectsymbols,:,:), nvar, csisubbandsize);
    % pmi selection sinr loss
    % calculate the performance loss as a ratio of the sinr after precoding
    % with pmis selected using a practical channel estimate and the sinr
    % after precoding with pmis selected using a perfect channel estimate.
    % calculate perfect channel estimate for perfect pmi selection
    hestperfect = nrperfectchannelestimate(ue,pathgains,pathfilters,offset);
    % perfect noise estimate from noise realization
    noisegrid = nrofdmdemodulate(ue,noise(1 offset:end,:));
    nvarperfect = var(noisegrid(:));
    [loss(:,nslot 1),pmiperfect(:,nslot 1)] = hpmiselectionsinrloss(pmi(:,nslot 1), nlayers,  hestperfect(:,csiselectsymbols,:,:), nvarperfect);
    % save a copy of all transmitted ofdm grids and channel estimates for
    % display purposes
    thisslot = nslot*l   (1:l); % symbols of the current slot
    alltxgrid(:,thisslot,:) = txgrid;
    allhest(:,thisslot,:,:) = hest;
    allhestinterp(:,thisslot,:,:) = hestinterp;
    allhestperfect(:,thisslot,:,:) = hestperfect;
end

% create x-axis and y-axis vectors
symbols = 0:(ue.nslot 1)*ue.symbolsperslot-1;
slots = 0:ue.nslot;
subcarriers = 1:ue.nsizegrid*12;
resourceblocks = 1:ue.nsizegrid;

figure
imagesc(symbols,subcarriers,abs(alltxgrid(:,:,1,1)));
xlabel('ofdm symbol'); ylabel('subcarrier'); axis xy;
title('transmitted srs (port 1)');

% remove first ofdm symbols not containing srs to improve visualization of
% channel estimation error
hestinterp = allhestinterp;
idx = 1:(srs.srsperiod(2)*ue.symbolsperslot   srs.symbolstart);
hestinterp(:,idx,:,:) = nan;
hestinterp((srs.nrb*12 1):end,:,:,:) = nan;
hchannelestimationplot(symbols,subcarriers,allhestperfect,hestinterp);

% first expand loss from subbands into rbs for display purposes
pmirb = hexpandsubbandtorb(pmi, csisubbandsize, ue.nsizegrid);
pmiperfectrb = hexpandsubbandtorb(pmiperfect, csisubbandsize, ue.nsizegrid);
lossrb = hexpandsubbandtorb(loss, csisubbandsize, ue.nsizegrid);
hpmiplot(slots,resourceblocks,pmirb,pmiperfectrb,lossrb);

numlayers = min(size(allhestinterp,[3 4]));
if numlayers ~= 1
    pmierr = sum( abs(pmi - pmiperfect) > 0, [1 2])./ sum( ~isnan(pmi), [1 2]);
    totpmiest = sum(~isnan(pmi),[1 2]);
    fprintf('number of estimated pmi: %d \n', totpmiest);
    fprintf('    number of wrong pmi: %d \n', ceil(pmierr*totpmiest));
    fprintf('         relative error: %.1f (%%) \n', pmierr*100);
else
    fprintf('for a single layer, pmi is always 0.\n');
end

hbestsinrplot(slots,resourceblocks,sinrsubband,pmi,csisubbandsize);

% displays perfect and practical channel estimates and the channel
% estimation error for the first transmit and receive ports. the channel
% estimation error is defined as the absolute value of the difference
% between the perfect and practical channel estimates.
function hchannelestimationplot(symbols,subcarriers,allhestperfect,allhestinterp)
    figure
    subplot(311)
    imagesc(symbols, subcarriers, abs(allhestperfect(:,:,1,1)));
    axis xy; xlabel('ofdm symbol'); ylabel('subcarrier');
    colorbar;
    title('perfect channel estimate (txant=1, rxant=1)');
    subplot(312)
    imagesc(symbols, subcarriers, abs(allhestinterp(:,:,1,1)), ...
            'alphadata',~isnan(allhestinterp(:,:,1,1)))
    axis xy; xlabel('ofdm symbol'); ylabel('subcarrier');
    colorbar;
    title('srs-based practical channel estimate (txant=1, rxant=1) ');
    % display channel estimation error, defined as the difference between the
    % srs-based and perfect channel estimates
    subplot(313)
    hesterr = abs(allhestinterp - allhestperfect);
    imagesc(symbols, subcarriers, hesterr(:,:,1,1),...
            'alphadata',~isnan(hesterr(:,:,1,1)));
    axis xy; xlabel('ofdm symbol'); ylabel('subcarrier');
    colorbar;
    title('channel estimation error (txant=1, rxant=1)');
end
% displays the pmi evolution and pmi estimation sinr loss over time and
% frequency. the sinr loss is defined as a ratio of the sinr after
% precoding with estimated and perfect pmis. estimated pmis are
% obtained using a practical channel estimate and perfect pmis are
% selected using a perfect channel estimate.
function hpmiplot(slots,resourceblocks,pmirb,pmiperfectrb,lossrb)
    figure
    subplot(311)
    imagesc(slots,resourceblocks,pmiperfectrb,'alphadata',~isnan(pmiperfectrb)); axis xy;
    c = clim;
    cm = colormap;
    colormap( cm(1:floor(size(cm,1)/(c(2)-c(1)) -1):end,:) ); % adjust colormap to pmi discrete values
    colorbar
    xlabel('slot'); ylabel('resource block'), title('pmi selected using perfect channel estimates')
    subplot(312)
    imagesc(slots,resourceblocks,pmirb,'alphadata',~isnan(pmirb)); axis xy;
    colorbar,
    xlabel('slot'); ylabel('resource block'), title('pmi selected using srs')
    subplot(313)
    imagesc(slots,resourceblocks,lossrb,'alphadata',~isnan(lossrb));
    colormap(gca,cm)
    xlabel('slot'); ylabel('resource block'); axis xy; colorbar;
    title('pmi estimation sinr loss (db)')
end
% displays the sinr per resource block obtained after precoding with the
% pmi that maximizes the sinr per subband.
function hbestsinrplot(slots,resourceblocks,sinrsubband,pmi,csibandsize)
    % display sinr after precoding with best pmi
    bestsinrpersubband = nan(size(sinrsubband,[1 2]));
    % get sinr per subband and slot using best pmi
    [sb,nslot] = find(~isnan(pmi));
    for i = 1:length(sb)
        bestsinrpersubband(sb(i),nslot(i)) = sinrsubband(sb(i),nslot(i),pmi(sb(i),nslot(i)) 1);
    end
    % first expand sinr from subbands into rbs for display purposes
    bestsinrperrb = hexpandsubbandtorb(bestsinrpersubband, csibandsize, length(resourceblocks));
    figure
    sinrdb = 10*log10(abs(bestsinrperrb));
    imagesc(slots,resourceblocks,sinrdb,'alphadata',~isnan(sinrdb));
    axis xy; colorbar;
    xlabel('slot');
    ylabel('resource block')
    title('average sinr per subband and slot after precoding with best pmi (db)')
end
% expands a 2d matrix of values per subband in the first dimension into a
% matrix of values per resource block.
function rbvalues = hexpandsubbandtorb(subbandvalues, bandsize, nrb)
    lastbandsize = mod(nrb,bandsize);
    lastbandsize = lastbandsize   bandsize*(lastbandsize==0);
    rbvalues = [kron(subbandvalues(1:end-1,:),ones(bandsize,1));...
                subbandvalues(end,:).*ones(lastbandsize,1)];
end