# How to graph BER with perfect results?

1. ## How to graph BER with perfect results?

Hello

i have a problem with matlab. suppose after completing a simulation with this results :
snr=0:5:35
for qpsk i have this vector: ber_qpsk=[1e-3 1e-4 1e-5 0 0 0 0]
for 16-qam i have this vector: ber_16qam=[1e-2 1e-3 1e-4 1e-5 0 0 0]
for 64-qam i have this vector: ber_64qam=[0.5e-2 1e-2 1e-3 1e-4 1e-5 0 0]

now, i want to graph them all in one graph. tries berfit and i got '[]'.

•

2. ## Re: How to graph BER with perfect results?

anyone can help me with that ?

3. ## Re: How to graph BER with perfect results?

Originally Posted by WUID
anyone can help me with that ?
It seems that there is some problem in simulation as you are getting 0 BER just after 1e-5. Just plot them on one graph using plot but it wont be a water fall curve.

4. ## Re: How to graph BER with perfect results?

berfit is doing more than you want to do. To plot those three vectors in one graph just use:

figure(1);
subplot(311), plot(snr, ber_qpsk);
subplot(312), plot(snr, ber_16qam);
subplot(313), plot(snr, ber_64qam);

Then, you can use the bertool command to create vectors of the theoretical BER curves for psk and qam, and plot them next to those.

- Ryan

1 members found this post helpful.

5. ## Re: How to graph BER with perfect results?

Hi,

Try this

snr=0:5:30;
ber_qpsk=[1e-3 1e-4 1e-5 0 0 0 0];
ber_16qam=[1e-2 1e-3 1e-4 1e-5 0 0 0];
ber_64qam=[0.5e-2 1e-2 1e-3 1e-4 1e-5 0 0];
figure(1);
semilogy(snr,ber_qpsk,'k-o','LineWidth',1);
hold on
semilogy(snr,ber_16qam,'k-x','LineWidth',1);
hold on
semilogy(snr,ber_16qam,'k-d','LineWidth',1);
grid on
legend('qpsk','16qam','64qam',3);
xlabel('SNR (dB)');
ylabel('BER');
title(['BER versus SNR for different receivers'])
axis([min(snr) max(snr) 10^(-6) 1]);

Best regards

2 members found this post helpful.

6. ## Re: How to graph BER with perfect results?

hi , i'll try this solutions thx.
but in addition i asked my lecturer and he told that i need to increase the number of bits. so for my expected value for my simulation i want to be around ber = 10^-6 so this demanding to transmit at least 10^7 bits. i can do it, but this can take me a day i've checked it. i know that my sim' isn't so efficient , cause i have many FOR's.

so here comes your part :)

i can reduce the for loops in some cases like :
for 1:N_sym %number of symbols usually 1000 up to 4500
x_Tx_bits=randsrc(1,2256,[0 1]);
CodedBits=encode(henc,x_Tx_Bits); %using ldpc encoder with 2256*4512 matrix
end

i want create matrix of bits:

x_Tx_Bits=randsrc(N_sym,2256,[0 1]);

is there any way i can encode this giant matrix at once ?

•

7. ## Re: How to graph BER with perfect results?

Hi,

Yes to simulate ber=10^-6, you must send 10^7 bits but the incertitude is big. If you sen 10^8 you will have ber=10^-6 with 1% incertitude.

You can encode a matrix. This is an example :

N_sym=4512;
n = 6; k = 4; % Set codeword length and message length.
x_Tx_Bits=randsrc(N_sym,4,[0 1]); % Message is a binary matrix.
code = encode(x_Tx_Bits,n,k,'cyclic');

But, you don't give enough parameters of LDPC coder. henc ?

1 members found this post helpful.

8. ## Re: How to graph BER with perfect results?

I am posting my raw code. meaning i am doing all my tests on this code. If i am getting expected results I am changing my simulation files.

The code including awgn/LTI/LTV with all the blocks that needed from ofdm scheme(not including Interleaving , * exist but not straightforward )
The LDPC matrix taken from the net and it working very well for AWGN+LTI channel I'm getting with QPSK the BER=0 (i know it's wrong write it) issue just with SNR~1.5dB.
My goal is to shorten the LTV stage with OVERLAP&ADD method, but first i want to stabilize my results.

For now the LTI+awgn working well just need to fix BER.
In the LTV(section) channel i have hard to figure the math and to normalize my parameters good enough.

many thanks to u all, it is very appreciated . keep helping me :)

