How to graph BER with perfect results?

WUID · Mar 17, 2012

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 '[]'.

WUID · Mar 22, 2012

anyone can help me with that ?

capcas · Mar 27, 2012

WUID said:
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.

rhearty1 · Mar 29, 2012

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

ChakerQ · Mar 30, 2012

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

WUID · Mar 31, 2012

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 :
instead of
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 ? :???:

ChakerQ · Mar 31, 2012

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 ?

WUID · Mar 31, 2012

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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load 802_16e_4512_d5_H H -mat;

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
    f_m=10;  %Doppler spread
    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 
    
   


% Add Noise 
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 ----------

Mohamed Chaker said:
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)

ChakerQ · Apr 1, 2012

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

rhearty1 · Apr 1, 2012

rhearty1 said:
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:
corr = sum(received_sig.*template_sig);

This would do the same as:
corr = 0;
for ii=1:LENGTH
corr = corr + received_sig(ii)*template_sig(ii);
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

WUID · Apr 1, 2012

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.

ChakerQ · Apr 1, 2012

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

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;

load 802_16e_4512_d5_H H -mat;
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
f_m=10; %Doppler spread
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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(1arityLength);
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*SNR_corrected);
y=y+sigma_noise*(randn(size)+1i*randn(size))/sqrt(2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Receiver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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(1channel_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
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

WUID · Apr 1, 2012

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.

ChakerQ · Apr 1, 2012

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

WUID · Apr 1, 2012

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

Hope that's what your meant...

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 ?

ChakerQ · Apr 1, 2012

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;

load 802_16e_4512_d5_H H -mat;
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
f_m=10;  %Doppler spread
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); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Receiver
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
that is why i add

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

WUID · Apr 1, 2012

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

ChakerQ · Apr 1, 2012

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.

WUID · Apr 1, 2012

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

ChakerQ · Apr 1, 2012

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

I follow this tutorial
http://www.dsplog.com/2008/08/26/ofdm-rayleigh-channel-ber-bpsk/

Welcome to EDAboard.com

How to graph BER with perfect results?

Member level 2

Member level 2

Member level 3

Junior Member level 1

Member level 1

Member level 2

Member level 1

Member level 2

Member level 1

Junior Member level 1

Member level 2

Member level 1

Member level 2

Member level 1

Member level 2

Member level 1

Member level 2

Attachments

Member level 1

Attachments

Member level 2

Member level 1

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor