%% Compare Practical and Perfect Channel Estimates
% Generate physical broadcast channel (PBCH) demodulation reference signal (DM-RS) symbols for physical layer cell identity number 42. The time-dependent part of the DM-RS scrambling initialization is 0.

ncellid = 42;
ibar_SSB = 0;
dmrsSym = nrPBCHDMRS(ncellid,ibar_SSB);

% Obtain resource element indices for the PBCH DM-RS.

dmrsInd = nrPBCHDMRSIndices(ncellid);

% Create a resource grid containing the generated DM-RS symbols.

nTxAnts = 1;
txGrid = complex(zeros([240 14 nTxAnts]));
txGrid(dmrsInd) = dmrsSym;

% Modulate the resource grid using the specified FFT length and cyclic prefix length.

nFFT = 512;
cpLengths = ones(1,14) * 36;
cpLengths([1 8]) = 40;
nulls = [1:136 377:512].';
txWaveform = ofdmmod(txGrid,nFFT,cpLengths,nulls);

% Create a TDL-C channel model with the specified properties.

SR = 7.68e6;
channel = nrTDLChannel;
channel.NumReceiveAntennas = 1;
channel.SampleRate = SR;
channel.DelayProfile = 'TDL-C';
channel.DelaySpread = 100e-9;
channel.MaximumDopplerShift = 20;

% Obtain the maximum number of delayed samples from the channel path by using the largest delay and the implementation delay of the channel filter.

chInfo = info(channel);
maxChDelay = ceil(max(chInfo.PathDelays*SR)) + chInfo.ChannelFilterDelay;

% To flush delayed samples from the channel, append zeros at the end of the transmitted waveform corresponding to the maximum number of delayed samples and the number of transmit antennas. Transmit the padded waveform through the TDL-C channel model.

[rxWaveform,pathGains] = channel([txWaveform; zeros(maxChDelay,nTxAnts)]);

% Estimate timing offset for the transmission using the DM-RS symbols as reference symbols. The OFDM modulation of the reference symbols spans 20 resource blocks at 15 kHz subcarrier spacing and uses initial slot number 0.

nrb = 20;
scs = 15;
initialSlot = 0;
offset = nrTimingEstimate(rxWaveform,nrb,scs,initialSlot,dmrsInd,dmrsSym);

% Synchronize the received waveform according to the estimated timing offset.

rxWaveform = rxWaveform(1+offset:end,:);

%  Create a received resource grid containing the demodulated and synchronized received waveform.

rxLength = sum(cpLengths) + nFFT*numel(cpLengths);
cpFraction = 0.55;
symOffsets = fix(cpLengths * cpFraction);
rxGrid = ofdmdemod(rxWaveform(1:rxLength,:),nFFT,cpLengths,symOffsets,nulls);

%  Obtain the practical channel estimate.

H = nrChannelEstimate(rxGrid,dmrsInd,dmrsSym);

% Obtain the perfect channel estimate.

pathFilters = getPathFilters(channel);
H_ideal = nrPerfectChannelEstimate(pathGains,pathFilters,nrb,scs,initialSlot,offset);

%  Compare practical and perfect channel estimates.

figure;
subplot(1,2,1);
imagesc(abs(H));
xlabel('OFDM Symbol');
ylabel('Subcarrier');
title('Practical Estimate Magnitude');
subplot(1,2,2);
imagesc(abs(H_ideal));
xlabel('OFDM Symbol');
ylabel('Subcarrier');
title('Perfect Estimate Magnitude');

% Copyright 2019 The MathWorks, Inc.