Code:
```function NumErrors=OurLDPCEncandDec_without_MMSE(itertion,Rate,Modulation,SNRdB)
% close all; clc ; clear
% Rate=2;Modulation=4;SNRdB=3

% Rate: 2 for 1/2, 4 for 3/4, 6 for 5/6, 8 for 7/8
NumErrors=zeros(1,itertion);
for iter=1:itertion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters Section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

H1=H(:,1:2256);
H2=H(:,2257:4512);
H=[H2 H1];

if Rate==2
ParityLength=2256;
elseif Rate==4
ParityLength=2256/3;
elseif Rate==6
ParityLength=464;
elseif Rate==8
ParityLength=320;
elseif Rate==10
ParityLength=240;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create Signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xt_long_with_CP=[];
xt_long_bits=[];
N_sym=100;
henc = fec.ldpcenc(H);
for k=1:N_sym

InfoBits=randsrc(1,2256,[1 0]); % Generate 2256 information bits
xt_long_bits=[xt_long_bits InfoBits];
% Encode and Modulate
CodedBits= encode(henc,InfoBits);
% Rate Matching
ParityLocations=2256+randperm(2256);
ParityLocations=ParityLocations(1:ParityLength);
TotalLocations(:,k)=[1:2256,ParityLocations]; % Should be an output of the encoder
CodedBits=CodedBits(TotalLocations(:,k));

Symbols=ModulateQAM(CodedBits,Modulation);

FFTsize=2^(ceil(log2(length(Symbols))));
ActiveSCs=length(Symbols);
PilotLocations=1:4:FFTsize; %pilot loc /for mod=4,16 spacing = 3. /for mod=64 spacing =4

SC_loc=setdiff(1:FFTsize,PilotLocations); SC_loc=SC_loc(1:ActiveSCs); %info location

%%%%pilots%%%%
x_f=zeros(FFTsize,1); x_f(SC_loc)=Symbols;
Xp=randsrc(1,length(PilotLocations),[-1 1]);    %Pilot generation
Pilots_PRBS(:,k)=Xp.'; % PRBS patterns
x_f(PilotLocations)=Xp.'; %insertion of the pilots
xt=ifft(x_f);
xt_withCP=[xt(end-1023:end);xt];
xt_long_with_CP=[xt_long_with_CP;xt_withCP];
end
%%%%%%%%%%%%%%%%%%%%%%%%Channel%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs=10e6; %Sampling rate
N_samples=length(xt_long_with_CP);
T=N_samples/Fs; %Sample time
SNR=10^(SNRdB/10);

% Taking channel parameters from ITU table
TapPowerdB=[0 -9.7 -19.2 -22.8];   % in dB
Delay=[1 110 190 410]*1e-9; % in Sec
%     ChannelLength=length(TapPower);

a_samples=[];
%%%%Tap power in dB to linear%%%%
TapPower=10.^(TapPowerdB/10);
%         TapPower=TapPower/norm(TapPower); % Normalize to unit channel power

F_low=f_m*10; % usaully 10*f_m the low sampling
% freq of the channel coefficients

N_low=ceil(T*F_low);

%%%%Gausssian random process pass throw filter(flat/U-shaped)%%%%

h=firpm(64,[0 2*f_m/F_low 2*f_m/F_low*1.2 1],[1 1 0 0 ]);
h=h/norm(h);
%        channel_length=length(h);
for k=1:length(TapPower)

a_low_temp=conv((randn(1,N_low+64)+1i*randn(1,N_low+64))/sqrt(2),h);%was ...,h)*TapPower(k);
a_low_temp=a_low_temp(65:N_low+64); % Generate first tap sampled at low freq
a_low(k,:)=a_low_temp;
channel_length=6;
%%%%Upsample To Fs%%%%

a_samples(k,:)=resample(a_low(k,:),Fs,F_low);
Sigvar=var(a_samples(k,:));
a_samples(k,:)=a_samples(k,:)*sqrt((TapPower(k)/Sigvar));
end
%%%%Initilaize Delay Vector%%%%
DelayInSamples=ceil(Delay*Fs);
DelayInSamples=DelayInSamples-DelayInSamples(1);

y=zeros(1,length(xt_long_with_CP));

%%%%Convolution%%%%
for k=max(DelayInSamples)+1:length(xt_long_with_CP)
for j=1:length(TapPower)
a(j)=xt_long_with_CP(k-DelayInSamples(j))*a_samples(j,k);
end
y(k)=sum(a);

end

SNR_corrected=SNR*(1366+2256)/FFTsize; %need to generic fit for qpsk 2256 infobits + 1366 pilots
sigma_noise=norm(xt_long_with_CP)/sqrt(length(y)*SNR_corrected);
y=y+sigma_noise*(randn(size(y))+1i*randn(size(y)))/sqrt(2);

yBeforeCP_RemoveMatrix=reshape(y.',FFTsize+1024,N_sym);
yCP_RemoveMatrix=yBeforeCP_RemoveMatrix(1025:end,:);

Y_matrix_PostFFT=fft(yCP_RemoveMatrix);

for kk=1:N_sym,
%     H_est = MMSE_CE(Y_matrix_PostFFT(:,kk),Pilots_PRBS(:,kk),PilotLocations,FFTsize,4,h,SNR);
%     h_est = ifft(H_est); h_DFT = h_est(1:channel_length);
%     H_DFT = fft(h_DFT,FFTsize);
Y_Pilots=Y_matrix_PostFFT(PilotLocations,kk);

EstTimeDomainChannel=ifft(Y_Pilots.*Pilots_PRBS(:,kk));%Taking only pilots
EstTimeDomainChannel=EstTimeDomainChannel(1:(channel_length)); % Kill contibutions outside the CP duration
H_est=fft(EstTimeDomainChannel,FFTsize); % est. channel

demodulatedSCs=Y_matrix_PostFFT(:,kk)./H_est;

Y_equalized_withoutPilots(:,kk)=demodulatedSCs(SC_loc);
ppSNR(:,kk)=abs(H_est(SC_loc)).^2/(sigma_noise^2*FFTsize);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Decode Signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DecodedBits_long=[];
for kk=1:N_sym,
LLRs=SISO_LLR(Y_equalized_withoutPilots(:,kk),ppSNR(:,kk),Modulation);
LLR_rate_matched=zeros(1,2*2256);
LLR_rate_matched(TotalLocations(:,kk))=LLRs;

hdec=fec.ldpcdec(H);
DecodedBits=decode(hdec, LLR_rate_matched);
DecodedBits_long=[DecodedBits_long DecodedBits];
end
NumErrors(1,iter)=sum(DecodedBits_long ~=xt_long_bits); % Count errors
end

%-------------------------------------------------------

function QAM_Symbols=ModulateQAM(CodedBits,Modulation)

% QPSK mapping
MapQPSK= ([1+i 1-i -1+i -1-i])/sqrt(2);

% 16QAM mapping
Map_16QAM= ((2/sqrt(10))* ...
[ 0.5+0.5*i 0.5+1.5*i 0.5-0.5*i 0.5-1.5*i ...
1.5+0.5*i 1.5+1.5*i 1.5-0.5*i 1.5-1.5*i ....
-0.5+0.5*i -0.5+1.5*i -0.5-0.5*i -0.5-1.5*i ....
-1.5+0.5*i -1.5+1.5*i -1.5-0.5*i -1.5-1.5*i]);

% 64QAM mapping
Map_64QAM= ((2/sqrt(42))* ...
[ 1.5+1.5*i 1.5+0.5*i 1.5+2.5*i 1.5+3.5*i ...
1.5-1.5*i 1.5-0.5*i 1.5-2.5*i 1.5-3.5*i   ...
0.5+1.5*i 0.5+0.5*i 0.5+2.5*i  0.5+3.5*i  ...
0.5-1.5*i 0.5-0.5*i 0.5-2.5*i 0.5-3.5*i   ...
2.5+1.5*i 2.5+0.5*i 2.5+2.5*i  2.5+3.5*i  ...
2.5-1.5*i 2.5-0.5*i 2.5-2.5*i 2.5-3.5*i   ...
3.5+1.5*i 3.5+0.5*i 3.5+2.5*i  3.5+3.5*i  ...
3.5-1.5*i 3.5-0.5*i 3.5-2.5*i 3.5-3.5*i   ...
-1.5+1.5*i -1.5+0.5*i -1.5+2.5*i  -1.5+3.5*i  ...
-1.5-1.5*i -1.5-0.5*i -1.5-2.5*i -1.5-3.5*i   ...
-0.5+1.5*i -0.5+0.5*i -0.5+2.5*i  -0.5+3.5*i  ...
-0.5-1.5*i -0.5-0.5*i -0.5-2.5*i -0.5-3.5*i   ...
-2.5+1.5*i -2.5+0.5*i -2.5+2.5*i  -2.5+3.5*i  ...
-2.5-1.5*i -2.5-0.5*i -2.5-2.5*i -2.5-3.5*i   ...
-3.5+1.5*i -3.5+0.5*i -3.5+2.5*i  -3.5+3.5*i  ...
-3.5-1.5*i -3.5-0.5*i -3.5-2.5*i -3.5-3.5*i]);

if Modulation==4
Map=MapQPSK;
elseif Modulation==16
Map=Map_16QAM;
elseif Modulation==64
Map=Map_64QAM;

end;

BitsPerCarrier=log2(Modulation);
CodedBitMat=reshape(CodedBits.',BitsPerCarrier,length(CodedBits)/BitsPerCarrier);
SymbolIndices=bi2de(fliplr(CodedBitMat.'))+1;
QAM_Symbols=Map(SymbolIndices);```

