Improved noise and SNR estimation

master
ismagom 10 years ago
parent 4429d45761
commit d81f8dfbcd

@ -117,13 +117,13 @@ int chest_dl_init(chest_dl_t *q, lte_cell_t cell)
}
/* Set default time/freq filters */
//float f[3]={0.1, 0.8, 0.1};
//chest_dl_set_filter_freq(q, f, 3);
float f[3]={0.1, 0.8, 0.1};
chest_dl_set_filter_freq(q, f, 3);
float f[5]={0.05, 0.15, 0.6, 0.15, 0.05};
chest_dl_set_filter_freq(q, f, 5);
//float f[5]={0.05, 0.15, 0.6, 0.15, 0.05};
//chest_dl_set_filter_freq(q, f, 5);
float t[2]={0.15, 0.85};
float t[2]={0.1, 0.9};
chest_dl_set_filter_time(q, t, 2);
q->cell = cell;
@ -247,6 +247,10 @@ static void average_pilots(chest_dl_t *q, uint32_t port_id)
memcpy(&pilot_tmp(0), &pilot_est(0), nref * sizeof(cf_t));
}
}
#ifdef NOISE_POWER_USE_ESTIMATES
q->noise_estimate[port_id] = estimate_noise_port(q, port_id, q->tmp_freqavg);
#endif
for (l=0;l<refsignal_cs_nof_symbols(port_id);l++) {
/* Filter in time domain. */
@ -268,9 +272,6 @@ static void average_pilots(chest_dl_t *q, uint32_t port_id)
memcpy(&pilot_avg(0), &pilot_tmp(0), nref * sizeof(cf_t));
}
}
#ifdef NOISE_POWER_USE_ESTIMATES
q->noise_estimate[port_id] = estimate_noise_port(q, port_id, q->pilot_estimates_average[port_id]);
#endif
}
@ -345,20 +346,20 @@ int chest_dl_estimate_port(chest_dl_t *q, cf_t *input, cf_t *ce, uint32_t sf_idx
/* Get references from the input signal */
refsignal_cs_get_sf(q->cell, port_id, input, q->pilot_recv_signal[port_id]);
/* Compute RSRP for the references in this port */
q->rsrp[port_id] = chest_dl_rsrp(q, port_id);
if (port_id == 0) {
/* compute rssi only for port 0 */
q->rssi[port_id] = chest_dl_rssi(q, input, port_id);
}
/* Use the known CSR signal to compute Least-squares estimates */
vec_prod_conj_ccc(q->pilot_recv_signal[port_id], q->csr_signal.pilots[port_id/2][sf_idx],
q->pilot_estimates[port_id], REFSIGNAL_NUM_SF(q->cell.nof_prb, port_id));
/* Average pilot estimates */
average_pilots(q, port_id);
/* Compute RSRP for the channel estimates in this port */
q->rsrp[port_id] = chest_dl_rsrp(q, port_id);
if (port_id == 0) {
/* compute rssi only for port 0 */
q->rssi[port_id] = chest_dl_rssi(q, input, port_id);
}
/* Interpolate to create channel estimates for all resource grid */
if (ce != NULL) {
interpolate_pilots(q, ce, port_id);
@ -381,16 +382,18 @@ int chest_dl_estimate(chest_dl_t *q, cf_t *input, cf_t *ce[MAX_PORTS], uint32_t
}
float chest_dl_get_noise_estimate(chest_dl_t *q) {
return vec_acc_ff(q->noise_estimate, q->cell.nof_ports)/q->cell.nof_ports;
float noise = vec_acc_ff(q->noise_estimate, q->cell.nof_ports)/q->cell.nof_ports;
#ifdef NOISE_POWER_USE_ESTIMATES
return noise*sqrtf(lte_symbol_sz(q->cell.nof_prb));
#else
return noise;
#endif
}
float chest_dl_get_snr(chest_dl_t *q) {
// Uses RSRP as an estimation of the useful signal power
#ifdef NOISE_POWER_USE_ESTIMATES
return chest_dl_get_rsrp(q)/chest_dl_get_noise_estimate(q)/sqrt(2*lte_symbol_sz(q->cell.nof_prb));
#else
return chest_dl_get_rsrp(q)/chest_dl_get_noise_estimate(q)/sqrt(2);
#endif
return chest_dl_get_rsrp(q)/chest_dl_get_noise_estimate(q)/sqrt(2)/q->cell.nof_ports;
}
float chest_dl_get_rssi(chest_dl_t *q) {
@ -406,10 +409,8 @@ float chest_dl_get_rsrq(chest_dl_t *q) {
}
float chest_dl_get_rsrp(chest_dl_t *q) {
// return linear average from port 0 only
//return q->rsrp[0];
// return linear average from all ports
return vec_acc_ff(q->rsrp, q->cell.nof_ports)/q->cell.nof_ports;
// return sum of power received from all tx ports
return vec_acc_ff(q->rsrp, q->cell.nof_ports);
}

@ -195,11 +195,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
sf_idx = sf%10;
}
/* Loop through the 10 subframes */
if (chest_dl_estimate(&chest, input_signal, ce, sf_idx)) {
mexErrMsgTxt("Error running channel estimator\n");
return;
}
if (cell.nof_ports == 1) {
predecoding_single(&cheq, input_signal, ce[0], output_signal2, nof_re, chest_dl_get_noise_estimate(&chest));
} else {
@ -244,8 +244,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
}
if (nlhs >= 4) {
plhs[3] = mxCreateDoubleScalar(chest_dl_get_noise_estimate(&chest));
plhs[3] = mxCreateDoubleScalar(chest_dl_get_snr(&chest));
}
if (nlhs >= 5) {
plhs[4] = mxCreateDoubleScalar(chest_dl_get_noise_estimate(&chest));
}
if (nlhs >= 6) {
plhs[5] = mxCreateDoubleScalar(chest_dl_get_rsrp(&chest));
}
chest_dl_free(&chest);
precoding_free(&cheq);
return;
}

@ -55,7 +55,7 @@ void prb_cp_ref(cf_t **input, cf_t **output, int offset, int nof_refs,
int ref_interval = ((RE_X_RB / nof_refs) - 1);
memcpy(*output, *input, offset * sizeof(cf_t));
print_indexes(*output, offset);
print_indexes(*input, offset);
*input += offset;
*output += offset;
for (i = 0; i < nof_intervals - 1; i++) {
@ -65,7 +65,7 @@ void prb_cp_ref(cf_t **input, cf_t **output, int offset, int nof_refs,
(*input)++;
}
memcpy(*output, *input, ref_interval * sizeof(cf_t));
print_indexes(*output, ref_interval);
print_indexes(*input, ref_interval);
*output += ref_interval;
*input += ref_interval;
}
@ -76,7 +76,7 @@ void prb_cp_ref(cf_t **input, cf_t **output, int offset, int nof_refs,
(*input)++;
}
memcpy(*output, *input, (ref_interval - offset) * sizeof(cf_t));
print_indexes(*output, ref_interval-offset);
print_indexes(*input, ref_interval-offset);
*output += (ref_interval - offset);
*input += (ref_interval - offset);
}
@ -84,7 +84,7 @@ void prb_cp_ref(cf_t **input, cf_t **output, int offset, int nof_refs,
void prb_cp(cf_t **input, cf_t **output, int nof_prb) {
memcpy(*output, *input, sizeof(cf_t) * RE_X_RB * nof_prb);
print_indexes(*output, RE_X_RB);
print_indexes(*input, RE_X_RB);
*input += nof_prb * RE_X_RB;
*output += nof_prb * RE_X_RB;
}
@ -92,7 +92,7 @@ void prb_cp(cf_t **input, cf_t **output, int nof_prb) {
void prb_cp_half(cf_t **input, cf_t **output, int nof_prb) {
memcpy(*output, *input, sizeof(cf_t) * RE_X_RB * nof_prb / 2);
print_indexes(*output, RE_X_RB/2);
print_indexes(*input, RE_X_RB/2);
*input += nof_prb * RE_X_RB / 2;
*output += nof_prb * RE_X_RB / 2;
}

@ -240,7 +240,7 @@ int main(int argc, char **argv) {
printf("Could not decode PCFICH\n");
exit(-1);
} else {
if (cfi_corr > 3 && cfi == 1) {
if (cfi_corr > 2.8 && cfi == 1) {
exit(0);
} else {
exit(-1);

@ -44,6 +44,8 @@ void help()
("[decoded_ok, llr, rm, bits, symbols] = liblte_pdsch(enbConfig, pdschConfig, trblklen, rxWaveform)\n\n");
}
//extern int indices[2048];
/* the gateway function */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
@ -61,7 +63,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
uint32_t rv;
uint32_t rnti32;
if (nrhs != NOF_INPUTS) {
if (nrhs < NOF_INPUTS) {
help();
return;
}
@ -150,7 +152,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
prb_alloc.slot[0].prb_idx[i] = false;
for (int j=0;j<prb_alloc.slot[0].nof_prb && !prb_alloc.slot[0].prb_idx[i];j++) {
if ((int) prbset[j] == i) {
prb_alloc.slot[0].prb_idx[i] = true;
prb_alloc.slot[0].prb_idx[i] = true;
}
}
}
@ -193,7 +195,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
} else {
noise_power = chest_dl_get_noise_estimate(&chest);
}
if (pdsch_harq_setup(&harq_process, mcs, &prb_alloc)) {
mexErrMsgTxt("Error configuring HARQ process\n");
return;
@ -219,6 +221,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
if (nlhs >= 4) {
mexutils_write_cf(pdsch.pdsch_d, &plhs[3], harq_process.prb_alloc.re_sf[sf_idx], 1);
}
if (nlhs >= 5) {
mexutils_write_cf(ce[0], &plhs[4], 12*14*cell.nof_prb, 1);
}
chest_dl_free(&chest);
lte_fft_free(&fft);

@ -4,8 +4,10 @@
clear
plot_noise_estimation_only=true;
SNR_values_db=linspace(0,30,8);
Nrealizations=1;
Nrealizations=1 ;
preEVM = zeros(length(SNR_values_db),Nrealizations);
postEVM_mmse = zeros(length(SNR_values_db),Nrealizations);
@ -25,7 +27,7 @@ rng(1); % Configure random number generators
cfg.Seed = 2; % Random channel seed
cfg.NRxAnts = 2; % 1 receive antenna
cfg.DelayProfile = 'EVA'; % EVA delay spread
cfg.DopplerFreq = 5; % 120Hz Doppler frequency
cfg.DopplerFreq = 120; % 120Hz Doppler frequency
cfg.MIMOCorrelation = 'Low'; % Low (no) MIMO correlation
cfg.InitTime = 0; % Initialize at time zero
cfg.NTerms = 16; % Oscillators used in fading model
@ -53,6 +55,7 @@ P = gridsize(3); % Number of transmit antenna ports
for nreal=1:Nrealizations
%% Transmit Resource Grid
txGrid = [];
RSRP = [];
%% Payload Data Generation
% Number of bits needed is size of resource grid (K*L*P) * number of bits
@ -77,7 +80,7 @@ for sf = 0:10
subframe = lteDLResourceGrid(enb);
% Map input symbols to grid
subframe(:) = inputSym;
%subframe(:) = inputSym;
% Generate synchronizing signals
pssSym = ltePSS(enb);
@ -121,12 +124,10 @@ cfg.SamplingRate = info.SamplingRate;
% Pass data through the fading channel model
%rxWaveform = lteFadingChannel(cfg,txWaveform);
rxWaveform = txWaveform;
rxWaveform = txWaveform;
%% Additive Noise
channel_gain = mean(rxWaveform(:).*conj(rxWaveform(:)))/mean(txWaveform(:).*conj(txWaveform(:)));
% Calculate noise gain
N0 = 1/(sqrt(2.0*enb.CellRefP*double(info.Nfft))*SNR);
@ -144,68 +145,79 @@ rxWaveform = rxWaveform(1+offset:end,:);
%% OFDM Demodulation
rxGrid = lteOFDMDemodulate(enb,rxWaveform);
rxGrid = rxGrid(:,1:140);
addpath('../../debug/lte/phy/lib/ch_estimation/test')
%% Channel Estimation
[estChannel, noiseEst(snr_idx)] = lteDLChannelEstimate(enb,cec,rxGrid);
output=[];
snrest = zeros(10,1);
nulls = rxGrid([1:5 68:72],6:7);
noiseEst(snr_idx) = var(nulls(:));
%for i=0:9
i=0;
% if (SNR_values_db(snr_idx) < 25)
[d, a, out, snrest(i+1)] = liblte_chest(enb.NCellID,enb.CellRefP,rxGrid(:,i*14+1:(i+1)*14),[0.1 0.8 0.1],[0.1 0.9],i);
% else
% [d, a, out, snrest(i+1)] = liblte_chest(enb.NCellID,enb.CellRefP,rxGrid(:,i*14+1:(i+1)*14),[0.05 0.9 0.05],[],i);
% end
snrest_liblte = zeros(10,1);
noise_liblte = zeros(10,1);
rsrp = zeros(1,10);
for i=0:9
[d, a, out, snrest_liblte(i+1), noise_liblte(i+1), rsrp(i+1)] = liblte_chest(enb.NCellID,enb.CellRefP,rxGrid(:,i*14+1:(i+1)*14),[0.1 0.8 0.1],[0.1 0.9],i);
output = [output out];
%end
SNRest(snr_idx)=mean(snrest);
disp(10*log10(SNRest(snr_idx)))
end
RSRP = [RSRP rsrp];
meanRSRP(snr_idx)=mean(rsrp);
SNR_liblte(snr_idx)=mean(snrest_liblte);
noiseEst_liblte(snr_idx)=mean(noise_liblte);
if ~plot_noise_estimation_only
%% MMSE Equalization
% eqGrid_mmse = lteEqualizeMMSE(rxGrid, estChannel, noiseEst(snr_idx));
%
% eqGrid_liblte = reshape(output,size(eqGrid_mmse));
%
% % Analysis
%
% %Compute EVM across all input values EVM of pre-equalized receive signal
% preEqualisedEVM = lteEVM(txGrid,rxGrid);
% fprintf('%d-%d: Pre-EQ: %0.3f%%\n', ...
% snr_idx,nreal,preEqualisedEVM.RMS*100);
%
%
% %EVM of post-equalized receive signal
% postEqualisedEVM_mmse = lteEVM(txGrid,reshape(eqGrid_mmse,size(txGrid)));
% fprintf('%d-%d: MMSE: %0.3f%%\n', ...
% snr_idx,nreal,postEqualisedEVM_mmse.RMS*100);
%
% postEqualisedEVM_liblte = lteEVM(txGrid,reshape(eqGrid_liblte,size(txGrid)));
% fprintf('%d-%d: liblte: %0.3f%%\n', ...
% snr_idx,nreal,postEqualisedEVM_liblte.RMS*100);
%
% preEVM(snr_idx,nreal) = preEqualisedEVM.RMS;
% postEVM_mmse(snr_idx,nreal) = mean([postEqualisedEVM_mmse.RMS]);
% postEVM_liblte(snr_idx,nreal) = mean([postEqualisedEVM_liblte.RMS]);
eqGrid_mmse = lteEqualizeMMSE(rxGrid, estChannel, noiseEst(snr_idx));
eqGrid_liblte = reshape(output,size(eqGrid_mmse));
% Analysis
%Compute EVM across all input values EVM of pre-equalized receive signal
preEqualisedEVM = lteEVM(txGrid,rxGrid);
fprintf('%d-%d: Pre-EQ: %0.3f%%\n', ...
snr_idx,nreal,preEqualisedEVM.RMS*100);
%EVM of post-equalized receive signal
postEqualisedEVM_mmse = lteEVM(txGrid,reshape(eqGrid_mmse,size(txGrid)));
fprintf('%d-%d: MMSE: %0.3f%%\n', ...
snr_idx,nreal,postEqualisedEVM_mmse.RMS*100);
postEqualisedEVM_liblte = lteEVM(txGrid,reshape(eqGrid_liblte,size(txGrid)));
fprintf('%d-%d: liblte: %0.3f%%\n', ...
snr_idx,nreal,postEqualisedEVM_liblte.RMS*100);
preEVM(snr_idx,nreal) = preEqualisedEVM.RMS;
postEVM_mmse(snr_idx,nreal) = mean([postEqualisedEVM_mmse.RMS]);
postEVM_liblte(snr_idx,nreal) = mean([postEqualisedEVM_liblte.RMS]);
end
end
end
% subplot(1,2,1)
% plot(SNR_values_db, mean(preEVM,2), ...
% SNR_values_db, mean(postEVM_mmse,2), ...
% SNR_values_db, mean(postEVM_liblte,2))
% legend('No Eq','MMSE-lin','MMSE-liblte')
% grid on
%
if ~plot_noise_estimation_only
plot(SNR_values_db, mean(preEVM,2), ...
SNR_values_db, mean(postEVM_mmse,2), ...
SNR_values_db, mean(postEVM_liblte,2))
legend('No Eq','MMSE-lin','MMSE-liblte')
grid on
end
% subplot(1,2,2)
%SNR_liblte = 1./(SNRest*sqrt(2.0*enb.CellRefP*double(info.Nfft)));
SNR_liblte = 1./(SNRest*sqrt(2.0));
SNR_matlab = 1./(noiseEst*sqrt(2.0));
if plot_noise_estimation_only
SNR_matlab = 1./(noiseEst*sqrt(2.0)*enb.CellRefP);
subplot(1,3,1)
plot(SNR_values_db, SNR_values_db, SNR_values_db, 10*log10(SNR_liblte),SNR_values_db, 10*log10(SNR_matlab))
legend('Theory','libLTE','Matlab')
plot(SNR_values_db, SNR_values_db, SNR_values_db, 10*log10(SNR_liblte),SNR_values_db, 10*log10(SNR_matlab))
%plot(SNR_values_db, 10*log10(noiseTx), SNR_values_db, 10*log10(SNRest),SNR_values_db, 10*log10(noiseEst))
legend('Theory','libLTE','Matlab')
subplot(1,3,2)
plot(SNR_values_db, 10*log10(noiseTx), SNR_values_db, 10*log10(noiseEst_liblte),SNR_values_db, 10*log10(noiseEst))
legend('Theory','libLTE','Matlab')
subplot(1,3,3)
plot(1:10*length(SNR_values_db),RSRP,10*(1:length(SNR_values_db)),meanRSRP)
end

@ -4,25 +4,29 @@
% A structure |enbConfig| is used to configure the eNodeB.
%clear
recordedSignal=[];
Npackets = 4;
SNR_values = 1;%linspace(2,6,4);
SNR_values = linspace(5,6,4);
%% Choose RMC
[waveform,rgrid,rmccFgOut] = lteRMCDLTool('R.11',[1;0;0;1]);
[waveform,rgrid,rmccFgOut] = lteRMCDLTool('R.4',[1;0;0;1]);
waveform = sum(waveform,2);
rmccFgOut = struct('NCellID',1,'CellRefP',1,'CFI',1,'NDLRB',15,'SamplingRate',3.84e6,'Nfft',256,'DuplexMode','FDD','CyclicPrefix','Normal');
rmccFgOut.PDSCH.RNTI = 1234;
rmccFgOut.PDSCH.PRBSet = repmat(transpose(0:rmccFgOut.NDLRB-1),1,2);
rmccFgOut.PDSCH.TxScheme = 'Port0';
rmccFgOut.PDSCH.NLayers = 1;
rmccFgOut.PDSCH.NTurboDecIts = 5;
rmccFgOut.PDSCH.Modulation = {'64QAM'};
rmccFgOut.PDSCH.TrBlkSizes = [0 5992*ones(1,4) 0 5992*ones(1,4)];
rmccFgOut.PDSCH.RV = 0;
if ~isempty(recordedSignal)
rmccFgOut = struct('NCellID',1,'CellRefP',1,'CFI',1,'NDLRB',15,'SamplingRate',3.84e6,'Nfft',256,'DuplexMode','FDD','CyclicPrefix','Normal');
rmccFgOut.PDSCH.RNTI = 1234;
rmccFgOut.PDSCH.PRBSet = repmat(transpose(0:rmccFgOut.NDLRB-1),1,2);
rmccFgOut.PDSCH.TxScheme = 'Port0';
rmccFgOut.PDSCH.NLayers = 1;
rmccFgOut.PDSCH.NTurboDecIts = 5;
rmccFgOut.PDSCH.Modulation = {'64QAM'};
rmccFgOut.PDSCH.TrBlkSizes = [0 5992*ones(1,4) 0 5992*ones(1,4)];
rmccFgOut.PDSCH.RV = 0;
end
flen=rmccFgOut.SamplingRate/1000;
Nsf = 9;
%% Setup Fading channel model
@ -58,15 +62,19 @@ for snr_idx=1:length(SNR_values)
N0 = 1/(sqrt(2.0*rmccFgOut.CellRefP*double(rmccFgOut.Nfft))*SNR);
for i=1:Npackets
%% Fading
rxWaveform = lteFadingChannel(cfg,waveform);
%% Noise Addition
noise = N0*complex(randn(size(rxWaveform)), randn(size(rxWaveform))); % Generate noise
rxWaveform = rxWaveform + noise;
if isempty(recordedSignal)
%% Fading
rxWaveform = lteFadingChannel(cfg,waveform);
%rxWaveform = waveform;
%% Noise Addition
noise = N0*complex(randn(size(rxWaveform)), randn(size(rxWaveform))); % Generate noise
rxWaveform = rxWaveform + noise;
else
rxWaveform = recordedSignal;
end
rxWaveform = x;
%% Demodulate
frame_rx = lteOFDMDemodulate(rmccFgOut, rxWaveform);
@ -79,7 +87,7 @@ for snr_idx=1:length(SNR_values)
% Perform channel estimation
[hest, nest] = lteDLChannelEstimate(rmccFgOut, cec, subframe_rx);
[cws,symbols] = ltePDSCHDecode(rmccFgOut,rmccFgOut.PDSCH,subframe_rx,hest,nest);
[cws,symbols,indices,pdschSymbols,pdschHest] = ltePDSCHDecode2(rmccFgOut,rmccFgOut.PDSCH,subframe_rx,hest,nest);
[trblkout,blkcrc] = lteDLSCHDecode(rmccFgOut,rmccFgOut.PDSCH, ...
rmccFgOut.PDSCH.TrBlkSizes(sf_idx+1),cws);
@ -88,7 +96,7 @@ for snr_idx=1:length(SNR_values)
%% Same with libLTE
if (rmccFgOut.PDSCH.TrBlkSizes(sf_idx+1) > 0)
[dec2, llr, pdschRx, pdschSymbols2] = liblte_pdsch(rmccFgOut, rmccFgOut.PDSCH, ...
[dec2, data, pdschRx, pdschSymbols2, deb] = liblte_pdsch(rmccFgOut, rmccFgOut.PDSCH, ...
rmccFgOut.PDSCH.TrBlkSizes(sf_idx+1), ...
subframe_waveform);
else
@ -96,15 +104,17 @@ for snr_idx=1:length(SNR_values)
end
decoded_liblte(snr_idx) = decoded_liblte(snr_idx)+dec2;
end
x = x(flen*10+1:end);
if ~isempty(recordedSignal)
recordedSignal = recordedSignal(flen*10+1:end);
end
end
fprintf('SNR: %.1f\n',SNRdB)
fprintf('SNR: %.1f. Decoded: %d-%d\n',SNRdB, decoded(snr_idx), decoded_liblte(snr_idx))
end
if (length(SNR_values)>1)
semilogy(SNR_values,1-decoded/Npackets/(Nsf+1),'bo-',...
SNR_values,1-decoded_liblte/Npackets/(Nsf+1), 'ro-')
semilogy(SNR_values,1-decoded/Npackets/(Nsf),'bo-',...
SNR_values,1-decoded_liblte/Npackets/(Nsf), 'ro-')
grid on;
legend('Matlab','libLTE')
xlabel('SNR (dB)')

@ -2,7 +2,7 @@
clear
blen=1008;
SNR_values_db=linspace(-1,0.5,6);
Nrealizations=10000;
Nrealizations=1000;
addpath('../../debug/lte/phy/lib/fec/test')

Loading…
Cancel
Save