Matlab UMTS/LTE Turbo Code
This Matlab code plots the BER, EXIT chart and iterative decoding trajectories for the UMTS and LTE turbo codes, when using BPSK modulation for transmission over an AWGN channel. Functions are provided to generate the UMTS and LTE interleavers, as well as to perform component encoding and decoding. Furthermore, functions are provided for generating a priori LLRs having a particular mutual information, as well as for measuring the mutual information of some LLRs.
main_ber.m draws a BER plot for the UMTS turbo code.
main_exit.m draws an EXIT chart for the UMTS turbo code.
main_traj.m plots the iterative decoding trajectories of the UMTS turbo code.
get_UMTS_interleaver.m provides a function for generating the UMTS interleaver.
component_encoder.m provides an encoder function for the UMTS component codes.
component_decoder.m provides a BCJR decoder function for the UMTS component codes.
jac.m provides a function for performing the exact, lookup-table-aided and approximate Jacobian logarithms.
generate_llrs.m provides a function for generating Gaussian distributed a priori LLRs.
measure_mutual_information_histogram.m measures the mutual information of some LLRs using the histogram method.
measure_mutual_information_averaging.m measures the mutual information of some LLRs using the averaging method.
You can download the Matlab code here. You can also download the LTE interleaver. The operation of the UMTS turbo decoder is described in Section 2.2 of Liang Li’s nine month report and this document.
Copyright © 2010 Robert G. Maunder. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
August 12th, 2010 at 8:29 am
Dear Rob,
In your main_exit.m file in Matlab UMTS Turbo code, the exit chart is based on the figure 2.13. Is this means that the exit chart plotted is only for one decoder? I understand that by swapped the axis, we can display the second curve since the output is the input of the second decoder. I just wonder if we need to repeat the process of figure 2.13 in figure 2.11. If not, the function of interleaver or deinterleaver is not considered, am I right?
August 12th, 2010 at 8:34 am
Hello again Annena,
In an EXIT chart, each EXIT function depends only on one of the decoders - it is completely independent of the other decoder and of the interleaver. Since the upper and lower codes in a UMTS turbo code are identical, their EXIT functions are also identical. For this reason, there is no need to run two simulations. Instead, you can just take the EXIT function for the upper code and invert its axes to get the EXIT function of the lower code. Please note that while the EXIT functions are completely independent of the interleaver, the iterative decoding trajectories do depend on the interleaver design. The trajectory will only match the EXIT functions if the interleaver is able to randomise the order of the LLRs sufficiently well.
Hope this helps, Rob.
August 12th, 2010 at 8:39 am
Dear Rob,
Thank you so much, now I clearly understand. I think, I need to check back the main_traj.m file first.
August 13th, 2010 at 7:35 am
Hi Rob,
May I know, do you consider channel reliability (Lc) in your BCJR decoding?
Annena
August 13th, 2010 at 8:47 am
Hi again,
In case I want to use your component_decoder.m function for different generator, what is the modification should I consider? Is it OK by only changes the trellis matrix?
Annena
August 13th, 2010 at 9:22 am
Hi Annena,
In the main_ber.m, main_exit.m and main_traj.m files, I use lines like the following to perform soft BPSK demodulation…
a_c = (abs(a_rx+1).^2-abs(a_rx-1).^2)/N0;
If you expand the squares, you can rearrange this to a form including the channel reliability…
a_c = Lc*real(a_rx);
… where…
Lc = -4/N0;
If you want to change the generator polynomials, all you need to do is change the transitions matrix in component_decoder.m and the equations in component_encoder.m.
Hope this helps, Rob.
August 18th, 2010 at 2:51 am
Hello Rob,
I want to know about the main_exit.m file on the encoder part. The data is transmitted using rate 1/2. in the file, you are transmitting systematic bit and parity bit from the first encoder only, right? This means that you did not apply the puncturing method in the encoder.
Please correct me if I’m misunderstood your program.
TQ.
August 18th, 2010 at 7:53 am
Hi Annena,
You are correct. Only the information that is used by the upper decoder is transmitted in this simulation. This allows the upper decoder to be characterised in isolation from the lower decoder. The simulation required to characterise the lower decoder is identical to this one. That is why the results are plotted twice, once for the upper decoder and inverted for the lower decoder. By the way, the UMTS turbo code does not use puncturing.
Take care, Rob.
August 19th, 2010 at 3:03 am
Hi,
Tq so much for your help.
I want to know,
1. if I run the main_exit.m file using my encoder and decoder, and the result is decreasing for an increasing IA_index, what is it means? I believe it must be something wrong with my program in the decoder part, but is it explain anything?
2. if we referring to standard parallel concatenated code in paper stephen ten brick fig1, 1st BCJR decoder require input of Z1(systematic and p1) and A1.
a) If I want to use your main_exit.m file, is it means that A1=a_a?
b) the length of Z1 is equal to the length of (a_c+c_c + e_c(terminated bits)). Since a_a is generated from a, which is without terminated bits, the length of a_a is not equal to Z1. If A1=a_a, can I just add bit zeros to a_a so that it will have the same length as Z1?
I’ll really appreciate your help.
August 19th, 2010 at 9:32 am
Hi Annena,
1. It sounds like there is a problem with either your encoder or your decoder. You can confirm this by comparing the EXIT functions that you get using the averaging method of measuring mutual information with the EXIT functions you get using the histogram method. If these EXIT functions are significantly different then it implies that there is a bug in your encoder, modulator, demodulator or your decoder.
2a. That’s right. Here, I’m assuming that you mean this paper…
http://scholar.google.co.uk/scholar?cluster=15764897033271139190&hl=en&as_sdt=2000
2b. Note that the length of A1 is not equal to the length of Z1 in that paper - it is half the length of Z1. The operation of main_exit, main_ber and main_traj is illustrated in Figures 2.11 and 2.13 of…
http://users.ecs.soton.ac.uk/rm/wp-content/liang_li_nine_month_report.pdf
Here, if the interleaver has a length of N bits, then a_c has a length of N, c_c and d_c have lengths of N+3, while e_c and f_c have lengths of 3. Instead of zeros as you suggest, e_c is concatenated onto a_a+a_c in order to give y_a, which has the same length as c_c of N+3.
Hope this helps, Rob.
August 23rd, 2010 at 2:02 am
Hi Rob,
Thank you so much for your help, I feel so happy because finally I manage to plot my first EXIT chart using my encoder and decoder.
Now, I’m learning to plot the trajectory part.
1) May I know, why is the file main_traj.m plot trajectory at different frame? Should we just average the value so that we can get just a smooth trajectory plot?
2) From the trajectory figure, can we plot the EXIT curve by taking the peak of each trajectory point or do we need to map it into figure from main_exit.m?
Really appreciate your help.
Regards,
Annena
August 23rd, 2010 at 9:00 am
Hi Annena,
1) You can average the trajectories, but you sometimes find that the result doesn’t give a very good match with the EXIT functions. Also, if you take the average of the trajectories, you don’t get a feel for how much variation there is from frame to frame. It is this variation that explains why short frame lengths give a poorer performance than longer frames.
2) You can do this, but the result will typically not be very smooth and will also depend on your interleaver design. Furthermore, the whole point of EXIT chart analysis is that is can be done without requiring lengthy iterative decoding simulations.
Hope this helps, Rob.
August 27th, 2010 at 12:35 am
Hi Rob,
I noticed that your BPSK encoder modulate the bit 1 to -1 and the bit 0 to 1. In your component_decoder.m file, extrinsic_uncoded_llrs(bit_index) = prob0-prob1, does it mean prob0 is for bit 1 and prob1 is for bit -1? in case I modulate it otherwise, 1–> 1 and 0 –>-1, can I still used your decoder file?
Appreciate your advise.
August 27th, 2010 at 8:48 am
Hi Annena,
I find that this causes a lot of confusion for people. It all comes down to whether LLRs are defined as LLR=ln(P0/P1) or LLR=ln(P1/P0). J Hagenauer uses the former definition and C Berrou uses the latter one. My code uses the definition LLR=ln(P0/P1). If another part of your system uses the other definition (such as your BPSK demodulator), all you need to do is multiply your LLRs by -1 when you pass them from part of the system to the other.
Hope this helps, Rob.
November 11th, 2010 at 8:29 am
Hi Rob,
I develop my own encoder file where only my first encoder is terminated and the second encoder is not. To develop the decoder, how can I modify your component_decoder.m file?
Please advise.
TQ
November 11th, 2010 at 12:13 pm
Hi Mila,
All you need to do for the unterminated decoder is change the line…
betas(1,length(apriori_uncoded_llrs))=0;
…to…
betas(:,length(apriori_uncoded_llrs))=0;
Take care, Rob.
November 22nd, 2010 at 2:46 am
Hi Rob,
Thank you for your respond.
May I know, is it possible to plot the EXIT chart using your function when the first decoder is terminated and the second decoder is not?
As I know, for parallel turbo code, we only used one decoder only to plot the EXIT chart.
Please advise
TQ
November 28th, 2010 at 5:35 am
Hi Mila,
You should fnd that termination makes only a negligible difference to the EXIT function, particularly if you have a long interleaved. My advice would be to simply mirror one or other of the EXIT functions - either the terminated one, or the unterminated one.
Take care, Rob.
February 14th, 2011 at 8:40 pm
Hi Bob,
1. Thank you so much for your code. It is really hard to find something so helpful about the BCJR algorithm in MATLAB on the Internet.
2. I have a question. Now I want to implement a conv code decoder with BCJR based on your code, but in my scenario the rate of the code is different, like 2/3.
(1) Are the equations in Liang Li’s report still applicable?
(2) I know I need to change the “transitions” matrix; for the “UncodedBit” and “EncodedBit” column, do I need to expand them into several binary column, or simply use decimal representations of the bits, such as 7 for 111?
(3) Other things I need to care about?
Thank you very much for your help.
February 16th, 2011 at 8:59 pm
Hello Lai,
The equations in Liang Li’s report are still applicable, but they need to be generalized. You can see generalized versions at…
http://users.ecs.soton.ac.uk/rm/wp-content/LogBCJR.pdf
You need to insert an additional column into the transitions matrix for each additional input and output bit. For example, a k=2, n=3 CC would have a transitions matrix with 7 columns: two for the states, two for the input bits and three for the output bits. You also need five sets of gamma calculations, one for each input and output bit. Finally, you need two sets of gamma and extrinsic LLR calculations, one for each input bit.
Hope this helps, Rob.
February 23rd, 2011 at 10:04 am
Somebody has asked me how to concatenate the UMTS turbo decode with an additional inner decoder, such as a SISO MIMO detector. I have drawn a schematic for this and uploaded it to…
http://users.ecs.soton.ac.uk/rm/wp-content/concat.png
Take care, Rob.
June 15th, 2011 at 7:31 pm
Hi Rob,
Fine algorithm, I’ve been using it to run some tests within a larger system. However, ECC are not my domain and I need to obtain the LLRs of the coded bits (lc+ld+le+lf) and not just the LLRs from the data bits (la) at the end of the turbo decoding step. Howmay I achieve it? Thanks in advance.
June 16th, 2011 at 9:26 am
Hi José,
This is not too difficult to achieve. You just need to modify component_decoder.m so that it has an additional output called extrinsic_encoded_llrs. You can generate this using these lines of code
% Calculate encoded extrinsic transition log-confidences. This is similar to
% Equation 2.18 in Liang Li’s nine month report or Equation 4 in the BCJR paper.
deltas2=zeros(size(transitions,1),length(apriori_encoded_llrs));
for bit_index = 1:length(apriori_encoded_llrs)
for transition_index = 1:size(transitions,1)
deltas2(transition_index, bit_index) = alphas(transitions(transition_index,1),bit_index) + uncoded_gammas(transition_index, bit_index) + betas(transitions(transition_index,2),bit_index);
end
end
% Calculate the encoded extrinsic LLRs. This is similar to Equation 2.19 in
% Liang Li’s nine month report.
extrinsic_encoded_llrs = zeros(1,length(apriori_encoded_llrs));
for bit_index = 1:length(apriori_encoded_llrs)
prob0=-inf;
prob1=-inf;
for transition_index = 1:size(transitions,1)
if transitions(transition_index,4)==0
prob0 = jac(prob0, deltas2(transition_index,bit_index));
else
prob1 = jac(prob1, deltas2(transition_index,bit_index));
end
end
extrinsic_encoded_llrs(bit_index) = prob0-prob1;
end
June 16th, 2011 at 8:58 pm
Hi Rob,
Thanks for your help, it works just fine. Much apprettiated.
September 20th, 2011 at 3:48 pm
Dear Rob
Can you plase tell somthing about turbo codes, that as why is it necessary that the information feeded into encoder should be unocrrelated as much as possible.
What if they are correlated, apart from turbo decoding how does it affect. Can you please explain it
thanks
swap
September 20th, 2011 at 4:33 pm
Hi Swap,
If the occurrence of two events A and B is correlated then their joint probability is given by P(A,B) = P(A|B)*P(B). If they are uncorrelated then the joint probability is given by P(A,B) = P(A)*P(B). The second situation is much easier to work with because we only need one P(A) value for each outcome of A. By contrast, the first situation requires a P(A|B) value for each possible combination of outcomes of A and B. To save on complexity, the BCJR assumes that P(A,B) = P(A)*P(B), i.e. that the LLRs are uncorrelated. If they are correlated, then this assumption is false and the performance of the BCJR suffers.
Take care, Rob.
October 30th, 2011 at 2:32 am
Hi Rob,
First of all, thanks for the code, it’s really helpful. For component turbo code, you have separated the information, parity as well as the termination bits when input to channel. I wonder whether the result will be the same if they are combined into one before input to the channel?
The reason for me asking the question because I’m not sure whether the Randn function in Matlab is independent statistically when you call it for the second or subsequent time. Let say it’s independent, therefore for your case, the AWGN noise is independent for each the systematic, parity and termination sequence (a_rx,c_rx_d_rx,e_rx and f_rx)?
Looking forward for your explanation.
Many Thanks,
Mory
October 31st, 2011 at 9:52 am
Hi Mory,
I kept the different bit sequences separate because it is very easy to make a mistake when programming their concatenation and de-concatenation. If they were concatenated, all the results would be exactly the same.
Take care, Rob.
November 3rd, 2011 at 10:08 pm
Hello,
Thanks for your great work & handy reference. I think I found a tiny typo/bug in the code: in main_exit.m, the value 0.3037 should be 0.3073 ( from http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1499068 ).
Regards,
-Colin
November 4th, 2011 at 10:03 am
Hi Colin,
Thank you for this. However, I’m not sure where in my code or the paper you mean - I can’t find any values of 0.3037 or 0.3073 in my code or that paper.
Take care, Rob.
November 8th, 2011 at 4:03 am
oops… I meant generate_llrs.m
November 8th, 2011 at 11:34 am
Hi Colin,
Thanks for noticing this typo - I have updated all the occurrences of this mistake in the code on my website (I think). This shouldn’t make any difference to the results though, because of the while loop in generate_llrs.m.
Take care, Rob.
March 6th, 2012 at 3:57 am
Dear Rob,
I have some questions in component_decoder.m,
why there are two parts: uncoded_gammas and encoded_gammas, while in bcjr_decoder is just gammas . And why uncoded_gammas is computed in alphas and betas,but not deltas. By the way , in the computaion of alphas or betas , why the bit_index -1 is used ?(encoded_gammas(transition_index, bit_index-1)) Many thanks.
March 6th, 2012 at 6:37 am
alphas(transitions(transition_index,2),bit_index) = jac(alphas(transitions(transition_index,2),bit_index),alphas(transitions(transition_index,1),bit_index-1) + uncoded_gammas(transition_index, bit_index-1) + encoded_gammas(transition_index, bit_index-1));
Based on equ. 5 in BCJR , uncoded_gammas(transition_index, bit_index) + encoded_gammas(transition_index, bit_index). Is it right?
But I don\’n understand why just the encoded_gammas is used to compute the deltas in component_decoder and both encoded_gammas and uncoded_gammas are used in bcjr. Wait for your answer . Thanks a million.
March 6th, 2012 at 8:09 am
Dear Rob,
In main_ber :a_p = a_a + a_c + a_e;
In main_traj :a_p = a_a + a_e;
So why not the same?
March 6th, 2012 at 11:21 am
Hi Lingjun,
This is the difference between a BCJR decoder that generates extrinsic LLRs and one that generates aposteriori LLRs. A decoder that generates aposteriori LLRs has the advantage of not needing to make a distinction between uncoded_gammas and encoded_gammas, therefore giving a lower complexity. However, there is a small disadvantage, which relates to the fact that the extrinsic LLRs must be calculated by subtracting the aproiri LLRs from the extrinsic LLRs - if the apriori and extrinsic LLRs are both infinite in magnitude, then subtracting one from the other gives an undetermined result. This is a very minor disadvantage and so in practice, decoders that generate aposteriori LLRs are preferred, owing to their lower complexity.
By the way, the line a_p = a_a + a_e; in main_traj.m should say a_p = a_a + a_c + a_e; I have fixed this bug and uploaded a new version of the code.
Take care, Rob.
April 12th, 2012 at 12:12 am
Dear Rob,
I want to implement puncturing using LDPC code in two parallel (AWGN) channels. I have the concept of puncturing but I am feeling difficulties in implementing.
could you suggest me any help in this matter please please
April 12th, 2012 at 9:12 am
Hi Ideal,
To puncture a bit sequence, you just remove some of the bits in the transmitter. To perform the corresponding depuncturing operation in the receiver, you just replace the punctured bits with zero-valued LLRs.
Take care, Rob.
April 12th, 2012 at 12:29 pm
Hi Rob,
That is possible for single channel. Actually I want to implement the puncturing in two channel (parallel optical and RF channel).
I think there would be a way to generate a puncturing pattern using single encoder/decoder.
I can assume the channel state information at the transmitter and can puncture the codeword accordingly but what can be the ensemble and how to send the punctured bits in different channel is a challeng for me. Thats why I am looking such a technique which can give me inside help about the puncturing. Thanks
April 12th, 2012 at 12:31 pm
Hi Ideal,
I think that you are entering new ground with what you are talking about and so I think that you will have to investigate how best to do the puncturing yourself.
Take care, Rob.
April 13th, 2012 at 2:12 am
Hi Rob,
I got few ideas but one thing,
If i puncture some bits at the trasmitter side then for the depuncturing on the receiver side, how can I add zeros at the exact locations, i.e.
let x= [1 0 1 1 1 1 0 1] is a codeword and I puncture either random or few bits from x.
let few bits punctured xp = [p 0 1 1 1 p 0 p], it means i puncture 1, 6 and 8 bit, how the receiver knows about it, I mean in matlab how can I write a command for this.
secondly If I have a (half rate, r=0.5) codeword of length 1024 and I punctured 171 bits randomly out of 1024 to have a rate of rp=0.6. now how can i introduced zeros at the particular values of punctured codeword.
The simple question is in random puncturing, how can i write a command in matlab which can puncture the bits at the transmitter and add zeros at those locations (punctured) at the reciever side. please send me this type of matlab script please. thanks
April 13th, 2012 at 9:39 am
Hi Ideal,
The receiver needs to have knowledge of the puncturing pattern. This code does everything you are talking about - it also interleaves the bits, which is useful in a Rayleigh fading channel.
temp = randperm(1024);
pattern = temp(1:853);
% In the transmitter
punctured_bits = bits(pattern);
% In the receiver
llrs = zeros(size(bits));
llrs(pattern) = punctured_llrs;
Take care, Rob.
April 15th, 2012 at 6:30 am
Dear Rob,
I really very much appreciate your help.
one things more please,
1. bits are the coded bits after the encoder?
2. In line, llrs(pattern) = punctured_llrs; you didn\’t defined punctured the punctured_llrs, I think, punctured_llrs = (2/sigma^2)*(punctured_bits + noise); for awgn channel?
Am I right in the above two questions?
Again I am really very appreciate your help. Thanks
April 16th, 2012 at 9:05 am
Hi Ideal,
Both of your assumptions are correct.
Take care, Rob.
April 17th, 2012 at 4:05 am
Dear Rob,
Thanks for the answer, One more question please
I just want to understand the concept of capacity approaching low density parity check codes, (some researcher are calculating the bit error rate curve for certain LDPC code and claim that its capacity approaching code)
can you have any suggestion/advise for that things?
Thanks for that
April 17th, 2012 at 10:16 am
Hi Ideal,
You should look into Irregular LDPC codes - these can operate extremely closely to the Shannon capacity.
Take care, Rob.
April 18th, 2012 at 3:23 am
Dear Rob,
Currently I am interested in low rate regular structure LDPC code and then introduce ouncturing to in increase the code rate. using simple bpsk channel.
But I am confuse with the threshold of code and capacity. Like I don’t understand, how to calculate the threshold of the code and how to prove that it is capacity approaching code for simple BPSK AWGN channels and even the punctured code with higher rate is performed well. These are the questions I have in my mind please
Thanks for the response
April 18th, 2012 at 8:30 am
Hi Ideal,
You can draw the EXIT chart of your LDPC codes and find the lowest SNR at which the EXIT chart tunnel remains open. Then, you can compare this SNR with that at which the capacity of the channel becomes equal to the throughput of your LDPC code. If the SNRs are similar, then your code is a near-capacity code.
Take care, Rob.
April 18th, 2012 at 11:55 pm
Dear Rob,
Thanks now I understand the threshold. Thanks
What is the relation between the throuhgput of LDPC code and the channel capacity (sorry for silly questions).
According to me, throughput = symbol rate x capacity, but I am not sure, how to calculate the throughput of LDPC code and compare with channel capacity?
I can calculate the EXIT chart, bpsk channel capacity but not sure how to compare the throughput of LDPC code with channel capacity at particular snr?
Thanks for help
April 19th, 2012 at 8:57 am
Hi Ideal,
The throughput is measured in bits per symbol (or equivalently bits per second per hertz) - it is given by R*log2(M), where R is the LDPC coding rate and M is the number of constellation points in your modulation scheme, e.g. M=16 in 16QAM. The capacity is plotted at http://users.ecs.soton.ac.uk/rm/resources/matlabcapacity/
Take care, Rob.
April 20th, 2012 at 12:22 am
Dear Rob,
I understand the calculation of capacity. Thanks
I understand the threshold will be that value of snr where the check node curve and the variable node curve just touches.
Also As we know that I_A = J(sigma), where sigma = sqrt(8*R*Es/No)?
sorry for not getting the point
1. EXIT chart (snr) = threshold (where two curves touches)
2. channel capacity (snr(Threshold)),
you said (Then, I can compare this SNR(threshold) with that at which the capacity of the channel becomes equal to the throughput (R*log2(M) of your LDPC code. If the SNRs are similar, then your code is a near-capacity code.)
It means for half rate code and bpsk modulation, the LDPC throughput will be 0.5. Then for the near capacity code, the channel capacity of bpsk channel should be near 0.5 at snr(threshold). I think I am right now
April 20th, 2012 at 5:04 pm
Hi Ideal, You’ve got it. Rob.
April 22nd, 2012 at 12:49 am
Dear Rob,
Can you please tell me any code of generating the ldpc parity check matrix (H) and the corresponding generator matrix (G).
I know about the PEG algorithm for generating H matrix but I want some other script if you know please. Thanks
April 23rd, 2012 at 9:00 am
DEAR ROB,
IS THERE ANY EASY WAY TO CALCULATE THE DECODER TRAJECTORY (STAIRS) IN MATLAB? THANKS
April 23rd, 2012 at 11:28 am
Hi Ideal,
I have uploaded some code for generating G and H matrices to…
http://users.ecs.soton.ac.uk/rm/wp-content/generate_H.m
The trouble with this code is that it tries too hard to find matrices that have exactly the specified degrees, and so it takes a long time to run. You should probably modify it to make it quicker to run.
You can download main_traj.m from the top of this page - it shows you how to plot the trajectory.
Take care, Rob.
April 23rd, 2012 at 12:01 pm
Dear Rob,
Its really great help. I will look in the code. Thanks dear Sir,
April 24th, 2012 at 8:02 pm
Hi Rob,
I’m trying to use your BCJR to decode my 1/3 Turbo coded stream.
I don’t understand the way you’re calculating “gammas” (Equation 2.14 in the 9 month report). Specificly, how you derive that (1-y(t)) .. formula and why don’t you directly use the conditional pdf instead? Is this specific to UMTS code?
Could you please explain?
Thanks,
Elnaz
April 25th, 2012 at 12:11 am
Hi Rob,
I have another question:
In calculating “betas”, according to the report you cited, shouldn’t we calculate each “betas” using the “gammas” value of the same stage i.e. to follow backward recursive move? In your code you are using the “gammas” value for the next stage. Am I missing something here?
I see the same thing when you caculate “deltas” i.e. “alphas” and “gammas” and “betas” are all indexed with “bit_index” but I think “betas” should be indexed by “bit_index”+1 (one stage ahead).
Thanks,
Elnaz
April 25th, 2012 at 1:08 am
Dear Rob,
I try to understand the code of main_traj.m but couldn’t. I am using the LDPC code and its trajectory will be the information exchange between the variable and check nodes. So i think there should be simple program in order to calculate the trajectory without using the sum-product algorithm.
At the moment I am using the Sum-Product algorithm and calculating the trajectory based on the information exchange when the SP algorithm runs. I was thinking to calculate the trajectory without running the decoder. I am not sure if it is possible? Thanks
April 25th, 2012 at 10:58 am
Hi Ideal,
If you have the EXIT chart then it is possible to predict the trajectory, without running the iterative decoder. You just need to manually draw a staircase that bounces between the two EXIT functions. The trouble is that this predicted trajectory will only be valid in the case of a very long block length. This approach will not let you consider the variation in the trajectory that results when using shorter block lengths.
Take care, Rob.
April 25th, 2012 at 11:05 am
Hi Elnaz,
Equation 2.14 is saying that if a transition is labelled with a zero, then the corresponding gamma should take on the value of the corresponding a priori LLR. If the transition is labelled with a one, then the corresponding gamma should take on the value of zero.
This approach has a slightly lower complexity than the one originally proposed for the Log-MAP. This says that if a transition is labelled with a zero, then the corresponding gamma should take on the value of the corresponding a priori LLR divided by 2. If the transition is labelled with a one, then the corresponding gamma should take on the value of the corresponding a priori LLR divided by -2.
Betas should be calculated using the gammas of the transitions to the right of the corresponding states. I think that the notation in the original explanation of the Log-MAP may be confusing you. It is because this notation is so confusing that we decided to explain things in this different way.
Take care, Rob.
April 25th, 2012 at 1:35 pm
Hi Rob,
Thanks for your explanation.
As for my first question, I don’t see you divide LLRs by 2? The way I explained it to my self is that you’re normalizong LLRs with respect to the probability of a bit value being 1 and zero (if you know what I mean).
Actually, this brings another question because in this normalization you have LLR= P(bit=1)/P(bit=0), but at the end of the program you calculate extrinsic LLR with P(bit=0)/P(bit=1) ratio?
Also, why are you omitting “uncoded-gammas” in calculating “deltas” ?
Thank you for answering my questions,
Elnaz
April 25th, 2012 at 1:58 pm
Hi Elnaz,
There’s no need to dived by two if you do it the way I do here - this is because the absolute value of the gammas doesn’t matter - it’s the difference between the gammas associated with one and zeros that matters.
I think that I’m using LLR=ln(p0/p1) throughtout this code, but I’ll check this.
The reason why uncoded_gammas is omitted from the delta calculation is because this causes the algorithm to generate extrinsic LLRs directly, rather than a posterior LLRs.
Take care, Rob.
April 25th, 2012 at 3:18 pm
Thanks Rob,
So for the “gammas” to take the value of LLR when input bit is zero, we should have: P1 = 1/(1+exp(LLR)) and P0= exp(LLR)/(1+exp(LLR)) which means LLR is log(P1/P0). Because with this we’ll have log(P1) = -log(1+exp(LLR)) and log(P0) = LLR - log(1+exp(LLR)) where by omitting the common term we’ll end up with what you do in the code. Is this correct?
Another point of confusion for me is where in calculating “deltas” you’re using “bettas of bit_index” while according to text books we should use the next stage “betas” here. Shouldn’t we?
Thanks,
Elnaz
April 25th, 2012 at 3:25 pm
Hi Elnax,
Actually, P1 = 1/(1+exp(LLR)) and P0= exp(LLR)/(1+exp(LLR)) implies that LLR = log(P0/P1), not LLR = log(P1/P0).
To put the calculation into words for you, the delta of a transition is equal to the sum of
* the encoded gamma of that transition,
* the alpha of the state that the left end of the transition connects to and
* the beta of the state that the right end of the transition connects to.
This is what the text books say, although they often use confusing notation to do it.
Take care, Rob.
April 25th, 2012 at 4:29 pm
You’re right; maybe the original explanation of the Log-MAP is confusing me. I think there is only one “gammas” for each transition between two states, however, the direction of movement changes in calculating “alphas” vs. “betas”. For example, gammas of 1 changes the state from 1 to 2. In calculating bettas of 1, I think, we should do: (Bettas of 2) X (gammas of 1) not (gammas of 2) because gammas of 2 takes state 2 to state 3. But, you’re doing the gammas of 2 in the code!
Elnaz
April 25th, 2012 at 5:49 pm
Hi Elnaz,
I’m afraid I’m not sure what you mean here. You should note that the gamma and delta values correspond to particular transitions. By contrast, the alphas and betas correspond to particular states. I’d encourage you to look closer at the code - I’m quite confident that there is no mistake in it.
Take care, Rob.
April 25th, 2012 at 7:59 pm
Hi Rob,
Surely, the code is correct. I’m trying to resolve my points of confusion which is partly the way you’re indexing “gammas” in calculating “betas”. In my previous note I wrote one of the iterations in your nested loops (where you calculate “betas”). So, for example you get bettas(..,bit_index=1) from betas(…, bit_index=2) (next state) and gammas(…, bit_index=2). My confusion is that I think it should be using gammas(…, bit_index=1) instead of gammas(…, bit_index=2) because it is gammas(…, bit_index=1) that takes us from stage 1 to stage 2. gammas (…, bit_index=2) will take us from stage 2 to stage 3.
Did I succeed explaining what I don’t understand?
Thanks,
Elnaz
April 25th, 2012 at 8:15 pm
Hi Elnaz,
I see what you mean now - you are looking at the following line…
deltas(transition_index, bit_index) = alphas(transitions(transition_index,1),bit_index) + encoded_gammas(transition_index, bit_index) + betas(transitions(transition_index,2),bit_index);
Since bit_index is the same for alphas and beta, it looks like the states that correspond to these are in the same horizontal position within the trellis. If you look closely though, you’ll see that bit_index refers to a transition’s horizontal position, not a state’s horizontal position. I’ve set things up so that a particular value in alphas corresponds to a state at the left end of the transition indexed by bit_index. By contrast, betas is set up so that each of its values corresponds to a state at the right end of the transition indexed by bit_index.
Take care, Rob.
April 25th, 2012 at 8:33 pm
Hi Rob,
Actually, I’m referring to where you calculate “betas”, line:
betas(…, bit_index) = jac(betas(…, bit_index),betas(…, bit_index+1) + uncoded_gammas(…, bit_index+1) + encoded_gammas(…, bit_index+1));
I think we should have:
betas(…, bit_index) = jac(betas(…, bit_index),betas(…, bit_index+1) + uncoded_gammas(…, bit_index) + encoded_gammas(…, bit_index));
Because it is gammas(…, bit_index) that moves us horizontally from bit_index to bit_index+1 position in trellis.
According to text books:
Beta (stage=k, state=p) = sum over different states of (gamma(stage=k, state=p to state=q) * Beta(stage=k+1, state=q).
i.e. in calculating beta in stage=1(horizontal position in trellis) we use beta of the position to the right (next stage, stage=2) and gammas of the same position (same stage, stage=1) which takes us from stage=1 to stage=2. Translating this to your code, it means we should use gammas(…, bit_index) instead of gammas(…, bit_index+1).
Elnaz
April 25th, 2012 at 8:42 pm
Hi again Elnaz,
Taking my previous reply a step further, betas(…, bit_index) refers to a state that is to the left of the transition referred to by uncoded_gammas(…, bit_index+1). Similarly, betas(…, bit_index+1) refers to a state that is to the right of that transition.
Take care, Rob.
April 25th, 2012 at 9:14 pm
I see that Rob; and that is exactly my point of confusion because there is only one gammas for each transition and we use it in calculating both betas and alphas.
The way you calculate alphas exactly matches with text book but betas does not. What makes it confusing to me is that we use say gammas(…, bit_index=2) to to go from betas(…, bit_index=2) to betas(…, bit_index=1)!
I’m sorry it seems I’m repeating my question here but it seems to me that each transition is marked with one gamma only and in the code we’re referring to the same transition by different gammas when caculating alphas versus betas.
Thank you,
Elnaz
April 25th, 2012 at 9:21 pm
Hi Elnaz,
That’s right - the gammas used for the alpha calculations are the same as the ones used for the beta calculations. The only difference between alpha calculations and beta calculation is that the former uses a recursion from left to right, while the latter uses a recursion from right to left.
Take care, Rob.
April 25th, 2012 at 9:38 pm
Yes, exactly! So when we move from time k (bit_index=k) to time k+1 (bit_index=k+1) in the trellis, we’re marking this transition with gammas(…, bit_index=k). Now, here, we move to the right from alphas(…, k) to alphas(…, k+1) and move to the left from betas(k+1) to betas(k) on the same transition marked with gammas(..,k).
The question is you’re not doing the above thing in the code. you have alphas(…,k+1) and betas(…, k) referring to the same node i.e. the same point in time!
This is where I don’t get it!
Elnaz
April 25th, 2012 at 9:42 pm
Ah, but you are assuming that alphas and betas are using the same time reference - i.e. that alphas(…, bit_index) and betas(…, bit_index) are referring to the same state. Actually the state that alphas(…, bit_index) refers to is one position to the left of the state that betas(…, bit_index) refers to.
Take care, Rob.
April 26th, 2012 at 12:42 am
Thank you Rob. I get it now. It is just the way you name your betas. I was thinking your theory may be different which is not.
Actually, I want to use two BCJR decoders, parallely concatenated, to decode a rate 1/3 Turbo code. I know that I should pass the output of each decoder (function here) as the apriori_uncoded_llrs input to the other function and iterate in a figure 8 shape. But where should I use the systematic output of the channel? I can’t feed both systematic and parity streams to the decoder function.
Thank you,
Elnaz
April 26th, 2012 at 10:04 am
Hi Elnaz,
Take a look at Figure 2.11 in the nine month report. The systematic LLRs can be added to the interleaved extrinsic LLRs.
Take care, Rob.
April 30th, 2012 at 5:31 am
Hi Rob,
I learned many thing just reading the above questions and answers, so thanks for that.
I just have a small question about puncturing LDPC codes.
In my matlab script, I puncture a portion of the information bits in a systematic codeword in the transmitter side.
My problem is with the receiver side. In an earlier comment you have mentioned that I need to replace the puncture bits with zero-valued LLR values. Did you mean the LLR values from the channel to the variable node which is use to generate the apostiriori LLR values?
Even thoug I perform the above, I get a BER curve exactly as that for a LDPC code without puncturing. Really appreciate if you can explain this to me.
Thanks in advance.
April 30th, 2012 at 9:53 am
Hi Anna,
Yes, I mean the LLRs from the channel to the variable nodes. I can’t think why your BER plot stays the same when you introduce puncturing - puncturing should degrade the performance. It sounds like you think you are puncturing the bits, but actually you are still transmitting them…
Take care, Rob.
May 4th, 2012 at 4:12 pm
Elnaz says…
“Hi Rob,
I am using two BCJR decoders, parallely concatenated, to decode a rate 1/3 turbo code. As you said earlier, I am adding the systematic LLR to the output of one decoder and feed them as apriori_uncoded_llr to the other decoder.
My question is with BER plot where I add a variance \"sigma2\" noise (\"randn\" command) to my sequence and simply plot the error rate vs. (Eb=1)/(2*sigma2). However, with small sigma2 values (smaller than 1.3), the waterfall discontinues and a series of up and down fluctuations begins. Fluctuations are about BER = 10^-4. Why this is happening?
Thank you,
Elnaz”
Hi Elnaz,
It may be that you have reached the error floor of your turbo code. The position of your error floor will depend on the design of your turbo code. In particular, it will depend on your interleaver length. An error floor of 10^-4 suggests to me that your interleaver length is quite small - does this sound correct to you?
Take care, Rob.
May 4th, 2012 at 6:15 pm
Hi Rob,
Thanks for your answer.
I am interleaving the whole sequence (about 4000 bits) at once with Matlab\’s \"randintrlvd\". What do you think?
Elnaz
May 4th, 2012 at 6:39 pm
Hi Elnaz,
That sounds okay to me. I normally use Matlab’s randperm function though…
interleaver = randperm(4000);
interleaved_bits = bits(interleaver);
deinterleaved_bits(interleaver) = interleaved_bits;
Take care, Rob.
May 4th, 2012 at 6:55 pm
Yes, Rob; doing exactly like above, the fluctuations around 10^-3 remain to be. They are happening with SNRs higher than -5 dB.
May 4th, 2012 at 7:03 pm
Hi again Elnaz,
My advice would be to try using a 10-times longer interleaver length, to see how that affects things.
Take care, Rob.
May 4th, 2012 at 7:11 pm
Hi Rob,
Can it be longer than the whole data length? Currently I am using
interleaver = randperm(4000); which covers my whole data length (entire stream of original bits).
Do you mean to somehow zero-pad the data?
Thanks,
Elnaz
May 7th, 2012 at 6:18 am
Dear Rob,
I was studying the paper \"Analysis of low-density parity-check codes based on EXIT functions\" by Sharon, E. \"Communications, IEEE Transactions on\". It was describing the EXIT trajectories in very good ways.
I tried to implement the EXIT trajectory using the update rules in equations (8,9,10), but I couldn\’t. Is anyone can implement these trajectory rules update for me please?
In this way, it is easy to calculate the trajectories. Thanks
May 7th, 2012 at 11:34 am
Hi Elnaz,
I would suggest to repeat your 4000-bit message ten times, to create a 40000-bit message.
Take care, Rob.
May 7th, 2012 at 11:49 am
Hi Ideal,
It is Equations 11 and 12 from that paper that you should use to draw the EXIT functions and trajectories.
Take care, Rob.
May 7th, 2012 at 12:00 pm
Dear Rob,
You are right, Equation 11 and 12, we can get the extrinsic mutual information for the variable node using 11 and we can get the extrinsic mutual information for the check node using equation 12.
But how can we plot the trajectory curve (stair curve). I don’t understand the trajectory curve. I can plot the variable node curve and check node curve but can’t able to plot the trajectory.
Secondly, Is the trajectory curve plot will also be same with the puncturing scheme? Thanks, I am asking so many questions, but I learn a lot with your supervision. Thanks
May 7th, 2012 at 1:34 pm
Hi Ideal,
You need to use Equation 11 and 12 alternately. You start with an input MI of zero, then use the output MI of each equation as the input to the other. This process is the same whether you use puncturing or not.
Tae care, Rob.
May 7th, 2012 at 3:37 pm
Hi Rob,
Even with increasing the length to 40000, it still jumps around 10^-4, not much of an improvement I guess. Could you guide me as to what is the standard basis of BER floor in a 1/3-rate Turbo decoder?
I don\\\’t know if 10^-4 is normal or there is something wrong with my code.
Thanks,
Elnaz
May 7th, 2012 at 5:28 pm
Hi Elnaz,
It sounds to me like you have a problem with a few bits in each frame. This suggests to me that there may be a problem with the way you are terminating the trellis. Perhaps you are not using termination in the encoder, but your decoder assumes otherwise…
Take care, Rob.
May 7th, 2012 at 5:49 pm
Hi Rob,
Actually, for the encoding part, I am the same encoder twice, once with systematic output and once without that:
trellis1 = poly2trellis(3,[7 5],7);
code1 = convenc(msg,trellis1);
Then I interleave:
intrlvd = msg(interleaver);
trellis2 = poly2trellis(3,5,7);
code2 = convenc(intrlvd,trellis2);
Using convenc, do I have control over termination?
After this I am using your BCJR decoder in a turbo-decoder and decode the whole stream at once (not frame-by-frame).
Could you please explain more as to what you mean by \\"decoder assumes otherwise…\\"?
Thanks,
Elnaz
May 7th, 2012 at 8:23 pm
Hi Rob,
Let me add this: in high SNRs I am having, as you said, a few bits of error.
Also, in my turbo decoder, I initialize the extrinsic LLR of the 2nd decoder with zeros before iterations start.
Does this mean that the encoders should be in state zero at the end of the stream? If yes, how can I force them to zero state while using \"convenc\". I appreciate if you could clarify this for me more.
Thanks,
Elnaz
May 7th, 2012 at 9:29 pm
Hi Elnaz,
I think my suspicion is correct. I don’t think that convenc uses termination. You can stop my component_decoder.m from using termination by changing the line…
betas(1,length(apriori_uncoded_llrs))=0;
…to…
betas(:,length(apriori_uncoded_llrs))=0;
As you say, termination forces the encoder into the state zero at the end of the stream. If you want to achieve this using convenc, I think you will need to add some additional bits to the end of your message, as in the UMTS turbo encoder.
Take care, Rob.
May 8th, 2012 at 5:18 am
Dear Rob,
Thanks for reply, I have one more question please.
In the capacity approaching codes, we usually compare the SNR_threshold of LDPC code with the corresponding channel capacity.
I found the channel capacity equation is the J function, i.e., C = J(sigma_ch), where the sigma_ch = sqrt(8 * R * Eb/No).
my question is that the equation I wrote is OK to compare with the SNR_threshold, if yes, then what should be the value of R? It will be the mother code rate or not? Thanks
May 8th, 2012 at 5:34 am
Hi Rob,
Changing the code, as you said, nulled out the error at high SNRs; Now I have zero errors at high SNRs. But BER still doesn’t fall smoothly. For example, with #iterations=20, it sharply falls from 10^-2 to -inf. And, with #iterations=2 it fluctuates between 10^-4 and -inf (step size is very small). I never see 10^-6, 10^-7 or…
Let me add that, here, since I just want to test the decoder, I only encode, then add noise, then decode. I’m thinking maybe this is causing that since I’m not doing any modulation and all.
Do you think the decoder is what it should be? Is not using termination causing this?
Thanks a lot,
Elnaz
May 8th, 2012 at 9:31 am
Hi Ideal,
I’m not sure where you have got that relationship between the channel capacity and J function from - I’m not sure that you can use it in this context. Presumably, you are trying to quantify the BPSK DCMC capacity. I would suggest doing this using the code at…
http://users.ecs.soton.ac.uk/rm/resources/matlabcapacity/
Take care, Rob.
May 8th, 2012 at 9:33 am
Hi Elnaz,
What you are describing sounds like a steep turbo cliff, which is usual for turbo codes. I think that the reason why you cannot observe low BERs is because your simulation is not long enough. To observe a BER of 10^-7, you need to be transmitting 10^9 or 10^10 bits, which requires a very long simulation (perhaps taking one week to run).
Take care, Rob.
May 8th, 2012 at 11:59 am
Dear Rob,
I got it the plot of channel capacity vs the coding threshold. Thanks again.
I am a bit poor in writing tricky MATLAB code. Like about the EXIT trajectory, You gave the idea that the two equations are used alternatively and in this way the information stairs end up at (1,1) but when I start writing the code, I am facing problem. Can you please give me a bit clue please.
Thanks
May 8th, 2012 at 12:25 pm
Hi Ideal,
Something like this should do the trick…
IA = 0;
IE = 0;
results = [IA,IE];
for iteration_index = 1:100
IE = variable_node(IA);
results = [results; [IA,IE]];
IA = check_node(IE);
results = [results; [IA,IE]];
end
plot(results(:,1), results(:,2);
May 16th, 2012 at 7:52 am
dear rob,
i want to test an image transmission over wireless sensor network using turbo codes can this be done using the code …?
THANK YOU
May 16th, 2012 at 8:28 am
Hi Sowmya,
This depends on if your wireless sensor network can run Matlab code. If not, you can build a turbo code in the C programming language using the code at…
http://users.ecs.soton.ac.uk/rm/resources/cfixedbcjr/
Take care, Rob.
May 16th, 2012 at 3:40 pm
Hi Rob,
Actually, I have a question on iterative timing recovery and decoding. Since, you helped me alot previously with BCJR decoder, I thought I try to ask this one as well.
I’m using GD (Gradient Decent algo) to estimate timing points and a Turbo decoder to estimate user bits. My question is that at each iteration when I receive the estimated user bits, how should I use the output of the channel to update the timing estimate?
Since the channel output is constructed using coded bits, I can’t use the same channel output in my GD, unless, I code the decoded bits of the Turbo decoder before feeding it to GD. Am I correct in assuming this? How does this generally work?
Thanks a lot,
Elnaz
May 16th, 2012 at 6:57 pm
Dear Rob,
In answer to my last comment, I think, I need to actually change BCJR blocks in my Turbo decoder to make decisions over coded bits as opposed to user bits. I know I should change the sums to sum over different terms in alfa, gamma, betta but I’m not sure how exactly. Could you please guide me how to change your code to get there?
Thanks,
Elnaz
May 17th, 2012 at 9:57 am
Hi Elnaz,
You are quite right - you need the BCJR decoders to generate extrinsic LLRs (or perhaps a posteriori LLRs would suit your application better) for the coded bits. You can do this by creating an additional set of deltas using the same for loops, but with the line…
deltas2(transition_index, bit_index) = alphas(transitions(transition_index,1),bit_index) + uncoded_gammas(transition_index, bit_index) + betas(transitions(transition_index,2),bit_index);
Then, you can create an additional set of extrinsic LLRs using the same for loops, with with the line…
if transitions(transition_index,4)==0
If you decide that you want the a posteriori LLRs, you can obtain them by adding the extrinsic LLRs to the a priori LLRs.
Take care, Rob.
May 17th, 2012 at 3:17 pm
Hi Rob,
Don’t I need to modify alfas, or betas, or anything else in BCJR?
Now, I expect that the extrinsic outputs that I get from my Turbo decoder to be equal to parity1 and parity2 (encoded bits) but they’re not.
Thanks,
Elnaz
May 17th, 2012 at 4:46 pm
I’m not sure I understand why replacing encoded_gammas with uncoded_gammas will give me extrinsic LLRs for the coded bits. Could you please explain more.
Thanks,
Elnaz
May 17th, 2012 at 5:14 pm
Hi Elnaz,
There is no need to modify the alphas, betas or gammas. In fact, you can see some code that gives extrinsic LLRs for the coded bits at…
http://users.ecs.soton.ac.uk/rm/resources/matlabexit/
The reason why replacing encoded_gammas with uncoded_gammas will give you extrinsic LLRs for the coded bits is because extrinsic information is obtained by considering all available information, except for the corresponding a priori information. So, we leave out encoded_gammas when we want encoded extrinsic information, and we leave out uncoded_gammas when we want uncoded extrinsic information.
Take care, Rob.
May 17th, 2012 at 5:27 pm
Hi Rob, Am I correct in assuming that sign of encoded extrinsic information (the output of say 1st BCJR decoder) after completing the iterations of Turbo decoder should be equal to parity1 i.e. the coded bits resulted from the first encoder in my 1/3 turbo decoder?
Currently, I’m not changing the tturbo decoder at all but only replacing the new BCJR blocks. Also, as you advised me before, I am adding systematic output to the extrinsic LLR of the second decoder and feed them as apriori_uncoded_llrs to the first BCJR decoder.
If all above is right, then whay am I not getting the same parity1 coded bits from the turbo decoder?
Thanks,
Elnaz
May 17th, 2012 at 5:30 pm
Hmm, I’m not sure. You should still be iteratively exchanging extrinsic LLRs pertaining to the uncoded bits. So you need to output extrinsic LLRs for both the uncoded and encoded bits. You should find that the sign of the LLRs correlates with the encoded bits, as you say. My recommendation would be to compare your code with that which I linked to in my previous reply…
Take care, Rob.
May 17th, 2012 at 5:43 pm
Oh, in that case, I was doing wrong. I have basically replaced the old BCJR with new BCJR that outputs extrinsic llrs which is made from delta2. So, turbo decoder was iteratively exchanging new extrinsic llrs.
But you are saying that turbo decoder should still be exchanging the old outputs (extrinsic llrs for uncoded bits). Here is where I’m not clear. Suppose I get both outputs old and new from each BCJR block, If the turbo decoder should still be exchanging old outputs, then what I should do with the new ones? exchange them too? what does the schematic of the new turbo decoder look like?
Thanks a lot,
Elnaz
May 17th, 2012 at 6:18 pm
No, you should only exchange the uncoded extrinsic LLRs. You should not use the encoded extrinsic LLRs until the iterative decoding process has finished. After this, you can provide these LLRs to your synchroniser.
Take care, Rob.
May 21st, 2012 at 9:16 am
Dear Rob,
Thanks for your thorough help. You suggest me the code for trajectory plot. I tried the following code. It works for variable and check node curve but doesn’t work for trajectory.
clear all; clc;
IA = 0;
IE = 0;
results = [IA,IE];
load newJ.mat; % loading J-function
Es = 2;
dv = 3; dc = 6;R=0.5;Ite_max=100;
snr_L = 10^(Es*0.1);
sig_Ch = sqrt(8 * R * snr_L);
% Variable Node curve
Iav=0:.1:1;
sigmaA = (interp1(J,SIGMA,Iav));
sigmaE = sqrt((dv - 1).*sigmaA.^2 + sig_Ch^2);
Iev=interp1(SIGMA,J,sigmaE);
Iev(isnan(Iev)) = 1;
% CNode EXIT curve
Iac=0:.1:1;
sigmaCn1 = interp1(J,SIGMA,(1-Iac));
sigmaEC = sqrt(dc-1) * sigmaCn1;
Iec = 1 - interp1(SIGMA,J,sigmaEC);
Iec(isnan(Iec)) = 0;
plot(Iav,Iev,’b-’,'LineWidth’,1.6),hold on,grid on
plot(Iec,Iac,’r-’,'LineWidth’,1.6)
% Trajectory
while iteration_index == Ite_max
sig_A = interp1(J,SIGMA,IA);
IE = interp1(SIGMA,J,sig_A);
results = [results; [IA,IE]];
sig_C = interp1(J,SIGMA,IE);
IA = interp1(SIGMA,J,sig_C);
results = [results; [IA,IE]];
end
plot(results(:,1), results(:,2));
I am really messed with this program.
May 21st, 2012 at 9:31 am
Hi Ideal,
It looks to me like your while loop for the trajectory will never be performed. You should replace it with a for loop, such as…
for iteration_index = 1:100
Take care, Rob.
May 21st, 2012 at 11:26 am
Dear Rob,
I am really very happy. I done that trajectory. Now its quite good. Thanks for your support. Thanks
May 22nd, 2012 at 11:28 am
Dear Rob,
Thanks for all the responses. One more question please
How can we find the area between the two EXIT curves? I want to find the minimum positive area (distance) when the tunnel open.
Is there any script for finding the minimum area or distance between the two curves please? Thanks
May 22nd, 2012 at 12:33 pm
Hi Ideal,
The code that I have posted at http://users.ecs.soton.ac.uk/rm/resources/matlabexit/ includes a calculation of the area beneath an EXIT curve. If you know the area beneath each EXIT curve, then you can subtract one from the other to find the area within the EXIT tunnel.
Take care, Rob.
May 23rd, 2012 at 7:20 am
Dear Rob,
I implemented that method. I was thinking that it will give negative area when the tunnel closes (Area = variable node curve area-check node curve area), but it still gives me the positive area for closed tunnel.
Can I find the minimum positive distance between the two curves. I thought may it will work but don’t know how to find the minimum positive distance between the two curves?
Thanks
May 23rd, 2012 at 8:03 am
Hi Ideal,
A greater variable node area than check node area does not guarantee an open tunnel. This is the reason why some LDPC codes cannot operate close to the channel capacity. To find the smallest gap between the curves, you’ll have to use interpolation to convert one of the EXIT curves into the domain of the other. You can then subtract each data point in one curve from the corresponding point on the other curve.
Take care, Rob.
May 25th, 2012 at 5:47 am
Dear Rob,
I got some of your idea but not the full one. can you please elaborate a bit more. How can we find the minimum positive distance between the two curve? I really need it. if possible please provide me some example script. Thanks a lot
May 25th, 2012 at 11:31 am
Hi Ideal,
Here is an example for you…
IA_outer = 0:0.1:1;
IE_outer = (0:0.1:1).^2;
IA_inner = 0:0.1:1;
IE_inner = 0.5:0.01:0.6;
IE_outer_new = IA_inner;
IA_outer_new = interp1(IE_outer, IA_outer, IE_outer_new);
plot(IA_inner, IE_inner, ‘-’, IE_outer_new, IA_outer_new, ‘–’)
Take care, Rob.
IE_inner - IA_outer_new
May 29th, 2012 at 12:47 pm
Dear Rob,
Thanks for reply,
I have one simple question. Usually, we implement puncturing in the codeword, for example something like this
Transmitter side (Uncoded Case)
temp = randperm(N); randomly permuted vector of size N
Np = N - N*p; length after puncturing with puncturing fraction p
pattern = temp(1:Np); Assignment of pattern
bits = rand(1,N)<0.5; Random binary vector
% Puncturing
punc_bits = bits(pattern);
% Calaculating LLR
tx = -2*(punc_bits - 0.5); (0->1, 1->-1)
y = tx + n;
Lch = 4 * SNR_Linear * y;
In the receiver side
llrs = zeros(size(bits));
llrs(pattern) = Lch;
My question is what about the coded bits, something like below?
code_bits = [p data]; (p=parity bits) (codeword)
tx_bits = -2*(code_bits - 0.5); (0->1, 1->-1)
% Puncturing
1. c_punc_bits = tx_bits(pattern); ?
OR
2. c_punc_bits = code_bits(pattern);
tx_bits = -2*(c_punc_bits - 0.5); (0->1, 1->-1)
Which of the above two way for coded system is right? Thanks
Please comment
May 29th, 2012 at 12:54 pm
Hi Ideal,
Those two ways are equivalent - you can use either.
Take care, Rob.
June 4th, 2012 at 12:29 pm
Dear Rob,
one quick question please. I have defined puncturing fraction p. now i want to distribute it into four different ratios whose overall results should be equal to p.
i.e. p = p1 + p2 + p3 + p4. how to define four variables based on one variable over a range of values for each variable. Not sure right question:(
Thanks
June 4th, 2012 at 2:33 pm
Hi Rob,
I have another question that is not related to the decoding here.
I want to implement a convolution of a stream of bits (with zero-padded intervals in between to avoid ISI) with a delayed pulse shape. However, delay is a variable of bit index; I want to implement something like this:
Sigma over k of (a(k) * g(t-kT-tau(k)),
where tau(k+1) = tau(k) + noise + some offset;
Do you know how to implement this?
Thanks,
Elnaz
June 6th, 2012 at 6:12 pm
Hi Ideal,
It seems to me that if you want p to be composed of four components, then it should be obtained as their product, not their sum…
p = p1*p2*p3*p4
This is because the overall coding rate of serially concatenated codes is given by the product of the individual coding rates.
Take care, Rob.
June 6th, 2012 at 6:15 pm
Hi Elnaz,
I’m afraid that I’ve never done this before, so I’m not sure how to do it. You may like to look at Matlab’s convolution functions, although I’m not sure that it will let you use a variety of delays…
Take care, Rob.
June 7th, 2012 at 5:09 am
I am running a simulation of umts turbo codes, with frame length 380, 8 iterations, using max-log-map , but BER results doesn\’t go below 10^-6.
I am using C code. I want this to go down to 10^-8. any suggestions
June 7th, 2012 at 9:04 am
Hi Nauman,
If you want your BER results to go down to 10^-6, then you will need to simulate the transmission of at least 10^8 bits. With your frame length, this corresponds to 264k frames…
Take care, Rob.
June 8th, 2012 at 7:39 am
Hi Rob,
I have sent these number of bits, I am only concerned about the error floor of
UMTS turbo codes, and can I achieve 10^-8 BER with the UMTS interleaver.
June 8th, 2012 at 8:37 am
Hi Nauman,
As you can see in my results above, the error floor of the UMTS turbo code is very steep. It may be that you can’t tell the error floor apart from the turbo cliff…
Take care, Rob.
June 11th, 2012 at 6:20 pm
Hi Rob,
I have one question about adding AWGN:
I need to plot my error versus (Eb/N0 (dB)) and I have about 4000 bit which after coding with 1/3 rate turbo code, I arrange them into a 111X111 matrix.
I calculate the noise variance as:
sigma2 = 1/(2*(10^(EbN0/10)));
r = R + sqrt(sigma2).*randn(size(R,1),size(R,2));
But I’m having second thoughts as if I have to somehow enter the message length 4000X3 bits or 4000 bits into the above equations. Should I?
Thanks,
Elnaz
June 12th, 2012 at 8:15 am
Hi Elnaz,
You should be adding the noise to the 12000 encoded bits. One thing to consider is that you are using purely real valued noise here. If you are using BPSK then this is okay. But otherwise, you should be using complex valued noise…
r = R + sqrt(sigma2)*(randn(size(R))+i*randn(size(R)));
Take care, Rob.
June 12th, 2012 at 4:12 pm
Hi Rob,
Yes, I am adding the noise to all the samples. The problem is with power of noise that I\’m adding. I\’m using a 1/3-rate turbo code and then I multiply my coded bits each with a 2D sinc function and then adding the noise.
Does the above formula: sigma2 = 1/(2*(10^(EbN0/10))) still hold here?
Thanks,
Elnaz
June 12th, 2012 at 4:18 pm
Hi Elnaz,
Your equation sigma2 = 1/(2*(10^(EbN0/10))) can be rearranged to give N0 = 2*sigma^2. This implies that your noise will have two dimensions (i.e. it will be complex). If you were using purely real-valued noise, then N0 = sigma^2.
Take care, Rob.
June 12th, 2012 at 5:17 pm
Shouldn’t I use something in line with: Eb/N0 = Fs/(2 x R x sigma2) ?
where
R = code rate 1/3 in my case
Fs = No. of samples per modulated symbol. Since my my system is two-dimensional, I guess I will have to take the no. of samples in my sinc(x).sinc(y) function.
Since I’m adding pure real-valued noise then: Eb/N0 = Fs/(R x sigma2) .
Actually, I think my error vs. EbN0 plot are shifted to the right compared to the correct plots as if I’m not adding enough noise power. That’s the reason I’m trying to increase the noise power somehow.
Thanks,
Elnaz
June 13th, 2012 at 3:31 pm
Hi Elnaz,
When I perform a baseband simulation, I use SNR = Es/N0, where Es is the energy per symbol. I then convert to Eb/N0 by using Es = R*log2(M)*Eb, giving Eb/N0 = Es/(N0*R*log2(M)). Here M is the number of constellation points in the modulation scheme.
Take care, Rob.
June 14th, 2012 at 9:59 pm
Hi Rob,
Yes, right. I depends, I think, where we’re adding the noise in the system. What I’m doing is somewhat different than a standard modulation and that is the reason I have difficulties applying the same Eb/N0 concept to my problem.
Actually, I am rearranging my stream of coded bits into a matrix and then convolve this matrix with a 2D sinc function ending up with an image (2D waveform). Therefore, I have an image wherein the coded bits are hidden. I add my noise to this matrix.
I want to know how exactly I should calculate the noise variance in terms of Eb/N0.
Any reference or suggestion would be great.
Thanks,
Elnaz
June 15th, 2012 at 5:04 pm
Hi Elnaz,
I’m afraid that I haven’t come across the approach you are describing here so I’m afraid that I can’t help you with the relationship between Eb/N0 and the noise variance. My suspicion is that N0 = 2sigma^2, since you are probably using complex noise. The bit I’m not sure about is the value of Eb, or Es…
Take care, Rob.
June 20th, 2012 at 11:05 am
Dear Rob,
In your post above (May 25th, 2012 at 11:31 am), you suggest me interpolation to convert one of the EXIT curves into the domain of the other. In this way we can then subtract each data point in one curve from the corresponding point on the other curve. This way can find the smallest gap between the curves.
I already implement this method and it is working fine but now I want to write these step of interpolation of two curves in a mathematical way but feeling difficulty. Please can you give me example of writing this?
June 20th, 2012 at 3:18 pm
Hi Ideal,
Suppose you have two data points (x1,y1) and (x2,y2). You wish to find the y-coordinate y3 that corresponds to a particular x-coordinate x3. You can interpolate along a line passing through the two known data points using the following steps…
m = (y2-y1)/(x2-x1)
c = y1 - m*x1
y3 = m*x3+c
Take care, Rob.
June 20th, 2012 at 6:28 pm
the exit chart graph is plotted in between Eb/No vs BER.i got the out put as after10 trial and 500 message bit passed i got the waveform.can i calculate the code rate?
June 21st, 2012 at 7:32 am
Hi Somya,
The code rate of a UMTS turbo code is given by N/(3*N+12), where N is the interleaver length. In the case of UMTS, the interleaver length N should in the range 40 to 5114 bits.
Take care, Rob.
June 25th, 2012 at 12:24 pm
Hello Rob,
I know you very competent in turbo codes and MIMO systems. I read some papers from L.Hanzo and you.
I work on turbo-MIMO matlab code: spatial multiplexing MIMO with MAP (or ZF, MMSE) and turbo-code (PCCC with MAP).
Could you please to suggest me any help in this field?
June 25th, 2012 at 2:19 pm
Hello Stas,
I’m afraid that I don’t have any Matlab code for spatial multiplexing MIMO, but the Matlab code that I have provided on this page can help you with your turbo-code - it is a PCCC with MAP.
Take care, Rob.
June 27th, 2012 at 1:27 am
Dear Rob,
Thanks for your help. Please one more question.
Do you have any script for plotting “zoom out a portion of plot within the same plot”?
Thanks
June 27th, 2012 at 9:24 am
Hi Ideal,
I’m afraid that I don’t, because I do all my final-version plotting in Gnuplot, rather than Matlab. I’m sure you will be able to find some code for this on Google though.
Take care, Rob.
June 29th, 2012 at 6:08 am
Dear Rob,
Yes I did that plot very well thanks
July 4th, 2012 at 5:00 am
Dear Rob,
You explain me last time the way of puncturing in MATLAB. I did that one and working very good. It was the way for binary channels.
But Now I am think if I am using QPSK or 4-PAM, in this case, each symbol consist of two bits. I am thinking to apply puncturing (pairwise bits), i.e., puncturing symbol treating all the bits uniformly but I don’t know how to implement pairwise puncturing. Please can you give me any idea, how to implement it. Thanks
July 4th, 2012 at 9:34 am
Hi Ideal,
My feeling is that puncturing bits in pairs places an unnecessary constraint on your design and that it may degrade your performance. My recommendation would be to puncture the bits in the way you have already been doing it, then pair the remaining bits up and modulate them.
Take care, Rob.
July 4th, 2012 at 12:02 pm
Dear Rob,
Thanks for your advise. But after puncturing, I need to calculate the LLR for each bit, i.e., after puncturing I need to perform QPSK modulation and then calculate the LLR(x1) and LLR(x2) in a symbol and then calculate the EXIT chart? I think your idea is very useful but..
I am still confuse please can you explain bit more. Thanks
July 4th, 2012 at 12:51 pm
Hi Ideal,
You just need get the LLRs from the demodulator, then insert a zero-valued LLR to replace each punctured bit. You should think of the puncturer as being independent from the modulator.
Take care, Rob.
July 5th, 2012 at 12:17 pm
Dear Rob,
I am thinking of first puncturing the codebits and then modulate them into QPSK and then transmit. On the receiver side, I will receive y=x + n, x is the QPSK symbol and n is the noise. Now LLR(y) = 2/sigma *y is the right thing? or I need to calculate the LLR for real bits and then for imaginary bits and combine the LLR into one LLR.
Then I can add the same pattern of punctured, I did on the transmitter side.
Am I right or any sample example will be very helpful please. I want to calculate the EXIT chart of the LDPC code. Thanks
July 5th, 2012 at 6:30 pm
Hi Ideal,
If you use Gray coded QPSK then you can use…
LLR1 = sqrt(2)*real(y)/sigma^2;
LLR2 = sqrt(2)*imag(y)/sigma^2;
I think…
Take care, Rob.
July 6th, 2012 at 12:41 am
Dear Rob,
Thanks but I am still confuse. Like I punctured some codebits, then I modulate using the Gray mapping QPSK, now each bit is converted into a symbol (two bits), Then similarly I calculate the bitwise LLR for real bit 1 and imaginary bit 2. It means that I have two LLRs on the receiver side, LLR1 and LLR2. How can I combine two LLRs into one LLR so that I can assign zeros to the punctured position(I did at the transmitter).
I am still not understanding the cooprations between bits and symbols? After all to calculate the extrinsic mutual information, I need one LLR vector (either bit 1 or bit 2 or combined) that is the question? Thanks for your help please
July 6th, 2012 at 3:26 pm
Hi Ideal,
There is no need for you to “combine two LLRs into one LLR”. Suppose you have 100 bits, then you puncture 40 of them to leave 60 bits. After that you combine two bits per QPSK symbol, to give 30 symbols in total. In the receiver, you can use the code I provided above to separate each QPSK symbol in to two LLRs. So you go from 30 symbols to 60 LLRs. After that you can insert 40 zero-valued LLRs to get you back to 100 LLRs.
Take care, Rob.
July 7th, 2012 at 1:22 pm
Dear Rob,
I am really very thankful to you for your detailed example. So it means that LLR1 is of 30 bits and LLR2 is of 30 bits too and total of 60 bits so LLR = [LLR1, LLR2]; to have LLR equal to 60 bits.
Secondly I need to use QPSK gray mapping to convert from 60 bits to 30 symbols. Can I use your script for QPSK gray mapping scheme, you provided for QPSK demapper exit chart?
I am very much thankful to you.
July 8th, 2012 at 8:36 am
Hi Ideal,
That’s right. You can use the QPSK code I have put on my website for this job.
Take care, Rob.
July 9th, 2012 at 4:10 am
Dear Rob,
Thanks, I implement all the steps and now my algorithm is working. But I have one question about the LLR, You wrote that
If I use Gray coded QPSK then I can use…
LLR1 = sqrt(2) * real(y) / sigma^2;
My question is why you write sqrt(2), I was thinking, it will be just 2 without square root, i.e..
LLR1 = 2 * real(y) / sigma^2;
Please can you clarify me or any reference will be Fine as well. Thanks
July 9th, 2012 at 8:46 am
Hi Ideal,
This is because QPSK puts half of its power into the I signal and the other half into the Q signal. By contrast, BPSK puts all of its power into the I signal. This is why different coefficients are needed.
Take care, Rob.
July 9th, 2012 at 5:11 pm
Hi Rob,
As you know, I am using a Turbo Decoder which includes your BCJR algo in my application. However, I am getting an inferior result which is traced down to my Turbo Decoder. I have a number of uncertainties which may be the culprit.
First is the format you defined the Trellis structure in your BCJR code. I have translated my Trellis (i.e. poly2trellis(3,[7 5],7); in MATLAB) aacording to your structure. But, I feel, I can never be 100% sure unless check it with you. So, could you please have a look at this and conform its correctness:
% FromState, ToState, UncodedBit, EncodedBit
transitions = [1, 1, 0, 0;
2, 3, 0, 0;
3, 4, 0, 1;
4, 2, 0, 1;
1, 3, 1, 1;
2, 1, 1, 1;
3, 2, 1, 0;
4, 4, 1, 0];
July 9th, 2012 at 6:00 pm
Hi Elnaz,
Actually, there is a simple way to know if you have programmed the turbo code correctly or not. You should start by drawing the turbo code EXIT chart by using the histogram method. After that, you should do it again, but using the averaging method. If there are no problems with your code, you will get the same result with both methods. If there are problems, the two EXIT charts will look very different to each other.
Take care, Rob.
July 9th, 2012 at 6:08 pm
Ok; meanwhile, I’m also suspicious if using convenc which doesn’t force termination may be the culprit also. So, I need to add zeros to the end of the input bits before convenc, I guess. The question is how many zeros for poly2trellis(3,[7 5],7) to ensure termination?
Thanks,
Elnaz
July 9th, 2012 at 6:45 pm
Hi Rob,
I need to correct my previous comment by changing the question to: What bits should I add to my input stream to ensure encoder termination? Is there any way to find out other than manually drawing the trellis?
Thanks,
Elnaz
July 9th, 2012 at 7:21 pm
Hi Elnaz,
Since your encoder includes two memory elements, you will need two extra bits to terminate it. These bits should be set equal to the feedback bits. If you take a look at the UMTS encoder schematic, you will see what I mean…
Take care, Rob.
July 9th, 2012 at 7:40 pm
Hi Rob,
Do you think that not using termination can contribute to noticable error when 4000 message bits are encoded with a 1/3-rate Turbo code?
Elnaz
July 10th, 2012 at 12:43 am
Hi Rob,
Thanks for your comments. I have worked out the things and well with gray mapped QPSK. Now please can you tell me a bit about the anti-gray mapped QPSK.
I mean, how to mapp from bits to symbols (ant-gray QPSK) and what will be the LLR expressions. Thanks a lot for your thorough help.
July 10th, 2012 at 4:20 am
Hi Rob,
I studied the parts in the report which deals with termination and also the way you have implemented it.
Actually, since my application is different than yours, I can not exactly do what you did. I am using a 1/3 rate turbo encoder. Therefore, I will be having two tail bits padded to the systematic output, and the parity 1 and parity 2 will each also be two bits longer than the original stream. These 3 stream of bits then pass through a channel. The discrepancy between my application and your application is that in the receiver site , I won’t have access to the termination bits of the second encoder to feed the second BCJR. But for the 1st BCJR there is no problem since termination bits of the 1st encoder is already padded to the systematic output and passed through the channel. So, I will have the channeled version of the termination bits of the 1st encoder but not the second one.
What would you suggest?
Thanks,
Elnaz
July 10th, 2012 at 7:54 am
Hi Elnaz,
Getting the termination wrong can give you a few bit errors per frame. So this would give you an error floor at about 10^-3, using a frame length of 4000 bits.
The UMTS turbo code solves the problem you’ve described by transmitting the bits that are padded onto the interleaved bit sequence in order to terminate the second encoder. The alternative is to pad the apriori LLRs with two zero-valued LLRs.
Take care, Rob.
July 10th, 2012 at 7:56 am
Hi Ideal,
If you use anti-Gray mapping (also known as natural mapping) then you can’t use the simple demodulator equations any more. You can see how to implement this using the code that I provided at…
http://users.ecs.soton.ac.uk/rm/resources/matlabexit/#comment-1545
Take care, Rob.
July 10th, 2012 at 2:59 pm
Hi Rob,
I implemented the alternative i.e. padding the apriori LLRs with two zeros. And my apriori encoded LLRs already include two encoded and channeled termination bits (output of component_encoder are passed through channel).
I am using the encoded outputs of BCJR decoder at each iteration.
Is it correct to use their entire length or the last two bits are erroneous?
Thanks,
Elnaz
July 11th, 2012 at 6:45 am
Hi Elnaz,
All of the LLRs should be good to use. The last tow will pertain to the two termination bits.
Take care, Rob.
July 11th, 2012 at 11:21 am
Dear Rob,
Thanks for QPSK demapper exit chart code. Can you advise me a bit its implementation for the LDPC EXIT chart, i.e., I want to simulation for Variable node and check node (3,6) half rate. In that case, we will have two a-priori LLRs from check nodes and one channel LLR in order to calculate the variable node extrinsic information. Can you advise me about its implementation like this. please
July 11th, 2012 at 3:37 pm
Hi Ideal,
As you say, each variable node takes LLRs from it’s connected check nodes and one LLR from the QPSK demodulator. If you use natural QPSK mapping, there is a benefit to iterating between the demodulator and variable nodes, in addition to between the check and variable nodes. In which case, each variable node can generate one extrinsic LLR to give to the demodulator.
Take care, Rob.
July 11th, 2012 at 5:44 pm
Hi Rob,
I don\’t see any performance gain after adding termination. I changed BCJR decoder line from
betas(:,length(apriori_uncoded_llrs))=0;
back to ….
betas(1,length(apriori_uncoded_llrs))=0;
And, I pad two zeros to apriori uncoded LLRs before feeding them to the decoders.
Also, previously I changed your BCJR code to also output encoded LLRs and I am only using those.
Is there anything I might be missing?
Thanks,
Elnaz
July 12th, 2012 at 3:06 am
Dear Rob,
I want to know, is there a way to test turbo decoder performance independant of modulation and channel scheme?
My observation is that in my application the turbo decoder performs inferior than what is expected. The EXIT charts that you advised me to draw, I think, test the encoding, modulation, and one decoder but not the turbo decoder.
What would you suggest?
Many thanks,
Elnaz
July 12th, 2012 at 5:07 am
Dear Rob,
Thanks for reply,
I think still I can’t understand the operation of iterations between the demodulator and variable node, and also between the variable node and check node. May be I don’t still have a clear overview of these iterations. Any clarification, how these iterations will work in simulation setup will be appreciated please.
July 12th, 2012 at 8:27 am
Hi Ideal,
You can see what I mean in Figure 2 of this paper…
http://gps-tsc.upc.es/comm2/publications/C_2006_ISTC06_Lahuerta_1274.pdf
Here you can see that LLRs are flowing from demodulator to variable node, but also from variable node to demodulator.
Take care, Rob.
July 12th, 2012 at 8:31 am
Hi Elnaz,
Termination doesn’t make any significant difference to the performance in the turbo cliff region of the BER plot. It only really makes a difference in the error floor region. Probably your error floor is at such a low BER that you can’t see it (and its improvement) in your BER plots.
You can test that the turbo decoder is operating property by plotting decoding trajectories into the EXIT chart. If there is not a good match between the trajectories and the EXIT functions, then you know something is going wrong.
Take care, Rob.
July 12th, 2012 at 1:56 pm
Hi Rob,
I plotted the EXIT charts and the two are different. The one with histogram method is more square-like but the other is more elliptical.
I know this means there is something wrong but I don’t know how should I further debug the problem.
All I did is that I replaced your encoding and modulation with mine. I don’t have demodulator and my frame-length is about 4000 bits and I used termination but Eb/N0 value was -4 same as yours.
Can you help me find out what is wrong here? I don’t exactly understand what these charts are measuring.
Thanks,
Elnaz
July 12th, 2012 at 2:33 pm
Hi Elnaz,
The first thing I would do is try to narrow down the problem to either the turbo code or to the demodulator. You can do this by seeing how the mutual information of the LLRs provided by the demodulator vary with SNR. You can plot this using both the averaging and histogram methods. If your demodulator is working, then these plots will match each other, suggesting that the problem is in the turbo code. If the plots don’t match then it suggests that the problem is in your demodulator.
Take care, Rob.
July 12th, 2012 at 3:02 pm
Hi Rob,
But I do not have demodulator only modulator. Did you mean modulator?
Without demodulator should I still expect to get similar plots (histogram and averaging) with different SNRs?
By the way, I ran your code to see what should I expect with my code. I’m not even sure if those plots I’m getting are correct. They are not exactly matching either. Is there any way I can send you the figures from your code?
Thank you so much,
Elnaz
July 12th, 2012 at 3:25 pm
Hi Elnaz,
You must have some kind of demodulator - this is the part of your code that converts the received signal into LLRs. The histogram and averaging methods will never give identical results, but the results should be close to each other. You may like to compare the results you get from running my code with the results that I’ve provided at the top of this page. Can you provide a link to a .jpg of your results?
Take care, Rob.
July 12th, 2012 at 4:08 pm
Rob, thank you so much for your help.
I see that by increasing SNR in your code, the discrepancy between histogram and averaging charts as well as EXIT functions in each graph decreases.
This is the graph I get from running my code with averaging method:
http://www.sendspace.com/file/epsckt
And, this is the one with histogram method:
http://www.sendspace.com/file/69eon8
(please use the download link at the bottom of the page)
Elnaz
July 12th, 2012 at 4:10 pm
Hi Elnaz,
These EXIT functions are significantly different and so I suspect that there is a bug somewhere in your code. This could be in your turbo encoder, modulator, demodulator or turbo decoder. You should try to narrow it down…
Take care, Rob.
July 12th, 2012 at 5:12 pm
Dear Rob,
I just want to make sure I understand correctly. When you say “If your demodulator is working, then these plots will match each other, suggesting that the problem is in the turbo code” do you mean that e.g. the histogram graphs of different SNR values should match each other? Or that for each SNR value the histogram graph should match the averaging graph?
What I see with running your code is that by changing SNR value graphs significantly change as well but by increasing the SNR the similarity between histogram and averaging graphs grows too.
Am I correct?
Elnaz
July 13th, 2012 at 7:19 am
Dear Rob,
can you please extend your MATLAB code for EXIT demapper (natural QPSK mapping) to half rate regular LDPC EXIT chart (natural QPSK mapping). Please I really need it. Please help me in this regards. I will be thankful to you.
July 13th, 2012 at 5:18 pm
Hi Elnaz,
I mean that for each SNR value, the histogram graph should match the averaging graph.
Take care, Rob.
July 13th, 2012 at 5:22 pm
Hi Ideal,
I’m afraid that I don’t have any Matlab code that combines LDPC with QPSK. The QPSK demodulator has an input and an output pertaining to the LDPC-encoded bits - in main_exit.m from QPSKEXIT.zip, the input is called a_a and the output is called a_e. Likewise, the variable node decoder has an input and an output pertaining to the LDPC-encoded bits - in main_exit_vnd.m the input is called c_a and the output is called c_e. It’s just a matter of connecting these inputs to these outputs - you should also use an interleaver and deinterleaver though.
Take care, Rob.
July 14th, 2012 at 3:10 am
Hi Rob,
I have some results. At higher SNR values e.g. SNR=2, my graphs almost match. Those graphs that you saw earlier were from SNR=-4.
Then, according to your advice above, this means that there is something wrong with the decoder. However, I replaced my encoder and decoder which are your codes exactly only with different transition matrix, and yet I get almost same graphs as of when using my encoder and decoder!
At this point, I guess, I’m left without a clue as to what is wrong with my turbo decoding.
What would you suggest?
Thanks,
Elnaz
July 14th, 2012 at 12:48 pm
dear Rob,
thanks for the suggestion. I will try it to combine the ldpc with qpsk. but how can I add interleaver and deinterleaver? is there any matlab command or function which perform interleaving deinterleaving. thanks for your support
July 16th, 2012 at 7:19 am
Dear Rob
I find the QPSK Exit script written by you but I cannot find the main_exit_vnd.m file in you resources. please can you update for me thanks
July 16th, 2012 at 6:49 pm
Hi Elnaz,
I’m afraid that you will need to track down the bug in your code and fix it. I would suggest starting from my code and gradually converting it into your code, checking that the averaging and histogram methods give the same result after each stage of the conversion. In particular, you may like to double check that your trellis and component encoder correspond to each other. Something you may like to do is rewrite the encoder so that it is trellis-based, rather than shift-register based. This way, you can copy-and-paste the transitions from the decoder to the encoder, allowing you to be sure that they match…
Take care, Rob.
July 16th, 2012 at 6:52 pm
Hi Ideal,
An interleaver and a deinterleaver can be implemented in Matlab as follows…
interleaver = randperm(length(bits));
interleaved_bits = bits(interleaver);
deinterleaved_bits(interleaver) = interleaved_bits;
You can find main_exit_vnd.m at…
http://users.ecs.soton.ac.uk/rm/wp-content/LDPC.zip
Take care, Rob.
July 16th, 2012 at 7:48 pm
Hi Rob,
I’m not clear on one point as to why you think there is a bug in the encoder/decoder part of my code? since when I replaced my encoder and decoder with your “component_decoder” and “component_encoder” codes I still get almost the same graphs at SNR=-4.
Having this am I correct in assuming that it is the modulation/demodulation part rather than the coding/decoding part?
Elnaz
July 17th, 2012 at 8:37 am
Hi Elnaz,
You are right - this suggests that the problem is in the modulator or demodulator.
Take care, Rob.
July 17th, 2012 at 12:26 pm
Hi Rob,
Thanks; Therefore, my takeaway point would be, correct me if I’m wrong, that depending on the modulation/demodulation used Histogram and averaging graphs can look different on low SNRs (depending on modulation).
Since what I do before turbo decoding (I call it modulation) is not a standard modulation technique but the model of some sort of channel, I guess it is burdening the decoding part and cause this diiference at low SNR.
Here is my BER plot. Could you please take a look at it and let me know what you think?
http://www.sendspace.com/file/d7brlw
Thanks,
Elnaz
July 17th, 2012 at 6:23 pm
Hi Elnaz,
Your BER is heading towards zero as the SNR increases, so your system is as least partially working. In my experience, the averaging and histogram methods should agree with each other, even at low SNRs. It sounds to me like your method for converting the received signal into LLRs is not quite right. You may like to use my display_llr_histograms.m Matlab code for checking this. You can see lots of discussion about this in the comments of…
http://users.ecs.soton.ac.uk/rm/resources/matlabexit/
Take care, Rob.
July 18th, 2012 at 6:28 am
Hi Rob,
I am trying to combine the QPSK (natural) demaper exit with with LDPC exit but still cannot get the clue. Please help me in this regards a bit more. Thanks for all your help
July 18th, 2012 at 9:13 am
Dear Rob,
Now I understand the concept of QPSK demapper plus ldpc exit chart but still confusion in implementing it it matlab. I have initiallt two questions.
1. I am not understanding the iterations between demapper + variable node decoder + check node decoder.
2. Do I need to implement the soft demodulation bitwise (symbol-wise) or at first I calculate the channel LLR a vector and then apply demapping. and again don’t know how to implement iterations?
Please help me in this regards. Thanks a lot
July 18th, 2012 at 5:14 pm
Hi Ideal,
You can see an example of implementing iterations in main_ber.m, which you can download from this page. When you perform demodulation, you will need a for loop that demodulates each symbol individually.
Take care, Rob.
July 19th, 2012 at 2:44 am
Dear Rob,
I did some kind of flow for the transmitter and reciever for first iteration.
Transmitter:
codebits –> interleaver –> group/mapping(QPSK) –> adding noise and txt.
Receiver:
Demapper input: a_a(zeros at first iteration) + received symbols
ouputs: a_e and then perform the deinterleaving
variable node decoder inputs: a_e + a_a(from check node)
outputs: c_e and then perform interleaving
now 2nd iteration start at the input of demapper with a_a = c_e(interleaved) with channel signal.
But I am thinking to introduce the soft values at the input of the demapper as was done by ten Brinks or the above procedure is OK? Thanks
and I think I need to have for loop for the above process.. Thanks
July 19th, 2012 at 4:36 pm
Hi Ideal,
Your flow looks good to me. I’m not sure what you mean when you say “I am thinking to introduce the soft values at the input of the demapper”. It seems to me that you are already doing this when you say “2nd iteration start at the input of demapper with a_a = c_e(interleaved) “.
As you say, you can use a for loop to perform the iteration. One thing to think about is what decoder activation order you will use within each iteration of your loop. Examples of potential activation orders include:
1) demodulator - variable node - check node
2) demodulator - variable node - check node - variable node
3) demodulator - variable node - demodulator - variable node - check node
4) demodulator - variable node - check node - variable node - check node
and so on…
I would recommend option 1, since it activates each of the decoders equally often…
Take care, Rob.
July 20th, 2012 at 6:46 pm
Hi Rob,
Are you familiar with Cramer-Rao bound analysis?
Thanks,
Elnaz
July 23rd, 2012 at 6:38 am
Dear Rob,
Thanks for all your fruitfull suggestions. I appreciated it. I am really interested in calculating the EXIT function of the variable node decoder(VND). I did implement it (QPSK gray and natural mapping) and hopefully the results are good.
I just did the iterations between the demapper and the variable node decoder as,
ist iteration on transmitter side:
demapper input (chanel plus a priori LLR, (L_ma) (all zeros)) = L_me (out)
VND input (a priori LLR from check node plus L_da(L_me(Interleaved))) = L_de (out)
2nd iteration
demapper input= same channel + L_ma, where L_ma = L_de(de-interleaved) and same out L_me
VND input (a priori LLR from check node plus L_da(L_me(Interleaved))) = L_de
and so on.
EXIT chart
I can calculate the VND EXIT curve using L_de and the demapper curve using L_me.
Now one question:
I want to implement puncturing in this setup. please can you give me feedback about the puncturing implementation please. Thanks again
July 23rd, 2012 at 4:59 pm
Hi Elnaz,
I’m afraid that I have never used the Cramer-Rao bound…
Take care, Rob.
July 23rd, 2012 at 5:08 pm
Hi Ideal,
The simplest way to implement puncturing is to do it at the same time as interleaving.
You can design the interleaver/puncturer using…
interleaver = randperm(number_of_ldpc_encoded_bits);
interleaver_puncturer = interleaver(1:number_of_bits_after_puncturing);
You can perform puncturing in the transmitter using…
interleaved_punctured_bits = ldpc_encoded_bits(interleaver_puncturer);
You can perform puncturing in the receiver using…
apriori_interleaved_punctured_llrs = extrinsic_ldpc_encoded_llrs(interleaver_puncturer);
You can perform depuncturing in the receiver using…
apriori_ldpc_encoded_llrs = zeros(1,number_of_ldpc_encoded_bits);
apriori_ldpc_encoded_llrs(interelaver_puncturer) = extrinsic_interleaved_punctured_llrs;
Take care, Rob.
July 24th, 2012 at 11:31 am
Dear Rob,
I am really thankful for your responses. One question about the soft demodulation.
In your script of soft_demodulate.m, you are converting from symbol to bit extrinsic llr (a_e). Now my question is that do we need to take into account the two variable node as well, i.e., one symbol consist of two bits (x1×0), so for each bit one variable node of ldpc code, it means for one symbol, we need to use two variable nodes connecting with check nodes. I am confused here?
2ndly, If I want to plot the EXIT curve between the check node and demapper, I think, do I need to calculate the demapper curve after performing the complete iterations and using the llr (a_e) value from last iteration( using in you average MI script) so that the demapper EXIT curve shows some slop. Can you tell me the right way please. Thanks again
July 24th, 2012 at 5:45 pm
Hi Ideal,
There is no need to modify the soft demodulator to take the variable nodes into account. Taking care of the variable nodes is the job of the variable node decoder.
There are a few ways of plotting the EXIT charts of your 3-stage concatenated scheme. The way that I would recommend is to treat the LDPC decoder as a black box that the demodulator exchanges LLRs with. You can then draw a single EXIT function for this black box. You can achieve this by using the following process for each value of Ia in the set 0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0:
1) Generate some random bits
2) Encode these bits using the LDPC encoder, in order to generate the encoded bits
3) Use these encoded bits as the input to generate_random_llrs, together with the current value of Ia
4) Give the resultant LLRs to the LDPC decoder and perform 10 decoding iterations between the variable node and check node decoders, then obtain the extrinsic LLRs out of the variable nodes
5) Measure the mutual information Ie of these extrinsic LLRs
You can then plot Ie vs Ie.
Take care, Rob.
July 26th, 2012 at 1:42 am
Dear Rob,
Thanks for the simple example.
One thing please. I want to use your modulate.m and soft_demodulate.m script for the 4-PAM constellation, where 4-PAM symbol ={0,1/3,2/3,1} and the symbol energy I calculate E_s = 1/M(0+1/3+2/3+1)d=1/2d (optical),
can I use your scripts by just changing the constellation labels and bit labels and what should I need to multiply with symbols for normalization, i.e.,
constellation_points = [0; 1/3; 2/3; 1]/(which factor rather than sqrt(2)).
Thanks
July 26th, 2012 at 5:19 pm
Hi Ideal,
Yes - it is as simple as changing the constellation points and bit labels. You should normalise by using…
constellation_points = constellation_points/sqrt(mean(abs(constellation_points).^2))
Take care, Rob.
July 27th, 2012 at 1:50 am
Hi Rob,
Your code is very helpful. I need to generate LLR values of the coded bits at the decoder. I tried your code and it didn\’t work.
% Calculate encoded extrinsic transition log-confidences. This is similar to
% Equation 2.18 in Liang Li’s nine month report or Equation 4 in the BCJR paper.
deltas2=zeros(size(transitions,1),length(apriori_encoded_llrs));
for bit_index = 1:length(apriori_encoded_llrs)
for transition_index = 1:size(transitions,1)
deltas2(transition_index, bit_index) = alphas(transitions(transition_index,1),bit_index) + uncoded_gammas(transition_index, bit_index) + betas(transitions(transition_index,2),bit_index);
end
end
% Calculate the encoded extrinsic LLRs. This is similar to Equation 2.19 in
% Liang Li’s nine month report.
extrinsic_encoded_llrs = zeros(1,length(apriori_encoded_llrs));
for bit_index = 1:length(apriori_encoded_llrs)
prob0=-inf;
prob1=-inf;
for transition_index = 1:size(transitions,1)
if transitions(transition_index,4)==0
prob0 = jac(prob0, deltas2(transition_index,bit_index));
else
prob1 = jac(prob1, deltas2(transition_index,bit_index));
end
end
extrinsic_encoded_llrs(bit_index) = prob0-prob1;
end
Does it computes the LLRs of (lc+ld+le+lf) ? I need just lc.
Thanks,
Ushi
July 27th, 2012 at 1:59 am
Dear Rob,
Thanks for the response.
You also suggest me some iterative loop for QPSK but I am thinking, if I just generate the random LLRs and give them to the variable node input and then take out the extrinsic LLRs from the variable node and measure the MI of this extrinsic LLR. In this situation, I didn’t see any iterations between the demapper and the variable node. Thats why I was doing the following steps for 3,6 half rate code so dv=3.
Receiver side:
first iteration demapper:
1)Lm_e=perform soft demodulation with all zeros a priori LLRs (from VND).
2) Ld_e=input the Lm_e to the VND plus all zeros priori LLRs(incoming from CND (dv-1)).
2nd iteration demapper:
1)Lm_e=perform soft demodulation with a priori LLRs (ld_e from VND).
2) Ld_e=input the Lm_e to the VND plus a priori LLRs(from CND, this time not all zeros (random)).
this iterative process carry on until the iteration count end or MI from the demapper stops improving.
After that I got the final Ld_e from the above iterative process and I use that to measure the VND extrinsic MI for each point of the IAs and this is called VND EXIT chart?
July 27th, 2012 at 6:37 pm
Hello Ushi,
The code you have provided will calculate extrinsic LLRs that pertain to the output bits of the corresponding component encoder. Using the notation from my main_ber.m file, these LLRs would be c_e for the first component decoder and d_e for the second component decoder.
Take care, Rob.
July 27th, 2012 at 6:48 pm
Hi Ideal,
You are close getting it right, but your variable node has d_v+1 number of ports - one connected to the demapper and d_v connected to the check nodes. Therefore you should be generating d_v number of random apriori LLRs and you should be measuring the mutual information of d_v number of extrinsic LLRs. You currently have d_v-1 number of LLRs in both cases.
The LLR that is output on a particular one of the variable nodes ports should be equal to the sum of the LLRs that are input at the other d_v ports. So,
1) the LLR that is passed to the demapper by the variable node should be equal to the sum of the d_v number of randomly generated LLRs.
2) the LLR that the variable node outputs on a particular one of its check-node-connected ports should be equal to the LLR provided by the demapper and the d_v-1 randomly generated LLRs that are provided on the other check-node-connected ports.
Take care, Rob.
July 28th, 2012 at 1:59 am
Dear Rob,
Thanks, this is right now. One thing about the soft demodulation. I am wondering that you are using the same equation (2) in ten brinks paper “Iterative demapping for QPSK modulation” Electronic letters July 1998?
I think that equation (2) is calculating bitwise LLR value. If the soft_demodulate.m is the output of same equation then I am wondering where is the first additive term. If not then what is the possibility of using that equation. I mean how can we implement that equation?
Thanks for discussion
July 30th, 2012 at 5:45 am
Dear Rob,
One more thing in plotting the trajectory of QPSK. I am using the following code. It works well in Gray mapping but not for natural mapping.
IA = 0;
IE = 0;
results = [IA,IE];
snr_L = 10^(SNR*0.1);
sig_Ch = sqrt(8 * 0.5 * snr_L);
% Trajectory (dv=3,dc=6)
for iteration_index = 1 : 100
sig_A = (interp1(J,SIGMA,IA));
sig_D = sqrt(sig_Ch^2 + 3*(sig_A.^2)); % one channel + 3 apriori
sig_E = sqrt((3 - 1).*sig_A.^2 + sig_D.^2); % VND extrinsic
IE = interp1(SIGMA,J,sig_E);
results = [results; [IA,IE]];
sig_C = interp1(J,SIGMA,(1-IE));
sig_C(sig_C>20) = 20; sig_C(sig_C20) = 20; sig_EC(sig_EC<=0) = 0;
IA = 1 - interp1(SIGMA,J,sig_EC);
results = [results; [IA,IE]];
end
plot(results(:,1), results(:,2),’g-.’); hold on
July 31st, 2012 at 6:00 pm
Hi Ideal,
Equation (2) in the paper you mention calculates the aposteriori LLR. By contrast, my soft demodulator calculates the extrinsic LLR, which is equal to the aposteriori LLR minus the apriori LLR. The first additive term is the apriori LLR. Therefore, my soft demodulator is equivalent to the right hand term in equation (2) of that paper.
It looks to me that your trajectory code is using ten Brink’s equations for the VND and CND EXIT functions. These hold for Gray mapping, but not for natural mapping, as you have found. When using natural mapping, you will need to use your simulation to get the VND EXIT function, rather than an equation.
Take care, Rob.
August 1st, 2012 at 4:40 am
Dear Rob,
It means in order to calculate the EXIT chart, we can use the extrinsic LLR rather than the aposteriori LLR. So we can just implement the equation (2) without the a prioi LLR and fed that extrinsic LLR into the LDPC decoder after interleaving?
Or we still need the aposteriori LLR to fed to the LDPC decoder? If that is the case then how can we introduced the a priori LLR in simulations? Thanks
August 1st, 2012 at 4:42 pm
Hi Ideal,
Yes, it is the extrinsic LLRs that should be iteratively exchanged with the variable node decoder. The aposteriori LLRs should only be used once the iterative decoding process has finished, in order to provide the final output. If the aposteriori LLRs are used during the iterative decoding process, then this will become a positive-feedback loop and the performance will be severely degraded.
Take care, Rob.
August 3rd, 2012 at 2:45 am
Dear Rob
I see in your script of calculating the demappe exit curve that you wrote
N0 = 1/(SNR_linear), then in the channel model sigma^2 = (No/2) and in the demodulation No according to the pdf of n. But I am not sure where you define the code rate R?
I think, sigma^2 = (2RmRc(EbNo_Linear))^-1, where for QPSK Rm=2 and Rc is the code rate. Because ten Brinks define the demapper EXIT curve Ie = J(sqrt(8Rc(EbNo_Linear))) for gray mapping in his paper “Design of low density parity check code for modulation and detection” equation (19) and plot the results in Fig. 6 for 1X1. It gives different result from your demapper script. Can you please comment? Thanks
August 4th, 2012 at 2:54 pm
Hi Ideal,
There is a difference between SNR and Eb/N0 - I am using SNR. The relationship between them is Eb/N0 (dB) = SNR (dB) - 10*log10(R*log2(M)), where R is the channel coding rate and M is the number of constellation points, e.g. M=4 in QPSK.
Take care, Rob.
August 6th, 2012 at 5:38 am
Dear Rob,
Thanks, it means that I can replace N0 with (R*log2(M)*EbNo_Lin)^-1 in your script of soft_demodulate in order to represents it in Eb/No. I think I am right now? Thanks again
August 6th, 2012 at 5:45 pm
Hi Ideal,
That looks correct to me.
Take care, Rob.
August 15th, 2012 at 4:10 am
Dear Rob,
Thanks but I have again a basic question please.
Actually for symbol mapping like QPSK/4PAM, we write the symbol energy as Es = 1/M(\sum|x|^2) in electrical domain and Es = 1/M(\sum|x|) in optical domain and Eb/No = (1/R*log2(M)) * Es/No, My confusion is that, when we define No, we simply write No=(R*log2(M)*Eb/No)^-1, where is the value of Es?
I was thinking in BPSK Es = 1, that’s why we don’t need to write it but for the case of QPSK/4-PAM, Es will not be 1 or I am wrong. Please can you rectify me. Thanks a lot
August 15th, 2012 at 5:12 pm
Hi Ideal,
Typically, people design all of their modulation schemes to have Es = 1. For example, in QPSK we multiply the constellation points by 1/sqrt(2), in order to make Es = 1. In 16QAM, we multiply the constellation points by 1/sqrt(10) to achieve this.
Take care, Rob.
August 16th, 2012 at 6:17 am
Dear Rob,
Thanks for answering. I am thinking to do some derivations for natural mapping, in order to calculate the EXIT chart. can you give me any advice in matter please? Thanks
August 16th, 2012 at 4:34 pm
Hi Ideal,
The first step is to calculate the probabilities P(y|x) for each constellation point x, where y is the received signal. After that, you can convert this into a probability for each bit value.
Take care, Rob.
August 17th, 2012 at 12:26 am
Dear Rob,
It looks quite OK. I think if I calculation the conditional probabilities for each constellation point assuming the AWGN channel for each constellation point. the second point I didn’t understand, like why should I need to convert the conditional probabilities into a probabilities for each bit? can you elaborate it bit more please. or any reference, from where I can have a start-up please please. Thanks
August 17th, 2012 at 5:32 pm
Hi Ideal,
You need to calculate probabilities for each bit because you want to obtain LLRs for each bit. References on BICM-ID will give you the mathematics behind all this.
Take care, Rob.
August 20th, 2012 at 3:15 pm
Hi Rob,
I was wondering if there is a known-to-be best performing equalizer the same way as we have LDPC code as the best performing coding strategy?
Basically, I am convolving a matrix of coded bits with a 2-D pulse shape which introduces ISI in both dimensions and I need to use equalizer to mitigate the problem.
Thanks,
Elnaz
August 20th, 2012 at 5:14 pm
Hi Elnaz,
Non-linear equalisers are the best performing. A good example is a turbo equaliser, which uses the BCJR algorithm.
Take care, Rob.
August 20th, 2012 at 9:15 pm
Hi Rob,
I will be doing joint iterative decoding, timing recovery, and equalization. I have done the first two but I don\’t know how should I combine my Turbo-decoder (BCJRs) to include equalizer (Turbo with BCJR) as well. Are there any resources or codes that I can benefit from? What do you suggest?
Thanks,
Elnaz
August 21st, 2012 at 2:20 pm
Dear Rob,
I am still not clear how to start calculating the analytical derivation of EXIT chart for natural QPSK mapping. Please elaborate it bit more. Thanks
August 21st, 2012 at 5:59 pm
Hi Elnaz,
My advice would be to keep the turbo equaliser separate from the other system components. You can then connect the various system components through interleavers and then perform iterative decoding among them all.
Take care, Rob.
August 21st, 2012 at 6:04 pm
Hi Ideal,
I didn’t realise that you are trying to derive an analytic expression for the EXIT function of the QPSK demapper. I’m afraid that I’ve only ever considered the simplest of cases for analytic EXIT functions (namely LDPC variable and check node decoders) - so I don’t know where to get started with this. It looks like this reference may be of use to you though…
http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4114253&tag=1
Take care, Rob.
August 21st, 2012 at 6:49 pm
Hi Rob,
I know that there are 3 blocks of timing recovery, Turbo decoder, and Equalizer which should iterate together. What I don’t know is that how should I exchange the output llrs between the decoder and the equalizer. Should I exchange extrinsic llrs or full llrs? and how?
Also, I am under the impression that “Turbo Equalizer with BCJR” is the same “Turbo decoder” I already have with the only difference that in the equalizer we output the llrs for the coded bits instead of the llrs for the message bits. Am I correct? Or the equalizer is totally something else?
Thanks a lot,
Elnaz
August 22nd, 2012 at 4:34 am
Dear Rob,
This old paper,
read.pudn.com/downloads144/doc/comm/626085/a2.pdf, seems to be doing exactly that. But instead of keeping the turbo equalizer separate, it combines turbo decoder with an equalizer to build a turbo equaliser.
More important, my concern is that since I don’t have access to the channel transfer function (my channel is when I “2-D convolve” a matrix of coded bits with a 2-D pulse shape), how can I use BCJR in my equalizer since I don’t have the trellis of the channel?
Thanks,
Elnaz
August 22nd, 2012 at 2:42 pm
Hello Dear Rob,
I want to use a general turbo code (which uses M-QAM or M-PSK) in a system, when space-time coding is used too. Channel is Rayleigh fading.
Could you please guide me?
Best Regards
ALI
August 22nd, 2012 at 5:51 pm
Hi Ali,
The turbo code that is provided here is compatible with M-ary space-time coding, although I’m afraid that I don’t have any Matlab code for the latter. You just need to build a space-time encoder that accepts bits from the turbo encoder and a space-time decoder that provides corresponding LLRs for the turbo decoder.
Take care, Rob.
August 22nd, 2012 at 6:00 pm
Hi Elnaz,
“how should I exchange the output llrs between the decoder and the equalizer. Should I exchange extrinsic llrs or full llrs? and how?” You should exchange extrinsic LLRs that pertain to the turbo-encoded bit sequence. In order to obtain this, you will need to modify component_decoder.m so that it also outputs extrinsic_encoded_llrs. Also, you will need to obtain extrinsic LLRs for the systematic bits by adding the uncoded_extrinsic_llrs that are provided by both component decoders. Finally, the last three bits in the extrinsic_uncoded_llrs can be extracted to give you extrinsic LLRs for the termination bits.
“Or the equalizer is totally something else?” The equalizer should exploit knowledge about the dispersion in the channel to recover LLRs about the transmitted bits. It can represent the dispersion using a trellis and apply the BCJR algorithm to get the LLRs.
If your receiver does not have knowledge of the channel dispersion then you will not be able to use an equaliser unless you use channel estimation to obtain this knowledge.
Take care, Rob.
August 22nd, 2012 at 8:50 pm
Hi Rob,
I see. Therefore, turbo equaliser is/cen be an equaliser + turbo decoder combined with turbo principle.
My question, therefore, boils down to “how can I transfer my channel into a trellis?”. All my channel is doing is convolving a matrix of coded bits with a 2-D pulse shape. I have the pulse shape but I do not know how to convert it to a transfer function in a suitable form that can be translated into a trellis.
Is there a way to do this?
Thanks,
Elnaz
August 23rd, 2012 at 12:21 pm
Hello dear Rob
Thanks for your help. where is that file? Does it have adjustable rate?
Best Regards
Ali
August 23rd, 2012 at 5:58 pm
Hi Elnaz,
You need a trellis having M^(T-1) states, where M is the number of constellation points you are using and T is the number of consecutive received symbols that are affected by a particular transmitted symbol. There should be M number of transitions that emerge from each state and M number of transitions that merge into each state.
The trellis I have described is shown in the following paper, where L = T-1…
http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1267050
Take care, Rob.
August 23rd, 2012 at 6:01 pm
Hi Ali,
You can download my turbo code from…
http://users.ecs.soton.ac.uk/rm/wp-content/matlabturbo.zip
This doesn’t have an adjustable rate, but you could modify it to do this by using puncturing or by repeating some of the encoded bit sequences.
Take care, Rob.
August 24th, 2012 at 7:58 pm
Hi Prof. Rob,
Thanks for your attention
Best Regards
Ali
August 26th, 2012 at 4:05 pm
Hi Prof.Rob
Could you please give me a Code that achieves to capacity (or closely near to capacity, like LDPC code). If MATLAB code is available, please give me.
Best Regads
Ali
August 26th, 2012 at 7:10 pm
Excuse me, which algorithm does your code follow? log-MAP, max log-MAP or another one?
Thanks
Ali
August 26th, 2012 at 8:43 pm
hi Dear Prof.
Excuse me, I have another question:
Let, we have two space-time code. First code achieves better capacity than the second one., but in a state, the second code has better BER. Now, we want to use a strong outer channel coder. Can we say that because first code has better capacity in all SNRs, we can use a good channel encoder to proof this claim?
August 28th, 2012 at 6:35 pm
Hi Rob,
why does a turbo encoder work worse than a LDPC code?
thanks
August 29th, 2012 at 8:32 am
Hi Ali,
I have some Matlab code for drawing the EXIT chart of an LDPC at…
http://users.ecs.soton.ac.uk/rm/resources/matlabturbo/#comment-10965
…you can modify this to build an LDPC decoder.
My Matlab code for the turbo decoder can use the Log-MAP, Appox-Log-MAP (which uses a lookup table to approximate the Jacobian logarithm) and the Max-Log-MAP. You can choose between these by modifying jac.m
The space-time code that offers the greatest channel capacity is the one that can give the best performance when combined with a convolutional code.
Take care, Rob.
August 29th, 2012 at 8:39 am
Hi Jack,
The relative advantages and disadvantages of LDPC and turbo codes is a matter of debate. My opinion is that because LDPC codes have longer interleavers than the equivalent turbo codes, the LDPC designer has more design choices to make, yielding a greater opportunity to find a good design.
Take care, Rob.
August 29th, 2012 at 6:37 pm
Hi Rob,
I have a question about how to use your BCJR decoder as an equalizer.
Say, the channel model (pulse shape) is a vector of 5 coefficients which makes a 4-tap delay line. How should I change your transition matrix to include this coefficients in the tap-delay line?
Thanks,
Elnaz
August 30th, 2012 at 6:03 pm
Hi Elnaz,
The transition matrix should have M^T number of rows. Each row should one column for the from state, one column for the to state and T number of columns that list every possible combination of T constellation points. The from state depends on the first T-1 constellation points and the to state depends on the last T-1 constellation points. You also need log2(M) number of columns to list the bits that are mapped to the last of the T constellation points.
Take care, Rob.
August 30th, 2012 at 6:28 pm
Hi Rob,
I know how your code and transition matrix works for the decoding job. What I was asking is how should I change your code to be an equalizer instead of a receiver? Could you please explain what parts of the code should be changed? For example, what are the uncoded_apriori_llrs here since all the bits are channeled and there are no systematic bits to use for this purpose.
Thanks,
Elnaz
August 31st, 2012 at 6:21 pm
Hi Elnaz,
The calculation of alphas, betas, deltas and extrinsic llrs should remain very similar to how they are now. The big change will be in the calculation of the gammas.
Take care, Rob.
September 3rd, 2012 at 2:57 am
Hello dear Rob,
Thank your for your constructive expressions. Excuse me, How can I change your matlab code for turbo code to M-ary modulations. you have presented these sentences:
% BPSK demodulator
% These labels match those used in Figure 2.11 of Liang Li\’s nine month report.
a_c = (abs(a_rx+1).^2-abs(a_rx-1).^2)/N0;
c_c = (abs(c_rx+1).^2-abs(c_rx-1).^2)/N0;
d_c = (abs(d_rx+1).^2-abs(d_rx-1).^2)/N0;
e_c = (abs(e_rx+1).^2-abs(e_rx-1).^2)/N0;
f_c = (abs(f_rx+1).^2-abs(f_rx-1).^2)/N0;
What will be new commands?
Thanks alot
September 3rd, 2012 at 6:12 pm
Hello Ali,
Anything higher than M=2 becomes a lot more complicated that these equations. You can download a soft demodulator for higher order modulation from…
http://users.ecs.soton.ac.uk/rm/wp-content/QPSKEXIT.zip
Take care, Rob.
September 4th, 2012 at 5:11 am
Dear Rob,
Please correct me if I’m wrong. Transferring BCJR decoder into an equalizer, I have my transition matrix for 16 states and the encoded bits here are numbers in 0:16 range (there is no feedback path; so generating the transfer matrix is a lot easier and can be done automatically). Therefore, I think, I no longer can do log(P0/P1) normalization when calculating encoded-gammas the way you’ve done.
However, I think, I calculate |channeloutput-encoded bit|^2 /(2*variance) for each row of my transition matrix separately; and then, add those values that correspond to uncoded bit=0 in prob0 and the rest in prob1 and that should do it. Everything else stays the same. What do you think? Correct?
Elnaz
September 4th, 2012 at 6:41 pm
Hi Elnaz,
Nothing that you have said here seems wrong to me - my suggestion would be to build it and then test it by comparing the EXIT curves obtained using the averaging and histogram methods - these should match each other.
Take care, Rob.
September 5th, 2012 at 3:21 am
Hi Rob,
I have changed:
\"if transitions(transition_index, 4)==0
encoded_gammas(transition_index, bit_index) = apriori_encoded_llrs(bit_index);
end\"
to
\"encoded_gammas(transition_index, bit_index) = -(encoded_input(bit_index)-transitions(transition_index,4))^2/(2*sigma2); \"
which is calculated for all the lines in the transitions matrix.
Also, I\’ve changed the transitions matrix according to the figure 4 in the paper you referred earlier. These are the only changes I\’ve made to convert the BCJR decoder into an equalizer. Also, I\’ve set apriori_uncoded-llrs to all zeros.
Now, I expect that if I pass an stream of bits(-,+1s) through the channel (convolve the bits with E2PR2 h = [1 4 6 4 1]) and then add the noise and then through the equalizer, the extrinsic_uncoded_llrs from the equalizer should have the same sign as of the original message bits. I think even before going to BER or EXIT plots, this simple sign check can show if the equalizer is working.
I don\’t get the result I expect i.e. the signs don\’t match.
Do you know what might be wrong? How can I debug? Have I done correctly so far you think?
Thank you,
Elnaz
September 5th, 2012 at 5:21 pm
Hi Elnaz,
Your encoded_gammas needs to be a function of not only the constellation point in transitions(transition_index,4), but also the constellation points that relate to the previous symbol periods.
The reason why the signs may not be matching could be because you have strong dispersion. To do your test properly, you may like to remove dispersion - this will effectively transform your equaliser into the soft demodulator of…
http://users.ecs.soton.ac.uk/rm/wp-content/QPSKEXIT.zip
Take care, Rob.
September 5th, 2012 at 5:36 pm
Dear Rob,
I did not understand the first part. Did you mean that my calculation for encoded_gammas is wrong? I am putting the power in the exponential term as my encoded_gammas for each line in the transitions matrix separately, i.e. I calculate |channeloutput-encoded bit|^2 /(2*variance) for each row of my transition matrix separately (channeloutput is simply the ISI-channeld output).
In your comment before the last one, you mentioned it was correct. Have I missed something?
Thanks,
Elnaz
September 5th, 2012 at 6:15 pm
Hi Elnaz,
You only have one encoded gamma, not one per column in the transitions matrix. This is because you only have one received symbol.
I’m afraid that I must refer you back to for a clearer explanation…
http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1267050
I don’t think that I will be able to explain it any better than the example given in Figure 5 of that paper.
Take care, Rob.
September 6th, 2012 at 4:21 pm
Hi Rob,
I am testing my BCJR equaliser with channel response E2PR2 (h=[1 4 6 4 1])).
The part that I do not understand is that I should discard the last 4 extrinsic_uncoded_llrs to get the signs match. If it is because of the zero-padding that the convolution in MATLAB does then I had to discard the first and last 2 numbers rather than the last 4.
I am trying this with both termination on and off cases and it does not make a difference (I am testing with no encoder just randn BPSK bits). Also, there is no additive noise here.
Also, even after discarding the last four, sometimes the first bits are incorrect. All the rest match.
Do you know why this is so?
Thank you,
Elnaz
September 7th, 2012 at 12:37 pm
Hi Elnaz,
I’m afraid that I don’t know why this is - I’ve never actually programmed a turbo equaliser myself. Sorry for not being much help.
Take care, Rob.
September 7th, 2012 at 6:57 pm
Hi Rob,
It’s alright. You’ve been a great help so far.
Do you know of a reference for BER plot for a BCJR equalizer for EEPR2 response? I’ve searched a lot and I was not able to find any.
Thanks,
Elnaz
September 8th, 2012 at 5:42 pm
Sorry Elnaz, I’m afraid that I don’t know where to find that BEE plot. Take care, Rob.
September 9th, 2012 at 10:29 pm
Dear Rob
Please if possible , can you provide me the code of the turbo code with DS-CDMA or WCDMA and the documentation of this work, because it is very important.
Thank you,
Belal
September 10th, 2012 at 7:08 pm
Hello Balal,
I’m afraid that I don’t have any code for CDMA.
Take care, Rob.
September 11th, 2012 at 9:10 am
Dear Dr Robe,
Thank you for all help,
And if possible, can i insert spreading code “walsh code” in to your code of UMTS with turbo code, becuase this work is a part of the requirements of my master thesis.
Thank you for all help,
Kind Regard,
Belal
September 11th, 2012 at 5:14 pm
Hi Belal,
It would certainly be possible to combine my turbo code with a Walsh spreading code. However, I’m afraid that I don’t have any Matlab code to do this.
Take care, Rob.
September 14th, 2012 at 12:51 am
Dear Rob,
I have one simple question about the QPSK demapper EXIT chart using the natural mappings (b_1,b_0) –> X(QPSK Symbol). Example
1. I have input bits for the QPSK mapper (b_1,b_0) mapp–> QPSK Symbol (using natural mapping).
2. add noise and transmit.
3. At the receiver I try to find the extrinsic llr of bit b_1 –>e_llr(b_1) which needs a priori information from bit b_0 –>a_llr(b_0).
But I see your QPSKEXIT chart script, you are generating an input random vector (a) and converting it into QPSK symbols and at the demapper you are demapping, i.e., calculating the extrinsic llr –> e_llr(a) with the a priori llr –>a_llr(a). I am confuse here that where can I relate the e_llr(b_1) and a_llr(b_0),
and on the other way around, you know that If I want to calculate the extrinsic llr of bit b_0, e_llr(b_0), I need a priori llr of bit b_1, a_llr(b_1), what will be the difference in this regards in your script in both way?.
can you please highlight this concept? Thanks and kind regards
September 14th, 2012 at 4:55 pm
Hi Ideal,
In soft_demodulate.m we have…
function extrinsic_llrs = soft_demodulate(apriori_llrs, rx, channel, N0)
Here apriori_llrs is a vector containing [a_llr(b_0), a_llr(b_1)] and extrinsic_llrs is a vector containing [e_llr(b_0), e_llr(b_1)].
Take care, Rob.
September 15th, 2012 at 11:32 am
Dear Rob,
Thats right and now if I want to calculate the extrinsic llr e_llr(b_1) with the a priori llr a_llr(b_0), then how can I separate these vectors conditioning on each other. Like I want the output extrinsic llr e_llr(b_1) given the a priori llr a_llr(b_0) or the extrinsic llr e_llr(b_0) given the a priori llr a_llr(b_1).
Thanks for your responses
September 15th, 2012 at 5:25 pm
Hi Ideal,
If you look closely at soft_demodulate.m, you will see that e_llr(b_1) is calculated based only on a_llr(b_0). Here, a_llr(b_1) does not affect this calculation. This is achieved using the following if statement…
if bit_index2 ~= bit_index
So, my Matlab code is already doing what you are asking for.
Take care, Rob.
September 17th, 2012 at 1:00 am
Dear Rob,
Thanks for that, now I got it. It means, if I want to see the extrinsic llr e_llr(b_0), I can extract the odd terms in extrinsic_llrs vector containing [e_llr(b_0), e_llr(b_1)] and even terms for e_llr(b_1). Am I right now?
and as you said the e_llr(b_1) is calculated based only on the a_llr(b_0) then in this way, I can see the histogram of e_llr(b_1) based on a_llr(b_0) and can calculate the extrinsic mutual information of the extrinsic llr.
I am really thankful for all your help. You are so great
September 17th, 2012 at 5:23 pm
Hi Ideal,
It looks like you’ve got it figured out now.
Take care, Rob.
September 19th, 2012 at 2:13 am
Dear Rob,
Thanks for the responses.
Now I want to compare the histogram of e_llr(b_1) given a_llr(b_0) with the normal pdf and want to see if these two fits.
I think, in order to perform these comparison, at first I can calculate the histogram of e_llr(b_1) using the hist function of matlab and then hold the plot and then using the normpdf(vector, mean(e_llr(b_1),0.5*mean(e_llr(b_1))) matlab function I can see if these two plots fit.
Am I right or give me any suggestion please?
Thanks
September 19th, 2012 at 8:05 am
Hi Ideal,
That sounds right to me.
Take care, Rob.
September 22nd, 2012 at 9:42 am
Hi Dear Prof.
Before two weeks ago, I wanted QPSK modulation for turbo coded system. and you proposed me this file:
http://users.ecs.soton.ac.uk/rm/wp-content/QPSKEXIT.zip
How can I see the BER performance of this system?
Thanks for your help
Ali
September 23rd, 2012 at 3:39 pm
Hi Ali,
For this QPSK scheme to fulfil its full potential, you should concatenate it with some kind of soft-in soft-out channel code. You could use the outer code from…
http://users.ecs.soton.ac.uk/rm/resources/matlabexit/
This will give you bit-interleaved coded modulation with iterative decoding (BICM-ID).
Take care, Rob.
September 24th, 2012 at 2:22 pm
Hi Rob,
I have some questions:
1. If I have the extrinsic LLRs of the decoder E=(E1,E2,…En), how to estimate the pdf p(E=Ei|x=-1) and p(E=Ei|x=+1) using a histogram of E?
2. How to determine the range and the bin width of the histogram?
3. After p(E=Ei|x=-1) and p(E=Ei|x=+1) have been calculated, how to get the mutual information Ie=I(E;X)?
Thank you very much!
Yue
September 24th, 2012 at 2:29 pm
Hi, I have some question.
If I select the different bins, the Ie has large different value?
If I select the smaller bin, Ie could be more accurate?
Thank you!
September 24th, 2012 at 6:08 pm
Hello Yue,
1. You can do this using my Matlab code at…
http://users.ecs.soton.ac.uk/rm/resources/matlabexit/#comment-246
2. My code will do this automatically for you. There is some discussion about how this is done at…
http://users.ecs.soton.ac.uk/rm/resources/matlabexit/#comment-158
3. You can do this using measure_mutual_information_histogram.m, which you can download from above.
Take care, Rob.
September 24th, 2012 at 6:10 pm
Hi Tui,
This requires a careful optimisation - if the bin width is too small or too big, the MI will be measured inaccurately. I do it by determining the bin width automatically, as discussed in my comment to Yue, immediately above.
Take care, Rob.
September 25th, 2012 at 9:50 am
Dear Rob,
I want to plot the histogram of the extrinsic llr and see either it is normal distributed or not?
Usually using your display_llr_histogram() function, I can see the distribution of the two bits in llr but when i want to plot histogram over this plot, it doesn’t plot well.
please give me advise what should I need to do? In the previous comments I told you about it but I see that doesn’t give me the right results.
Actually I want to plot like you script is giving the results of two bits in one llr but with histogram comparison. Thanks
September 26th, 2012 at 8:39 am
Dear Rob,
Thank you for your reply.
I still have some questions:
1. I notice that you use the numerical integration to calculate I_E which is an integration. Why don’t you multiply bin_width in accumulation for calculating I_E in program “measure_mutual_information_histogram.m”?
2. If the encoded sequence is an all zero sequence instead of round(rand(1,n)), how to use your program “measure_mutual_information_histogram.m” to calculate Ie giving E={E1,E2,…,En}?
Thank you again!
September 26th, 2012 at 12:45 pm
Hi Ideal,
I’m not sure why you want to plot a histogram over the plot that display_llr_histogram() gives you - this plot is already the histogram of the LLRs.
I think that you want to plot the histograms for the two QPSK modulated bits separately. In which case, you need to arrange the LLRs and the bits into two groups and then give each of these to display_llr_histogram() separately.
Take care, Rob.
September 26th, 2012 at 12:59 pm
Hi Tui,
1. I can’t remember the specific reason for this but I remember there an important one. Perhaps this is because the bin width gets cancelled out by something else. In any case, the histogram method doesn’t work if we multiply by the bin width…
2. measure_mutual_information_histogram.m assumes that the values of the bits are equiprobable, i.e. that the bit sequence contains an equal number of 0s and 1s. So it doesn’t work when all of the bits are 0s. Here are some versions that work when the bit sequence does not contain an equal number of 0s and 1s…
http://users.ecs.soton.ac.uk/rm/wp-content/measure_MI_averaging_generalised.m
http://users.ecs.soton.ac.uk/rm/wp-content/measure_MI_histogram_generalised.m
Note however, that when all the bits are zero, the MI will come out as zero. You should also note that these methods assume that the LLRs know that the bit values are not equiprobable. i.e. in the absence of any information, the LLRs should have values of ln(p0/p1), rather than values of 0.
Take care, Rob.
September 27th, 2012 at 7:19 am
Dear Rob,
I have a question to your comment 2 to Tui. Sometimes, we can sent all zero codewords, in that case the process of encoding can be cut down. According to the symmetric pdf of E, we have that p(E=Ei|x=+1)=p(E=-Ei|x=-1). Therefore, we can estimate p(E=Ei|x=+1) using histogram, and also get p(E=-Ei|x=-1) as well. Is it okay?
Besides, I want to ask you another questions:
1. If the length of code is too small, such as 11, I think that it is not accurate to calculate p(E=Ei|x=+1) and p(E=-Ei|x=-1) using histogram. How to deal with it?
Thank you very much!
September 28th, 2012 at 9:01 am
Hi Yue,
This will only work if a number of conditions are met:
- your code needs to be linear and symmetrical
- your modulation scheme needs to be symmetrical
In general, I would not recommend this approach. I would always generate random messages, rather than all-zero messages.
If the length of each frame short, you can simulate a number of frames and then concatenate them together. Then you can use the histogram method on the long concatenated sequence.
Take care, Rob.
September 30th, 2012 at 1:16 am
Dear Rob,
I want to plot the histogram after puncturing symbols and want to see that whether after puncturing, the pdf of LLR is still normal. That’s why I was trying to see the histogram plot and display_llr_histogram plot over each other and by this I can confirm that the llr is Normal after puncturing or not.
This was the objective behind. But I am not sure I was doing right or wrong? Thanks
September 30th, 2012 at 10:17 pm
Hi Rob,
Could you please explain how best I can list the error patterns in a certain order, given the generator matrix, for syndrome decoding. Basically, I want to do what \"syndtable\" in MATLAB does but in a simpler way.
Thanks,
Elnaz
October 1st, 2012 at 8:47 am
Hi Ideal,
I see - you want to compare before and after histograms. I think that this sounds reasonable. My expectation is that if you are using random puncturing, the two histograms will be very similar.
Take care, Rob.
October 1st, 2012 at 8:51 am
Hi Elnaz,
If the generator matrix is not too big, the simple way is to build a list of every possible input into the encoder and then use this to obtain every possible output. Then, for each possible pairing of outputs, you can do an XOR to get the error pattern.
Take care, Rob.
October 3rd, 2012 at 3:00 am
Hi Rob,
Do you know how can I generate a Hamming code parity check matrix (with parameter m) in systematic format? Basically, I’m trying to do what “hammgen” of MATLAB does in a simpler way.
Thanks,
Elnaz
October 3rd, 2012 at 6:29 am
Hi dear prof.
I want to use LTE intelreaver, how will be the MATLAB code after considering this matter, to find BER?
Best regards
October 3rd, 2012 at 5:37 pm
Hi Elnaz,
I’m afraid that I don’t know of any simpler way than hammgen…
Take care, Rob.
October 3rd, 2012 at 5:39 pm
Hi Ali,
All you need to do is replace the mentions of get_UMTS_interleaver in main_ber with get_LTE_interleaver. You can download get_LTE_interleaver.m from…
http://users.ecs.soton.ac.uk/rm/wp-content/get_LTE_interleaver.m
Take care, Rob.
October 3rd, 2012 at 7:18 pm
Hi dear Prof.
What is the range of LTE interleaver?
Are you sure that the other parts of MATLAB code are the same as the MATLAB code for UMTS interleaver?
Best Regards
Ali
October 5th, 2012 at 4:50 pm
Hi Ali,
The shortest LTE interleaver is 40 bits, the longest is 6144 bits. The other parts of the LTE turbo code are the same as in UMTS. To see this for yourself, you can take a look in the documentation for the standards in the links at the top of this page.
Take care, Rob.
October 8th, 2012 at 12:14 am
Dear Rob,
How can I use your modulate and soft_demodulate scripts for higher order mapping schemes like 16-QAM or 64-QAM?
and how can they be implemented in simulations, like for 16-QAM, each symbol has four bit, then how can I replace it in your script of calculating demapper exit curve, where you write modulate(a(start:stop)) and a_e(start:stop) = soft_demodulate(.) for QPSK (2 bits per symbol);
I appreciate your help. Thanks
October 8th, 2012 at 6:06 pm
Hi Ideal,
All you need to do is modify the constellation_points and bit_labels - every thing else is automatic.
In main_exit, you just need to change bit_count.
Take care, Rob.
October 11th, 2012 at 1:26 pm
Hi dear Prof.
Excuse me, I don’t know about channel Coding very much. Could you please help me to puncture your turbo code to rate 1/2.
I thank you for your lats helps.
Bests Ali
October 11th, 2012 at 6:07 pm
Hi Ali,
You can see how to implement puncturing at…
http://users.ecs.soton.ac.uk/rm/resources/matlabturbo/#comment-8929
Take care, Rob.
October 16th, 2012 at 4:00 am
Dear Rob,
Thanks for all your support. I again has one tricky question please.
Lets I got one demapper curve (y) for a particular SNR, i.e.
y = ax^3 + bx^2 + cx + d, SNR = 1dB, x=0:0.1:1 (a-priori MI), (a,b,c,d=constants).
It means I have three different variables, y,x,SNR. I need to store a large number demapper equations for a range of SNRs and want to utilize them somewhere else. How can I use interpulation in matlab function which can cover three parameters in which x and y are vector but SNR is number? I though I can use interp2 but not sure how to use it?
Please do you have any idea about this problem? Thanks a lot
October 16th, 2012 at 6:10 pm
Hi Ideal,
I’m afraid that I’ve never done that sort of thing before, so your guess is as good as mine…
Take care, Rob.
October 17th, 2012 at 5:53 pm
Hi Rob,
I have a question about the PLL with Mueller-Muller estimate. Do you know how it works?
Thanks,
Elnaz
October 18th, 2012 at 2:20 pm
Hi Elnaz,
I’m afraid that I have never worked on PLLs, so I can’t help you with this.
Take care, Rob.
October 22nd, 2012 at 5:25 am
Dear Rob,
Thanks for the fruitful discussion. I have very basic question please.
If I am considering the real valued 4-PAM constellation, i.e.,
constellation_points = [-1; -1/3; 1/3; 1]/sqrt(5/9); and the received signal is represented by y = Px + n; with P is the received signal power and n is the AWGN noise with zero means and unit variance noise. In that case, can I replace N0 = 1 in your soft_demodulate.m? Thanks
October 22nd, 2012 at 5:48 pm
Hi Ideal,
You should use…
constellation_points = sqrt(P)*[-1; -1/3; 1/3; 1]/sqrt(5/9);
y= x+n;
N0 = 1;
Take care, Rob.
October 23rd, 2012 at 10:56 pm
Dear Rob,
I don’t understand why did you write sqrt(P);
The symbol energy of the system considering optical to electrical domain is
Es = P^2(5/9), for the constellation_points = [-1; -1/3; 1/3; 1];, where P is the optical transmit power and considering the zero mean unit variance noise, The SNR=P^2(5/9);
October 23rd, 2012 at 11:12 pm
Hi Rob,
Could you please explain a stopping criteria for turbo decoder or actually for turbo equalizer when applied in a iterative receiver scheme?
Thanks,
Elnaz
October 24th, 2012 at 5:00 pm
Hi Rob,
Regarding my previous question, I’m using sum of LLRs for the stopping criteria. What do you think?
My new question:
I am simulating Turbo equalizer on a E2PR2 channel. I need to somehow convert the LLR outputs of the Turbo equalizer back to the channeld format. Right now I am taking the hard estimate (sign) of those LLRs and modulate them back again with the channel response which is not the optimum thing to do. How can I do this properly?
Thank you,
Elnaz
October 24th, 2012 at 5:56 pm
Hi Elnaz,
I normally stop the iterative decoding process when the mutual information of the extrinsic LLRs stops improving.
It would be possible to modify the turbo equaliser so that it not only outputs extrinsic information about the bits, but also extrinsic information about the modulated symbols. This would be more optimal than modulating a hard decision for the bits…
Take care, Rob.
October 24th, 2012 at 5:58 pm
Hi Ideal,
This square root is needed because amplitude is related to the square root of power. Here, P is relating to power, while x, y, n and sqrt(P) are all related to amplitude…
Take care, Rob.
October 24th, 2012 at 6:10 pm
Hi Rob,
Could you explain a bit more how should I modify the turbo equalizer to output information on the modulated bits?
Thanks,
Elnaz
October 24th, 2012 at 6:12 pm
Please add this to my previous comment:
Also, how should I change the LLRs into soft channel outputs?
October 25th, 2012 at 4:34 pm
Hi Elnaz,
This is easiest for BPSK, in which case the LLRs for the BPSK modulated symbols are equal to the LLRs for the bits. In the case of a higher order modulation scheme, you just need to output the symbol probabilities that are calculated during the turbo equalisation algorithm.
Take care, Rob.
October 26th, 2012 at 4:08 pm
Hi Rob,
What I meant to ask is that should I follow the same procedure that I used to get the encoded outputs from component_decoder in my BCJR equalizer as well i.e. to add the lines:
deltas2(transition_index, bit_index) = alphas(transitions(transition_index,1),bit_index) + uncoded_gammas(transition_index, bit_index) + betas(transitions(transition_index,2),bit_index);
And,
prob1 = jac(prob1, deltas2(transition_index,bit_index)); (I don’t think I can do this since the outputs are no longer 0s and 1s but some integer values)
Then, if correct, this will give me the extrinsic llrs for modulated bits?
Then, I can’t change these llrs to normal channel output format by simply dividing them by 2/sigma2 i.e. the way to go in case of decoder but not equalizer. My channel is E2PR2 response applied on BPSK symbols. Then I’m wondering how can I extract channel output from llr-outputs of my turbo equalizer.
Thanks,
Elnaz
October 26th, 2012 at 6:18 pm
Hi Elnaz,
I think that you are on the right track here. For each of the M symbol values (e.g. M = 4 in QPSK), you can use an equation that is similar to the prob1 equation that you have provided. The exponential of these prob values will give you values that are proportional to symbol values. You just need to normalise these values in order to get the correct probabilities.
Take care, Rob.
October 26th, 2012 at 6:35 pm
Hi Rob,
Do you mean summing over delta2s corresponding to each row in the transition matrix and then for example prob4 = exp(sum of delta2s corresponding to output=4)/sum of delta2s for all output values?
Thanks,
Elnaz
October 28th, 2012 at 9:59 am
Hi Elnaz,
That sounds right to me.
Take care, Rob.
October 29th, 2012 at 2:28 am
Hi Rob,
One additional question:
My simulations run verry slowly and according to MATLAB profiler, “conv” and “jac” commands are the major consumers. What can I do to make “jac”run faster?
Thanks,
Elnaz
October 29th, 2012 at 8:45 pm
Hi Elnaz,
If you look into jac.m, you’ll see that it has three possible modes of operation. I would suggest using the second mode, which approximates the exp and log functions using a look up table.
Take care, Rob.
November 16th, 2012 at 6:09 pm
Hi Rob,
As you know I am using your BCJR component_decoder as an equalizer after making some changes of course. However, I can’t use the code for long stream of bits e.g. 32 kbit in MATLAB. I get out of memory problem and it takes a long long time to equalize say 15 kbits. Is there anyway I can improve the code or the speed of its execution?
Thanks,
Elnaz
November 19th, 2012 at 10:47 am
Hi Elnaz,
The sliding window technique can be used to reduce the memory requirement of a BCJR decoder. You can read about this technique in…
http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4726112&tag=1
Take care, Rob.
November 20th, 2012 at 7:39 pm
Hi Rob,
Thanks for all your answers, they are really useful and helped me understand much more how you obtained the results.
To my understanding, your results are for a AWGN channel.
Assume we add Rayleigh fading as defined below :
n = 1/sqrt(2)*[randn(1,N) + j*randn(1,N)]; % AWGN o db variance
h = 1/sqrt(2)*[randn(1,N) + j*randn(1,N)]; % Rayleigh channel
% Channel and noise Noise addition
y = h.*s + 10^(-Eb_N0_dB(ii)/20)*n;
Is there anything we have to change in the BPSK demodulator ?
What about the BCJR decoder ?
Thanks again for your help, very useful.
Regards
November 21st, 2012 at 4:58 pm
Hello Sam,
Only the demodulator needs to be updated. The updated code should look like…
a_c = (abs(a_rx+h).^2-abs(a_rx-h).^2)/N0;
Take care, Rob.
November 23rd, 2012 at 5:48 am
Hi Rob,
I’ve tried that, and it’s not giving me very accurate results for BER vs ebno. I also changed frame length to 4000 instead of 40 and result remained the same.
Actually my curve is going up and down, it’s not a stable one. I wonder why.
As a reference, i get to a BER value of 0.01 for 10 dB.
I was wondering if i had to modify something else in the decoder, like for the calculation of gammas, alphas or betas ?
Thanks for your help,
Regards
Sam
November 23rd, 2012 at 4:32 pm
Hi Sam,
You don’t need to modify anything else. My suggestion is for you to test each component of your scheme separately. This will help you to determine if the problem is in your encoder, modulator, channel, demodulator, decoder, etc.
You can do this by drawing the EXIT curves for each component using both the averaging and the histogram methods. These two methods will give similar curves if everything is working okay. If not, then this will identify which component has the problem.
Take care, Rob.
November 27th, 2012 at 5:07 pm
Hi Rob,
In order to test it, i have used the same code as yours, i added the rayleigh and modified the part for the demodulator BPSK as you mentionned :
a_c = (abs(a_rx+h).^2-abs(a_rx-h).^2)/N0;
instead of
a_c = (abs(a_rx+1).^2-abs(a_rx-1).^2)/N0;
But that doesn\’t work neither.
I am still trying to figure out what isn\’t working,
What do you think ? Weird part is without Rayleigh, all works just fine.
Thanks again for your help,
November 27th, 2012 at 5:27 pm
Hello Sam,
You may like to try using the modulator and demodulator that can be downloaded from…
http://users.ecs.soton.ac.uk/rm/resources/matlabexit/#comment-1545
This code demonstrates how to use the demodulator in the case of a Rayleigh fading channel. You can convert the modulator and demodulator into BPSK by replacing the corresponding lines of code with…
constellation_points = [+1;-1};
bit_labels = [0;1];
Take care, Rob.
November 28th, 2012 at 12:28 am
Dear Rob,
one quick question about the variable node extrinsic mutual information (VNEMI) please.
1. we can calculate the VNEMI using the ten Brinks formula, i.e.,
I_ev = J(sqrt[(dv - 1)*sigma_A^2 + sigma_Ch^2]),
where
sigma_Ch^2 = 4/sigma_n^2 and
sigma_n^2 = 1/(2*log2(M)*(Es/No)) (noise variance, Es=symbol energy)
2. My question is that if we implement QPSK and 4-PAM, what will be the difference in sigma_n^2?. We know that the MI of both QPSK and 4-PAM are different. Please comment
November 29th, 2012 at 4:14 pm
Dear Rob,
I am facing the same problem with Sam and trying to demodulate BPSK after Rayleigh channel, I dont understand why it cant just be demodulated as normal e.g (rx+1)/2, where rx are the received values after transmitting through the Rayleigh and AWGN channel. I hope to gain better understanding in this. Just like you Sam, all works well in AWGN. I hope Rob could shed some light on this matter. Thanks
I think this is a great site!
November 29th, 2012 at 7:10 pm
Hello Ideal,
In this case, you need to use…
I_ev = J(sqrt[(dv - 1)*sigma_A^2 + [J^-1(I_demod)]^2]),
where I_demod is the mutual information of the LLRs provided by the demodulator.
Take care, Rob.
November 29th, 2012 at 7:13 pm
Hello Mcrave,
The simple answer is that because the channel is different, the received signal will be different and so the demodulator needs to be different. A more involved answer is to say that the demodulator needs to use h, in order to undo the phase rotation and fading that is imposed by the Rayleigh channel.
Take care, Rob.
November 29th, 2012 at 11:37 pm
Dear Rob,
Yes I do need to use, I_ev = J(sqrt[(dv - 1)*sigma_A^2 + [J^-1(I_demod)]^2]),
but still I_demod is the function of sig_Ch^2 and the a-priori values. So I am asking what will be the formula for sig_Ch^ in the case of gray maping QPSK and real 4-PAM modulation schemes.
I think in case of QPSK, we can use sig_Ch^2 = 4/sigma_n^2 = 8 * SNR but for the case of 4-PAM, what is this relationship? thanks
November 30th, 2012 at 4:13 am
Thanks Rob!
That clears up my concepts a little
so in other words if the transmitted BPSK signal is corrupted by a Rayleigh and AWGN channel, at the receiver, the first step to perform is equalization with rx_hat=rx./h then turbo decoding right? Also would just want to confirm that in practice, the channel coefficient estimated could be erroneous as it is an estimation right?
Thanks so much Rob
December 3rd, 2012 at 9:43 am
Hi Ideal,
As far as I am aware, there is no closed form equation that links I_demod to sig_Ch^2 in the general case. I think that the only thing you can do is simulate the demodulator and measure I_demod.
Take care, Rob.
December 3rd, 2012 at 9:47 am
Hi Mcrave,
That’s right. In practical scenarios, h would be obtained at the receiver using channel estimation. However, this can never be perfect and so a corrupted version of h must be used in practice. This will degrade the performance. In theoretical work however, it is typically acceptable to assume that perfect channel estimation is possible.
Take care, Rob.
December 3rd, 2012 at 2:45 pm
Thanks Rob. That’s very helpful. Another quick question, what is the difference of equalizing a_rx/h then taking the real part of the result as compared to what you have done for the demodulator with:
a_c = (abs(a_rx+h).^2-abs(a_rx-h).^2)/N0;
I would be very much thankful if you could point me in the right direction on this.
December 3rd, 2012 at 5:39 pm
Hi Mcrave,
If you divide a_rx by h, then you are not only amplifying the signal part of a_rx, but also amplifying the noise part. So, you would need to use a different value for N0 (although I can’t remember what value you should use)…
Take care, Rob.
December 3rd, 2012 at 5:58 pm
Hi Rob,
I see.. I think I have a better idea now. What a_rx/h does is zero forcing. I think that the codes provided might be implementing MMSE which takes into account of the noise power. I hope I’m not wrong, but I think the N0 is the noise power. N0=2/sigma^2. Thanks for the guide. Much appreciated!
Thanks again for your help.
December 3rd, 2012 at 6:04 pm
Hi Mcrave,
You are welcome.
Take care, Rob.
December 7th, 2012 at 2:18 am
Hello Rob,
I’ve just finished my number crunching for a Rayleigh channel with perfect channel knowledge with parameters below:
16 state CCSDS Turbo decoder
Modulation scheme:BPSK
No_of_FRAMES=2810;
INTERLEAVER_SIZE=1780;
ITER=5; %Number of iterations
RATE=1/2;
with different Rayleigh coefficient for each frame and have also changed the demodulator to a_c = (abs(a_rx+h).^2-abs(a_rx-h).^2)/N0; To take account of the Rayleigh channel. The results does not seem right however, it performs much poorly than expected. At 10db, a BER of 10^-2 was obtained. Which does not agree with other reported works. Any ideas what could have gone wrong? again, it all works fine in AWGN.
December 7th, 2012 at 6:09 pm
Hi Mcrave,
It sounds like you are using a quasi-static Rayleigh fading channel - i.e. you use only one value of h for every BPSK transmission in the frame. This will cause some frames to have low BERs, but other frames to have very high BERs, which complete eclipse the frames with low BERs - even when the SNR is 10 dB. My advice would be to try using an uncorrelated Rayleigh fading channel - i.e. you should use a different value of h for each BPSK transmission in a frame.
Take care, Rob.
December 8th, 2012 at 4:17 am
Dear Rob,
That makes sense, I did implement a sort of controlled seed of h values for each frame, but it is a different seed (i.e. different h for each frame). I think I’ll re-run the simulation again and let the randn function run loose without control. Fingers crossed! Thanks again Rob
You’ve been a great help 
December 12th, 2012 at 7:08 pm
Hi Rob,
Turns out that I didnt transfer correct noise sigma values for AWGN.
Just wanted to thank you for your advise and contributing to a working simulation 
Have a good day!
h*bpsk_tx + 10^(-ebnodb/20)*n;
Did the trick. Also, I made sure to have different h for every frame.
December 12th, 2012 at 7:09 pm
No problem! Take care, Rob.
December 30th, 2012 at 7:41 am
Hi Rob,
Hope all is going well. I just wanted to check whether my understanding is right on something. According to my simulation, for slow fading channels, i.e, the same h for each symbol in the frame yet different for each frame, the BER for turbo decoding is worse than a fast fading channel (different h for each symbols in the frame). After rechecking everything, I’m quite sure the problem earlier (with bad BER at high SNR~10dB) was because of this. Or that simply I would need to do a whole lot of simulation to obtain the correct BER result. From what I have observed, a channel with magnitude of ~0.1 has a very bad effect on BER and this seems to be the reason why I have been obtaining unexpectedly high BERs for high SNR. However, from the papers that I have been looking at, my result for fast fading channel seems to match the results reported in papers etc. When do I use fast fading or slow fading channels? Which do I stick to for publishing purposes? Thanks in advance for your advice
December 30th, 2012 at 8:47 am
Hi Rob,
Another thought, maybe its better if I plot the Frame error rate instead of BER? Thanks.
December 30th, 2012 at 1:40 pm
Hi Mcrave,
In general, I would recommend fast fading channels for physical layer work and slow fading channels for network layer work. This is because these are typically the simplest channels that sufficiently challenge these types of work. Of course, you may prefer to use more sophisticated channel models if this is your particular focus. For physical layer work, I typically use BER, because this more easily shows how the frame length affects the performance of an iterative decoding scheme. For network layer work, I recommend FER.
Take care, Rob.
December 31st, 2012 at 3:00 am
Hi Rob,
Thanks again for your advice
Much appreciated!
January 10th, 2013 at 6:25 pm
Hi Robert. Is it possible to generate bits (e.g 1000000 bits) instead frame (40 to 5114) in yours matlabturbo ?
January 10th, 2013 at 8:00 pm
Hi Luke,
It is possible to use any interleaver length, although the UMTS and LTE interleavers only support a limited range of lengths. You can instead use a random interleaver by setting random_interleaver=1. This will allow any interleaver length, such as 1000000 bits.
Take care, Rob.
January 11th, 2013 at 10:02 pm
Hi Rob,
I’m using BCJR-based turbo equalizer in my application and I have cross interference between different adjacent streams of BPSK message bits. In the receiver I somehow compensate for this interference but there is a weird behavior I see in my final BER curve. My BER curve is not always monolitically decreasing with increasing number of iterations and also with respect to increasing SNR values. I’m using sum reliability as my stopping criteria. However, if I use the error rate as my stopping criteria (i.e. if I cheat) to terminate iterations I do get monolitically decreasing curves. What do you think is the best terminating criteria to use here?
Let me also add that I think this is related to the high amount of interference I’m dealing with; and, for now, I should stick to my not-so-strong interference cancellation method and try to get smoother BER curves.
January 14th, 2013 at 3:43 pm
Hi Elnaz,
If successive iterations are not improving the BER, then it sounds to me like your equalizer (or some other part of your receiver) is not working properly. My advice would be to see if this is the case by comparing the EXIT curves of the equalizer using both the averaging and histogram methods of measuring the MI. If the curves obtained using these two methods are significantly different, then it suggests that your equalizer is not working properly…
Take care, Rob.
January 25th, 2013 at 8:25 pm
Hi Rob,
I have a question regarding the “eoor floor” problem. In order to get rid of this problem and have a nice, continuous falling down curve, do I increase the length of the stream or the number of averaging trials? Which one is more effective?
Thanks,
Elnaz
January 28th, 2013 at 1:29 pm
Hi Elnaz,
Increasing the number of trials in your simulation will not change the error floor - it will just give you a smoother and more accurate BER plot. Increasing the frame length will reduce the error floor. Another way of reducing the error floor is to use a better interleaver design - a good starting place is replacing a random interleaver with an S-random interleaver…
Take care, Rob.
January 31st, 2013 at 8:19 am
Hi Rob,
I’m trying to study the effects of SNR mismatch on turbo codes. In all literatures, the lowest BER will happen when there is perfect knowledge of
the SNR, i.e SNR offset=0. I find that my simulation fits the overall profile but the problem is that the lowest BER point is at offset -6dB. I think one of the reasons for this is probably the way I calculated the values for sigma etc. Here is a snippet of my program, hope you could point
me in the right direction
%% SNR is fixed
ebnodb=2;
ebno=10^(ebnodb/10);
N0=1/ebno;
sigma=sqrt(N0);
%% SNR Offsets
START=-8;
INTERVAL=1;
STOP=8;
for OFFSET=START:INTERVAL:STOP
ebno_off=10^((ebnodb+OFFSET)/10);
L_c_off=4*ebno_off;
%% Fast fading channel (COMPLEX)
h_u1=(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN))/sqrt(2); %ENC_OP_LEN is the length of the systematic bits u1, from RSC1.
h_u2=(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN))/sqrt(2);
h_v1=(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN))/sqrt(2);
h_v2=(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN))/sqrt(2);
%% AWGN noise (COMPLEX)
n_u1=sigma/sqrt(2)*(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN));
n_u2=sigma/sqrt(2)*(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN));
n_v1=sigma/sqrt(2)*(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN));
n_v2=sigma/sqrt(2)*(randn(1,ENC_OP_LEN)+1i*randn(1,ENC_OP_LEN));
%% Apply channel effect
rx_sig_u1= h_u1.*bpsk_tx_u1 + n_u1; % bpsk_tx are bpsk modulated bits
rx_sig_u2= h_u2.*bpsk_tx_u2 + n_u2;
rx_sig_v1= h_v1.*bpsk_tx_v1 + n_v1;
rx_sig_v2= h_v2.*bpsk_tx_v2 + n_v2;
%% Fix phase rotation and fading from channel with equalization
u1rx = real(rx_sig_u1.*conj(ho_u1));
u2rx = real(rx_sig_u2.*conj(ho_u2));
v1rx = real(rx_sig_v1.*conj(ho_v1));
v2rx = real(rx_sig_v2.*conj(ho_v2));
Then it goes on to perform turbo decoding. Would appreciate very much if you have any advice. Thanks!
January 31st, 2013 at 5:06 pm
Hi Mcrave,
I suspect that the mistake is in the lines that come after what you have provided. Your equaliser makes the signal appear to have been transmitted over an AWGN channel, rather than a Rayleigh fading channel. However, the equaliser will change the SNR - I guess that you are just using ebno_off.
My advice would be to use code of the sort…
tx = -2*(bits-0.5);
h = sqrt(1/2)*(randn(size(tx))+i*randn(size(tx)));
n = sqrt(N0/2)*(randn(size(tx))+i*randn(size(tx)));
rx = h.*tx + n;
llrs = (abs(rx+h).^2-abs(rx-h).^2)/estimated_N0;
Take care, Rob.
February 3rd, 2013 at 12:13 am
Hi Rob,
I am using the your BCJR algo as my equalizer in my application and the problem is that I see a “dip” in my BER curve in low SNR values.
I have tested one thing: when using the exact Jacobian logarithm in the code I have the “dip” in low SNRs but when I replace it with the approximation of mode=2 the problem is resolved and the curve is linear, slightly higher error at low SNRs but linearly decreasing with increasing SNR.
Do you know of any justification why this is happening? I mean why the exact definition gives me nonstandard behavior at low SNRs but the approximation doesn’t?
Thanks,
Elnaz
February 4th, 2013 at 6:45 pm
Hi Elnaz,
I’m afraid that I’ve never come across anything like that before. The only thing that springs to mind is perhaps one part of your system is assuming LLR=ln(P0/P1), while another part is assuming LLR=ln(P1/P0). I’m afraid that some of the code on my website assumes the former, while some of my code assumes the latter…
To track the bug down, you could try my usual recommendation, which is to plot the EXIT curve for each of your system components using both the averaging and histogram methods. If these don’t match for a particular one of your system components, then that’s where I think you’ll find the bug…
Take care, Rob.
February 9th, 2013 at 12:02 pm
Hi Rob,
Much thanks for your advice earlier. I did as advised but the results are not as expected. All of the mismatch results in 0 errors, only a severely underestimated channel with -8dB showed a slight error. There are a few things that I would need advice on.
1) the estimated_N0 in this line should be from ebno_off from my earlier post right? i.e estimated_N0= ebno_off=10^((ebnodb+OFFSET)/10); then this is applied to
llrs = (abs(rx+h).^2-abs(rx-h).^2)/estimated_N0;
2) I also read from a few sources now that Lc=4*a*Eb/No*R, where R=rate, and a=fading attenuation (a=1 when channel is gaussian I wonder if I need to change for Rayleigh, if so to what value?)
Again hoping you could point me in the right direction. It’s Chinese New Year from where I am, Happy Chinese New Year to all!
February 11th, 2013 at 2:47 am
Hi Rob,
I found something that might be useful in point 2). “In case of no fading a = 1. In the case of fading with no channel state information the fading value becomes the expected value a = E[a] of the underlying fading channel. When we have ideal channel state information available at the decoder then a takes the exact fading value.”
I am assuming a non CSI system. So a =E[h], for the codes suggested,
h = sqrt(1/2)*(randn(size(tx))+i*randn(size(tx)));
a=sqrt(1/2)?
Thanks.
February 11th, 2013 at 3:03 am
Think I might be wrong. Rayleigh channels are zero mean. Hmm…
February 11th, 2013 at 6:01 pm
Hi Mcrave,
What you are saying in point 1 looks correct to me.
With regard to point 2, the equation that I provided…
llrs = (abs(rx+h).^2-abs(rx-h).^2)/estimated_N0
…can be rearranged into the form…
llrs=4*real(rx*conj(h))/estimated_N0
I think that this is what you are talking about in point 2. Hopefully, my version of this equation shows you how to use it - although I don’t think that it will change your results…
Take care, Rob.
February 12th, 2013 at 3:06 am
Dear Rob,
Many thanks for your reply. I was able to get the expected results only when I apply
h = randn(size(t));
n = sqrt(N0)*(randn(size(tx));
rx = h.*tx + n;
rx = real(rx.*conj(h))/N0_off
Not sure why though. It is probably because of its modulation type BPSK and most papers focuses on the real part of the symbol sent? Most papers I see so far are all on BPSK.
February 12th, 2013 at 3:55 pm
Hi Mcrave,
In this code, you are using only real noise and fading - there is no imaginary component to your results. This makes things more difficult for your scheme because all of the noise and fading is in the same dimension as your BPSK modulation, namely the real part of the complex numbers. In my version, the noise and fading is shared between the real and imaginary parts, making things easier for the BPSK modulation, since it is immune to imaginary noise and fading. I guess that the results you are trying to reproduce are using only real-valued noise and fading, in which case I think you have got everything correct. However, it is normal practice to use complex-valued noise and fading…
Take care, Rob.
March 14th, 2013 at 5:47 am
Hi Rob,
I have a very basic question regarding BER performance. When no coding is applied we get the Uncoded BER. Then we use FEC code to bring down the BER.
Matlab 2012 has Turbo encoder- decoder demo (rate ~1/5.07). In this demo when I use BPSK (Bipolar) scheme, add AWGN noise and perform hard decision on received signal I got uncoded BER around 0.22; and Turbo codes were still giving coded BER ~10^-5. Is it possible?
Is there any coding technique (with code rate < 1/6) available to handle uncoded BER of 0.22??????
How much uncoded BER can be handled using present coding techniques?????
Thanks
March 15th, 2013 at 4:07 pm
Hello Mahesh,
Your results do not sound unreasonable to me. Turbo codes can have very steep BER curves, allowing low BERs to be achieved at relatively low SNRs. By contrast, in the absence of coding, the BER plot will have a very gradual gradient, requiring a high SNR in order to achieve a low BER. Turbo codes having long interleavers are near-perfect codes - they can achieve very low BERs when transmitting data at a rate (in bit/s) which approaches the capacity of the channel.
Take care, Rob.
March 15th, 2013 at 5:56 pm
Hi Rob,
When we say that the Ultimate Shannon limit is -1.59 dB, it means that no matter what code rate we use, we cannot perform reliable communication below -1.59 dB. Now for BPSK scheme, over AWGN channel, the uncoded BER at -1.59 dB is ~0.12. Means on an average 12% of bits (hard decoded) are in error. If turbo code can handle 0.22 uncoded BER isn’t that sound spooky. Because even in CCSDS 2012 report the turbo code with code rate 1/6 can go upto -0.2 dB Eb/No only. Which means it can handle only ~0.095 uncoded BER.
I am working on a channel coding technique. What should be my AWGN channel model ? I have used BPSK scheme with Bipolar (bit 0 is modulated using +1 and bit 1 is modulated with -1) scheme. Then I add AWGN using inbuilt matlab “awgn” function and set the SNR to 3 (SNR=3 so for BPSK scheme Eb/No=0 dB). I do hard decision at the receiving end and I got the uncoded BER ~0.0786, which matches with the theory.
I thought, channel coding is all about the amount of uncoded BER you can handle and everything in-between is just the process to achieve this goal.
Please correct me if there is anything wrong with my understanding:
i.e. Uncoded BER at -1 dB for BPSK is ~ 0.10
When we encode 100 bits with code rate 1/6 we get 600 coded bits.
When these 600 bits are transmitted at -1 dB Eb/No (or equivalent SNR) the received 600 bits will have ~0.10 * 600 bits in error
We need to correct these errors and bring down the BER to < 10^-5
I really appreciate your help.
Regards,
Mahesh Patel
March 18th, 2013 at 10:37 am
Hi Mahesh,
I suspect that you may be getting confused between SNR and Eb/N0. The relationship between these is…
Eb/N0 [in dB] = SNR [in dB] - 10*log10(eta),
where eta is the number of bits of information that are conveyed in each transmitted symbol. In uncoded BPSK, eta is 1, so Eb/N0 = SNR.
Take care, Rob.
April 18th, 2013 at 6:31 am
Hi Rob,
I have implemented my own Turbo Encoder and Interleaver. Now I want to use your decoder. But I am completely stuck in the receiver side programs. Your decoder takes LLR inputs and LLR generation program needs mutual information while the mutual information program needs LLR inputs.
So kindly will you please guide me from where will I start??
One more question, the component_decoder.m is one constituent decoder. So to implement the complete turbo decoder just feeding the output of one decoder to another will work? Or I not to make other modifications?
Waiting for your kind reply. Thanks in advance.
April 18th, 2013 at 4:33 pm
Hi Deep,
The functions for generating LLRs and measuring mutual information are only useful for drawing EXIT charts, to characterise the operation of the turbo decoder. These functions are not a part of the turbo decoder itself. Instead, your LLRs should be output by your demodulator. If you look in my main_ber.m script, you can see how I get a BPSK demodulator to provide LLRs to my turbo decoder. This script will also show you how to feed the LLRs produced by one component decoder into the other.
Take care, Rob.
April 19th, 2013 at 1:42 pm
Hi Rob,
Thanks a lot. I am unable to understand the code in two different places.
1. a_c = (abs(a_rx+1).^2-abs(a_rx-1).^2)/N0;
Is this the LLR output of the BPSK demodulator? If yes, then how is it in the LLR form? Is it in normal domain or logarithmic domain? Will you kindly please explain a bit.
2. errors = sum((a_p < 0) ~= a);
You have used this line to make a hard decision and count the number of errors. My question is how we will make the hard decision. This line will just calculate the number of errors.
Waiting for your kind reply. Thanks in advance.
April 19th, 2013 at 4:51 pm
Hi Deep,
1. This is the LLR output of the BPSK demodulator. It is in the LLR form - namely ln(P0/P1).
2. The hard decision is obtained by
decoded_bits = a_p < 0;
Take care, Rob.
April 19th, 2013 at 6:34 pm
Thanks Rob,
One more question, if I want to use my LLR as ln(P1/P0) then will the changing of the modulator and demodulator be sufficient? Or I have to make changes in the decoder? Will the alpha,beta,gamma,delta calculation be the same??
April 20th, 2013 at 6:47 am
Dear Rob,
Will you please explain a bit elaborately how
a_c = (abs(a_rx+1).^2-abs(a_rx-1).^2)/N0;
is in the ln(P0/P1) form?
April 22nd, 2013 at 6:10 pm
Hi Deep,
You will need to make changes in the decoder. You need to change
if transitions(transition_index, 3)==0
if transitions(transition_index, 4)==0
extrinsic_uncoded_llrs(bit_index) = prob0-prob1;
to
if transitions(transition_index, 3)==1
if transitions(transition_index, 4)==1
extrinsic_uncoded_llrs(bit_index) = prob1-prob0;
Note that the equation
a_c = (abs(a_rx+1).^2-abs(a_rx-1).^2)/N0;
gives an identical result to
a_c = 4*real(a_rx)/N0;
Perhaps you are more familiar with this version of the equation…
Take care, Rob.
May 6th, 2013 at 12:46 am
Hi Rob,
Thanks for the code. My question is why you use the rate 1/2 for SNR calculations in main_ber.m. To quote your comment, you said \"Only the information that is used by the upper decoder is transmitted in this simulation.\" but it seems like you use the parity bits from both encoders in main_ber.m as far as I understand. Then, the rate would be 1/3 as in standard UMTS turbo code. Could you please clarify this?
May 6th, 2013 at 1:24 am
Dear Rob,
I need to find the best LDPC code (parity check matrix) for my application. My block length =10 Kbits and rate=1/2 are fixed. Should I use PEG algorithm? Can you guide me to some codes for calculating the parity check matrix? I want to use the H with MATLAB’s LDPC encoder/decoder.
I tested the code you uploaded above but it takes a long time considering 10 Kbits block length. When I feed the resulted H to fec.ldpcenc I get the error : “the last (n-k) columns of the parity-check matrix must be invertible in gf(2).” Can you please explain why?
Thank you,
Elnaz
May 6th, 2013 at 10:09 pm
Hi Rob,
I would like to also know if this code is punctured to get rate 1/2, how should I modify the script to simulate the turbo code with original rate 1/3?
FYI, the reported UMTS BER vs Eb/No results on this website are significantly better than the ones I found here http://www.csee.wvu.edu/~mvalenti/documents/DOWLA-CH12.pdf and couple of other papers like this one: http://delivery.acm.org/10.1145/1340000/1330710/a105-iliev.pdf?ip=128.54.22.219&acc=ACTIVE%20SERVICE&key=C2716FEBFA981EF147DAD23CB8BFD2078D75470BC0F2371B&CFID=214026711&CFTOKEN=31039678&__acm__=1367864840_44a3b8b785b00ea67003e472fe6321ca.
Thanks in advance,
Selin
May 7th, 2013 at 8:02 pm
Hi Selin,
This is a rate-1/3 code, as you suspect. Please can you let me know where I say that it is a rate-1/2 code, so that I can fix it!
You can increase the rate to 1/2 by discarding some bits before passing them to the modulator - you can either discard a random selection of bits, or use a particular puncturing pattern. In the receiver, the discarded bits should be replaced, using zero-values LLRs.
Take care, Rob.
May 7th, 2013 at 8:06 pm
Hi Elnaz,
I’m afraid that I don’t know where to get a good algorithm for designing LDPC codes - in the past, I have always focused on the LDPC designs used in the standards. I think that Matlab’s decoder generates this error because it can’t work out the generator matrix from the parity check matrix you have provided it with - I’m afraid that I don’t know how you can fix this…
Take care, Rob.
May 7th, 2013 at 10:29 pm
Hi Rob,
Thanks for the reply. In main_ber.m file, you use this equation to calculate noise variance
% Convert from SNR (in dB) to noise power spectral density.
N0 = 1/(10^(SNR/10));
but σ^2 = 1/(2R *Eb/N0) and R = 1/3 in this case. So I think you are using a wrong variance value for your AWGN channel.
Do you happen to know about the MATLAB\’s new built-in Turbo Encoder and Decoder? I get worse performance in terms of BER with MATLAB\’s own implementation of this and wondered if you have any knowledge of this.
Thanks,
Selin
May 8th, 2013 at 5:10 pm
Hi Selin,
I agree that σ^2 = 1/(2R *Eb/N0) in the case of BPSK modulation. However, I don’t agree that this is inconsistent with N0 = 1/(10^(SNR/10)). I think that both equations are true.
I guess that you are referring to…
http://www.mathworks.co.uk/help/comm/ref/turboencoder.html
http://www.mathworks.co.uk/help/comm/ref/turbodecoder.html
I’m afraid that I don’t have the latest version of Matlab, so I haven’t been able to try these…
Take care, Rob.
May 8th, 2013 at 11:50 pm
Hi Rob,
I guess if you are defining SNR as Es/No, which is R*Eb/No, then yes, both equations would be correct. ( σ^2 = N0/2 = 1/ (2 SNR)) => N0=1/(SNR) = 1/(10^(SNRdB/10)) )
As for MATLAB\’s implementation, they provide the demo results for UMTS turbo code on this page http://www.mathworks.com/help/comm/examples/parallel-concatenated-convolutional-coding-turbo-codes.html#commpccc-13
Link for the BER curve:
http://www.mathworks.com/help/comm/examples/commpccc_03.png
These results are clearly worse than those provided here:
http://www.csee.wvu.edu/~mvalenti/documents/DOWLA-CH12.pdf
Do you have an opinion or sources you can refer me to help me figure which one is correct?
Thanks,
Selin
May 9th, 2013 at 12:51 pm
bsr,svp pouvez vous me donner les differrent code rate (pattern puncturing)des reseaux telephonie(LTE,UMTS,GPRS,CDMA2000,CDMA,HSDPA,HSUPA,HSP+…).j’ai besion de ces codes rates urgent et
merci d’avance
May 9th, 2013 at 5:53 pm
Hi Selin,
That’s right - I’m using SNR as Es/N0.
It looks like the Matlab turbo decoder provides sub-optimal performance for some reason. Both my results and the results in the book you have linked to are superior. If I were you, I would trust that book over Matlab…
Take care, Rob.
May 9th, 2013 at 5:55 pm
Hello Awatef,
I’m afraid that I don’t have any code for these different puncturing patterns - I normally focus on unpunctured codes…
Take care, Rob.
May 9th, 2013 at 5:55 pm
Hello Awatef,
I’m afraid that I don’t have any code for these different puncturing patterns - I normally focus on unpunctured codes…
Take care, Rob.
May 9th, 2013 at 6:50 pm
@ ROB
je vous remercie tout d’abord pour votre reponse.E nsuite
j’ai besion de connaitre les défferentes motifs de poiçonnage non les codes.c.a.d j’entrain de faire une recherche sur les les défferentes motifs de poiçonnage afin de les coder plus tard.
merci
May 9th, 2013 at 10:00 pm
Hi Rob,
Thanks for taking the time to look at those graphs. I decided not to use MATLAB’s turbo decoder since all the other turbo code performances seem to be matching except for MATLAB’s. Did you publish any of your BER results in a paper so I can refer to them as well? (I’m a PhD student at UCSD and using turbo codes as part of my research)
Awatef,
You might want to have a look at the CML library. They have puncturing methods there. http://www.iterativesolutions.com/download.htm
Also, for optimal puncturing patterns for certain constituent encoders (unfortunately not UMTS, but still gives you an idea), you can refer to this paper:
http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=848555&tag=1
Cheers,
Selin
May 10th, 2013 at 5:17 pm
un grand merci selin.
mais les fichiers sont (.mat)en exel.ils ne s\’ouvrent plus.
merci
May 10th, 2013 at 5:23 pm
Hi Selin,
We have published some LTE BER plots in…
http://eprints.soton.ac.uk/271820/
Thank you for offering some help to Awatef. Your French must be better than mine!
Take care, Rob.
May 10th, 2013 at 5:24 pm
Hi Awatef,
I’m afraid that my French is rather rusty and so I’m not sure what you are asking for. I wonder if you can restate the question in English?
Take care, Rob.
May 10th, 2013 at 7:41 pm
a big thank you Selin.
but the files that are in (CML Source Code) are files with extension (. Mat) (in Exel) so they s’ open more.
thank you
May 10th, 2013 at 7:42 pm
a big thank you rob
May 10th, 2013 at 8:06 pm
Hi Rob,
Thanks a lot. All the BER curves but Matlab’s are indeed matching.
Awatef,
If you download the cml.1.10.zip, you can see that under cml/mat folder, there are MATLAB files with .m extension. And depending on what you exactly want to simulate, you can just run the CmlSimulate.m with the already existing scenarios by only changing the puncturing pattern as you wish. .mat files are just outputs they simulated before and provided for convenience. Hope that was helpful.
Selin