---------- Post added at 22:50 ---------- Previous post was at 22:24 ----------

Originally Posted by Mohamed Chaker
Hi,

Yes to simulate ber=10^-6, you must send 10^7 bits but the incertitude is big. If you sen 10^8 you will have ber=10^-6 with 1% incertitude.

You can encode a matrix. This is an example :

N_sym=4512;
n = 6; k = 4; % Set codeword length and message length.
x_Tx_Bits=randsrc(N_sym,4,[0 1]); % Message is a binary matrix.
code = encode(x_Tx_Bits,n,k,'cyclic');

But, you don't give enough parameters of LDPC coder. henc ?
To simulate 10^8 bits i will need to seat for a week... i think... this is way i asking how can i get rid from most of the FOR loops.

With the code i sent last post i will have to generate (1e8/2256=44000 symbols=> for i=1:44000)(with QPSK code rate=0.5)

•

9. ## Re: How to graph BER with perfect results?

Hi,

You can began with optimizing your code and tested for 10^7 bits. After that simulate with 10^8 to get more precision for 10^-6 ber.

It is normal that BER simulation takes much time (more than one day) but it is possible to win time by changing code.

I have some question before i can propose some modifications

- Why you are using the iterative loop ? could be for iterative decoding of LDPC code ?!

- There a simple modification that you can do to save time like

. % Parameters Section must be before the iterative loop. it is useless to repeat loading 802_16e_4512_d5_H every time.

. hdec=fec.ldpcdec(H) instruction is repeated twice. Once at coding and one at decoding. Only one before the loop will be ok.

- encode(henc,InfoBits) doesn't support Matrix input. I will see what you can change after i analyse more the code.

Regards

Chaker

1 members found this post helpful.

10. ## Re: How to graph BER with perfect results?

Originally Posted by rhearty1
berfit is doing more than you want to do. To plot those three vectors in one graph just use:

figure(1);
subplot(311), plot(snr, ber_qpsk);
subplot(312), plot(snr, ber_16qam);
subplot(313), plot(snr, ber_64qam);

Then, you can use the bertool command to create vectors of the theoretical BER curves for psk and qam, and plot them next to those.

- Ryan
I must have been tired when I answered. This does not plot the three vectors in one graph. Mohamed's way is right.

Matlab is optimized for vector multiplication, so use the '.' operator as much as you can to avoid nested for loops. This is a godsend in communications simulation.

For instance if you want to implement a matched filter for your receiver, where received_sig and template_sig are the two signal vectors, you can multiply an integrate-and-dump in one line:

This would do the same as:
corr = 0;
for ii=1:LENGTH
end

In my experience this can really speed things up. I would like to quantify that statement, but I'm not on a computer with Matlab right now.

- Ryan

1 members found this post helpful.

11. ## Re: How to graph BER with perfect results?

Hi Mohamed
- Why you are using the iterative loop ? could be for iterative decoding of LDPC code ?!
I did the iterative loop because i got non consistent results with the LTV channel. e.g for 100 iter i getting numerror vector that most of it contain "zero" errors and part of it is with different number of errors. I believe is due to my misunderstanding the math completely and some of the vectors in the LTV section aren't normalized well.

- encode(henc,InfoBits) doesn't support Matrix input.
the encoder/decoder function receives 'henc' sparse matrix.

as for % Parameters Section your are very right, but this is the deal with my simulation(this code that is modified to function). I want to create Data matrix in one line:

InfoBits=randsrc(1000,2256,[0 1]);
and be able to encode it at once.

the aim right now with eliminating the FOR is in the "%Create Signal" part.

12. ## Re: How to graph BER with perfect results?

Hi,

I have made some change in your code by putting all fixed parameters before iteration loop and 1:N_sym loop. I was forced to add some parameters. Try the code if it give same results. I can't do this because i don't have SISO_LLR function.
function NumErrors=OurLDPCEncandDec_without_MMSE(itertion,R ate,Modulation,SNRdB)
% close all; clc ;
%clear all ;
% Rate=2;
% Modulation=4;
% SNRdB=3;
% itertion=5;
% Rate: 2 for 1/2, 4 for 3/4, 6 for 5/6, 8 for 7/8

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
% Parameters Section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
Symbols_length=2256;
CP_lenght=1024;
N_sym=100;

if Rate==2
ParityLength=2256;
elseif Rate==4
ParityLength=2256/3;
elseif Rate==6
ParityLength=464;
elseif Rate==8
ParityLength=320;
elseif Rate==10
ParityLength=240;
end;

H1=H(:,1:2256);
H2=H(:,2257:4512);
H=[H2 H1];
henc = fec.ldpcenc(H);
hdec = fec.ldpcdec(H);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
% IFFT/FFT Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
FFTsize=2^(ceil(log2(Symbols_length)));
ActiveSCs=Symbols_length;
PilotLocations=1:4:FFTsize; %pilot loc /for mod=4,16 spacing = 3. /for mod=64 spacing =4
SC_loc=setdiff(1:FFTsize,PilotLocations); SC_loc=SC_loc(1:ActiveSCs); %info location
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
%Channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
Fs=10e6; %Sampling rate
N_samples=N_sym*(FFTsize+CP_lenght);
T=N_samples/Fs; %Sample time

% Taking channel parameters from ITU table
TapPowerdB=[0 -9.7 -19.2 -22.8]; % in dB
Delay=[1 110 190 410]*1e-9; % in Sec
% ChannelLength=length(TapPower);

a_samples=[];
%%%%Tap power in dB to linear%%%%
TapPower=10.^(TapPowerdB/10);
% TapPower=TapPower/norm(TapPower); % Normalize to unit channel power

F_low=f_m*10; % usaully 10*f_m the low sampling
% freq of the channel coefficients
N_low=ceil(T*F_low);

%%%%Gausssian random process pass throw filter(flat/U-shaped)%%%%

h=firpm(64,[0 2*f_m/F_low 2*f_m/F_low*1.2 1],[1 1 0 0 ]);
h=h/norm(h);
% channel_length=length(h);
for k=1:length(TapPower)

a_low_temp=conv((randn(1,N_low+64)+1i*randn(1,N_lo w+64))/sqrt(2),h);%was ...,h)*TapPower(k);
a_low_temp=a_low_temp(65:N_low+64); % Generate first tap sampled at low freq
a_low(k,:)=a_low_temp;
channel_length=6;
%%%%Upsample To Fs%%%%

a_samples(k,:)=resample(a_low(k,:),Fs,F_low);
Sigvar=var(a_samples(k,:));
a_samples(k,:)=a_samples(k,:)*sqrt((TapPower(k)/Sigvar));
end

%%%%Initilaize Delay Vector%%%%
DelayInSamples=ceil(Delay*Fs);
DelayInSamples=DelayInSamples-DelayInSamples(1);

%%% Noise parameters
SNR=10^(SNRdB/10);
SNR_corrected=SNR*(1366+2256)/FFTsize; %need to generic fit for qpsk 2256 infobits + 1366 pilots

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
% Began iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%

NumErrors=zeros(1,itertion);
for iter=1:itertion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
% Transmitter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
xt_long_with_CP=[];
xt_long_bits=[];

for k=1:N_sym
InfoBits=randsrc(1,2256,[1 0]); % Generate 2256 information bits
xt_long_bits=[xt_long_bits InfoBits];
% Encode and Modulate
CodedBits= encode(henc,InfoBits);
% Rate Matching
ParityLocations=2256+randperm(2256);
ParityLocations=ParityLocations(1:ParityLength);
TotalLocations(:,k)=[1:2256,ParityLocations]; % Should be an output of the encoder
CodedBits=CodedBits(TotalLocations(:,k));
Symbols=ModulateQAM(CodedBits,Modulation);
%%%%pilots%%%%
x_f=zeros(FFTsize,1); x_f(SC_loc)=Symbols;
Xp=randsrc(1,length(PilotLocations),[-1 1]); %Pilot generation
Pilots_PRBS(:,k)=Xp.'; % PRBS patterns
x_f(PilotLocations)=Xp.'; %insertion of the pilots
xt=ifft(x_f);
xt_withCP=[xt(end-1023:end);xt];
xt_long_with_CP=[xt_long_with_CP;xt_withCP];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
% Channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
y=zeros(1,length(xt_long_with_CP));
%%%%Convolution%%%%
for k=max(DelayInSamples)+1:length(xt_long_with_CP)
for j=1:length(TapPower)
a(j)=xt_long_with_CP(k-DelayInSamples(j))*a_samples(j,k);
end
y(k)=sum(a);
end
%%%Noise%%%
sigma_noise=norm(xt_long_with_CP)/sqrt(length(y)*SNR_corrected);
y=y+sigma_noise*(randn(size(y))+1i*randn(size(y)))/sqrt(2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
yBeforeCP_RemoveMatrix=reshape(y.',FFTsize+1024,N_ sym);
yCP_RemoveMatrix=yBeforeCP_RemoveMatrix(1025:end,: );

Y_matrix_PostFFT=fft(yCP_RemoveMatrix);

for kk=1:N_sym,
% H_est = MMSE_CE(Y_matrix_PostFFT(:,kk),Pilots_PRBS(:,kk),P ilotLocations,FFTsize,4,h,SNR);
% h_est = ifft(H_est); h_DFT = h_est(1:channel_length);
% H_DFT = fft(h_DFT,FFTsize);
Y_Pilots=Y_matrix_PostFFT(PilotLocations,kk);

EstTimeDomainChannel=ifft(Y_Pilots.*Pilots_PRBS(:, kk));%Taking only pilots
EstTimeDomainChannel=EstTimeDomainChannel(1:(chann el_length)); % Kill contibutions outside the CP duration
H_est=fft(EstTimeDomainChannel,FFTsize); % est. channel

demodulatedSCs=Y_matrix_PostFFT(:,kk)./H_est;

Y_equalized_withoutPilots(:,kk)=demodulatedSCs(SC_ loc);
ppSNR(:,kk)=abs(H_est(SC_loc)).^2/(sigma_noise^2*FFTsize);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
% Decode Signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%
DecodedBits_long=[];
for kk=1:N_sym,
LLRs=SISO_LLR(Y_equalized_withoutPilots(:,kk),ppSN R(:,kk),Modulation);
LLR_rate_matched=zeros(1,2*2256);
LLR_rate_matched(TotalLocations(:,kk))=LLRs;
DecodedBits=decode(hdec, LLR_rate_matched);
DecodedBits_long=[DecodedBits_long DecodedBits];
end
NumErrors(1,iter)=sum(DecodedBits_long ~=xt_long_bits); % Count errors
end
I still don't understand the utility of iterative loop iter=1:itertion. Based on my understanding to the code, you calculate the BER 'itertion' time for the same configuration. Only random signal, noise and interleaving are modified witch supposed not affecting the BER.

Finely, for me the time taken by the code is related to FFTsize=4096 and Channel effects (LTV+noise) addition witch are applied to all signal xt_long_with_CP.
For IFFT/FFT size we can't change.
For Channel effect, why you don't apply effects to one symbol. It is long enough 2256.

For changing signal to matrix form. encode and decode function will not support matrix input. It is mentioned in matlab help (help fec.ldpcenc/encode)

Regards

Chaker

1 members found this post helpful.

13. ## Re: How to graph BER with perfect results?

For Channel effect, why you don't apply effects to one symbol. It is long enough 2256.
_
Look a variable T, if the symbols/samples(xt_long_with_CP) is to small T is small ->N_low is small -> the channel is very short and it will not affect the Tx signal.

As for the iteration, I did it just because i dont get constant results with the LTV channel. At the end of all the iteration i am getting Numerror vector length 100
with different results. and it's not what i am expecting...

***Thanks for the code, thought like that and implement what you did.

14. ## Re: How to graph BER with perfect results?

Hi,

There is so many constraint in your system.

I get final suggestion.

I propose to use frame notion where 1 frame = 10 * N_sym or more. By this, you can fix y length to 5120*10 for exemple.

Then you add while loop where simulation stop after FrameNb send (FrameNb=100frames) or stop after MaxErrorBits detected (MaxErrorBits=100 error bits).
while ((i_frame<FrameNb)&&(NumErrors<MaxErrorBits))
This loop will eliminate the send off 1000 symbol every simulation.

Good luck WUID

1 members found this post helpful.

15. ## Re: How to graph BER with perfect results?

Muhamed , i liked your suggestion.
Code:
```N_sym=100;
Frame=10*N_sym;
FrameNb=100;
i_frame=1;
ErrorBits=0;
MaxErrorBits=100;
while ((i_frame<FrameNb)&&(ErrorBits<MaxErrorBits))...
for k=1:Frame```

anyone else (or you :) ) on another subject...
On the first channel section where i'm generating the Gaussian random process. I'm having hard times to understand and implement the math. I mean to ask, about my figures and normalization of the filter, taps variance... is it make scene to any1 who is reading this post ?

16. ## Re: How to graph BER with perfect results?

I try to code the solution but i can't simulate.

Code:
```function NumErrors=OurLDPCEncandDec_without_MMSE(itertion,Rate,Modulation,SNRdB)
% % close all; clc ;
%  clear all ;
%  Rate=2;
% Modulation=4;
% SNRdB=3;
% itertion=5;
% Rate: 2 for 1/2, 4 for 3/4, 6 for 5/6, 8 for 7/8

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters Section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Symbols_length=2256;
CP_lenght=1024;
N_sym=100;
Frame=10*N_sym;
FrameNb=100;
MaxErrorBits=1000;

if Rate==2
ParityLength=2256;
elseif Rate==4
ParityLength=2256/3;
elseif Rate==6
ParityLength=464;
elseif Rate==8
ParityLength=320;
elseif Rate==10
ParityLength=240;
end;

H1=H(:,1:2256);
H2=H(:,2257:4512);
H=[H2 H1];
henc = fec.ldpcenc(H);
hdec = fec.ldpcdec(H);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IFFT/FFT Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FFTsize=2^(ceil(log2(Symbols_length)));
ActiveSCs=Symbols_length;
PilotLocations=1:4:FFTsize; %pilot loc /for mod=4,16 spacing = 3. /for mod=64 spacing =4
SC_loc=setdiff(1:FFTsize,PilotLocations); SC_loc=SC_loc(1:ActiveSCs); %info location
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs=10e6; %Sampling rate
N_samples=N_sym*(FFTsize+CP_lenght);
T=N_samples/Fs; %Sample time

% Taking channel parameters from ITU table
TapPowerdB=[0 -9.7 -19.2 -22.8];   % in dB
Delay=[1 110 190 410]*1e-9; % in Sec
% ChannelLength=length(TapPower);

a_samples=[];
%%%%Tap power in dB to linear%%%%
TapPower=10.^(TapPowerdB/10);
% TapPower=TapPower/norm(TapPower); % Normalize to unit channel power

F_low=f_m*10; % usaully 10*f_m the low sampling
% freq of the channel coefficients
N_low=ceil(T*F_low);

%%%%Gausssian random process pass throw filter(flat/U-shaped)%%%%

h=firpm(64,[0 2*f_m/F_low 2*f_m/F_low*1.2 1],[1 1 0 0 ]);
h=h/norm(h);
%        channel_length=length(h);
for k=1:length(TapPower)

a_low_temp=conv((randn(1,N_low+64)+1i*randn(1,N_low+64))/sqrt(2),h);%was ...,h)*TapPower(k);
a_low_temp=a_low_temp(65:N_low+64); % Generate first tap sampled at low freq
a_low(k,:)=a_low_temp;
channel_length=6;
%%%%Upsample To Fs%%%%

a_samples(k,:)=resample(a_low(k,:),Fs,F_low);
Sigvar=var(a_samples(k,:));
a_samples(k,:)=a_samples(k,:)*sqrt((TapPower(k)/Sigvar));
end

%%%%Initilaize Delay Vector%%%%
DelayInSamples=ceil(Delay*Fs);
DelayInSamples=DelayInSamples-DelayInSamples(1);

%%% Noise parameters
SNR=10^(SNRdB/10);
SNR_corrected=SNR*(1366+2256)/FFTsize; %need to generic fit for qpsk 2256 infobits + 1366 pilots

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Began iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NumErrors=zeros(1,itertion);
for iter=1:itertion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Transmitter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i_frame=1;
ErrorBits_all_frame=0;

while ((i_frame<FrameNb)&&(ErrorBits_all_frame<MaxErrorBits))

xt_long_with_CP=[];
xt_long_bits=[];

for k=1:N_sym
InfoBits=randsrc(1,2256,[1 0]); % Generate 2256 information bits
xt_long_bits=[xt_long_bits InfoBits];
% Encode and Modulate
CodedBits= encode(henc,InfoBits);
% Rate Matching
ParityLocations=2256+randperm(2256);
ParityLocations=ParityLocations(1:ParityLength);
TotalLocations(:,k)=[1:2256,ParityLocations]; % Should be an output of the encoder
CodedBits=CodedBits(TotalLocations(:,k));
Symbols=ModulateQAM(CodedBits,Modulation);
%%%%pilots%%%%
x_f=zeros(FFTsize,1); x_f(SC_loc)=Symbols;
Xp=randsrc(1,length(PilotLocations),[-1 1]);    %Pilot generation
Pilots_PRBS(:,k)=Xp.'; % PRBS patterns
x_f(PilotLocations)=Xp.'; %insertion of the pilots
xt=ifft(x_f);
xt_withCP=[xt(end-1023:end);xt];
xt_long_with_CP=[xt_long_with_CP;xt_withCP];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y=zeros(1,length(xt_long_with_CP));
%%%%Convolution%%%%
for k=max(DelayInSamples)+1:length(xt_long_with_CP)
for j=1:length(TapPower)
a(j)=xt_long_with_CP(k-DelayInSamples(j))*a_samples(j,k);
end
y(k)=sum(a);
end
%%%Noise%%%
sigma_noise=norm(xt_long_with_CP)/sqrt(length(y)*SNR_corrected);
y=y+sigma_noise*(randn(size(y))+1i*randn(size(y)))/sqrt(2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yBeforeCP_RemoveMatrix=reshape(y.',FFTsize+1024,N_sym);
yCP_RemoveMatrix=yBeforeCP_RemoveMatrix(1025:end,:);

Y_matrix_PostFFT=fft(yCP_RemoveMatrix);

for kk=1:N_sym,
%     H_est = MMSE_CE(Y_matrix_PostFFT(:,kk),Pilots_PRBS(:,kk),PilotLocations,FFTsize,4,h,SNR);
%     h_est = ifft(H_est); h_DFT = h_est(1:channel_length);
%     H_DFT = fft(h_DFT,FFTsize);
Y_Pilots=Y_matrix_PostFFT(PilotLocations,kk);

EstTimeDomainChannel=ifft(Y_Pilots.*Pilots_PRBS(:,kk));%Taking only pilots
EstTimeDomainChannel=EstTimeDomainChannel(1:(channel_length)); % Kill contibutions outside the CP duration
H_est=fft(EstTimeDomainChannel,FFTsize); % est. channel

demodulatedSCs=Y_matrix_PostFFT(:,kk)./H_est;

Y_equalized_withoutPilots(:,kk)=demodulatedSCs(SC_loc);
ppSNR(:,kk)=abs(H_est(SC_loc)).^2/(sigma_noise^2*FFTsize);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Decode Signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DecodedBits_long=[];
for kk=1:N_sym,
LLRs=SISO_LLR(Y_equalized_withoutPilots(:,kk),ppSNR(:,kk),Modulation);
LLR_rate_matched=zeros(1,2*2256);
LLR_rate_matched(TotalLocations(:,kk))=LLRs;
DecodedBits=decode(hdec, LLR_rate_matched);
DecodedBits_long=[DecodedBits_long DecodedBits];
end
ErrorBits_one_frame=sum(DecodedBits_long ~=xt_long_bits); % Count errors
ErrorBits_all_frame=ErrorBits_all_frame+ErrorBits_one_frame;
i_frame=i_frame+1;
end % End while loop (all frame send or 100 errors bits detected)
NumErrors(1,iter)=ErrorBits_all_frame;
all_transmitted_bits=(i_frame-1)*N_samples;
ber_values(1,iter)=ErrorBits_all_frame/all_transmitted_bits;
end % End iter loop```
Note that the number of transmitted bits is variable.
Code:
```NumErrors(1,iter)=ErrorBits_all_frame;
all_transmitted_bits=(i_frame-1)*N_samples;
ber_values(1,iter)=ErrorBits_all_frame/all_transmitted_bits;```
Try to simulate and give me feedback :)

For your question, concerning the AWGN noise and OFDM yes but other part no :/

best regards

1 members found this post helpful.

17. ## Re: How to graph BER with perfect results?

thx Mohamed!!!
for now I've dropped down the iteration part for now. as for
Code:
`(i_frame-1)*N_samples`
it's wrong. for all type of modulation I am transmitting 2256 InfoBits. so
Frame=10*N_sym
Totalbits=2256*Frame*FrameNb

can u tell me about the "other part"? what do you think ? or is it completely wrong ?

***attaching my simulation code with your suggestion and my modifications including SISO_LLR function :)

18. ## Re: How to graph BER with perfect results?

Hi WUID,

I have done some simulation.

I get ber=4 10^-2 for SNRdB=0 and ber=5 10^-5 for SNRdB=1.

According to ieee article (see files attached), this is not true.

I think there is a problem with awgn noise generation.

Began with verifying SNR value.

2 members found this post helpful.

19. ## Re: How to graph BER with perfect results?

ya i know it looks that the encoder doing some magic and in addition you can see that the SNR is being corrected in SNR_corrected to the number of fft bins i am using e.g
Code:
`SNR_corrected=SNRdB*(1366 pilots + 2256 data bins)/FFTsize`
-> so the SNRdB value decrease

•

20. ## Re: How to graph BER with perfect results?

There is no magic in digital communication lol
I think the ber you get it is not for the snr declared.
To have corresponding values, you must correct the SNR expression.

I used this code to simulate AWGN
Code:
```EbNo = 0:2:15;
EsNo= EbNo+10*log10(Used_subcarrier/len_fft)+ 10*log10(len_fft/Symbol_interval) + 10*log10(k);
snr=EsNo - 10*log10(len_fft/Symbol_interval);
chan_awgn = awgn(ser_data,snr(ii),'measured');```