The Consultative Committee for Space Data Systems (CCSDS) standard uses a Serially Concatenated convolutionally coded Pulse Position Modulation (SCPPM) that features a convolutional code with only four states but 7560 trellis stages for rate-1/2. Traditional decoders work sequentially 979-8-3503-5597-0/25/$31.00 ©2025 IEEE through these many stages, constraining throughput. This paper presents an alternative method that simultaneously processes all trellis stages allowing a dramatic increase in throughput and lowering the latency of decoding. Parallel trellis processing can be realized on both Graphics Processing Unit (GPU) and Field Programmable Gate Array (FPGA) implementations. GPU implementation can provide multiple benefits: easier development, scalability, and portability across GPU models. With the rapid advancement of technology, the same implementation remains functional on future GPU devices and compiler versions.
Convolutional codes are typically decoded via the Viterbi algorithm [1], [2] that minimizes the frame error rate (FER), or the algorithm introduced by Bahl, Cocke, Jelinek, and Raviv [3] (BCJR), that minimizes the symbol error rate. Low delay requirements have been a limiting factor for Viterbi decoders. Several works have achieved a level of parallelism to address this problem. Fettweis and Meyr [4] combined every P stages of an K-stage trellis using a smaller trellis per each of the M states, and achieved a linear speedup of P . Mohammadidoost and Hashemi [5] optimized the memory access for both the surviving path and the trace-back path on a GPU. Peng et al. [6] introduce a method to achieve a high throughput also using GPU. Their method is flexible enough for application on general block codes and is shown to attain a throughput of up to 1.8 Gb/s on convolutional codes. Zhihui et al. [7] and Hanif et al. [8] show that part of the computations can be done independently. In [7] a sequence is split into smaller segments and the processing of each segment is further split into an independent part and a dependent part. A matrix representation of the process is used in [8] that allows a GPU to exploit the independence of some operations. Seshadri and Sundberg [9] compare convolutional codes with Low Density Parity Check (LDPC) codes under delay constraints and show that convolutional codes outperform LDPC in the more restrictive regime. An FPGA version of a Viterbi decoder is also provided by Ben Asher et al. in [10], with parallelism that can be tuned via a parameter P . Concatenated convolutional codes like the SCPPM code are often decoder iteratively, see the parallel concatenated Turbo codes by Berrou et al. [11] and serially concatenated codes by Benedetto et al. [12]. In iterative decoding, two or more decoders alternate to produce a-posteriori probabilities that are used as input to the other decoders in the next
Information Source Append CRC Pulse Position Modulator Convolutional Encoder Accumulator Interleaver Channel PPM Symbol Log Likelihoods Soft BCJR Accumulator De-interleaver Soft BCJR CC ENCODER DECODER CRC Check Poisson Channel Interleave Figure 1: SCPPM Encoder model. The information generates a bit sequence to which a 32-bit error detecting CRC sequence and two zeros are appended, to obtain a K-bit, e.g. K = 7560 when the rate is 1/2. The K-bit sequence is encoded with an outer convolutional code and the output is interleaved with a quadratic interleaver. The interleaved sequence is encoded with a two-state recursive convolutional code or accumulator, and the output sequence is modulated with pulse-position modulation.
iteration. The process repeats until a stopping condition is met. The SCPPM code stops when no error is detected by an error detecting sequence or when the maximum number of iterations is exceeded. The serially concatenated (SC) pulseposition modulation (PPM) code consists of two serially concatenated convolutional codes with an interleaver between the two codes and is depicted in Fig. 1.
A convolutional code is a linear block code that encodes an information sequence into another coded sequence of symbols. The input sequence could be a sequence of bits, and the outputs could be bits or words comprised of two or more bits each. A convolutional code (CC) consists of a set of states {0, . . . , M -1} and two sets of linear functions
}, where f i maps an input symbol to one or more output symbols and h i defines a state transition that maps input symbol from source state i to a destination state j. The ratio of k input symbols to n output symbols per encoding operation defines the code rate k n . The model of the SCPPM encoder is depicted in the left side of Fig. 1. At the beginning, an error detecting 32-bit cyclic redundancy check (CRC) code sequence is appended to the information sequence. Then, the information sequence, with the 32-bit CRC sequence, is encoded with a four-state, nonrecursive, zero-terminated convolutional code, which will be called the "outer code." The output of the outer code is interleaved with a quadratic interleaver, and the interleaved sequence is encoded with a two state recursive convolutional Transmitter SCPPM Iterative Decoder SCPPM Encoder Receiver Free Space Optics The transmitter emits a pulse of light at each symbol's pulsed slot, and the receiver counts the number of photos received at each slot. The input to the decoder are the photon counts for every slot of every symbol.
code, which will be called the inner code or accumulator.
The output sequence of the inner code is split into symbols that are between two and eight bits long, and each symbol is modulated with a pulse position modulation.
Serially concatenated convolutional codes are traditionally decoded by message-passing decoders, with a separate BCJR decoder implemented for each component convolutional code. Each component BCJR decoder takes as input the output of the other component BCJR decoders and the channel symbol probabilities, according to the encoding graph.
Each decoder computes new a-posteriori probabilities that are passed along to the other component BCJR decoders.
The BJCR decoders traditionally compute the reliabilities by working sequentially through the trellis first in the forward direction and then in the backward direction [3]. This paper presents a method to break the BCJR algorithm [3] dependency on the values at each previous stage and instead simultaneously compute many stages to increase throughput and lower the latency of decoding.
The rest of the paper proceeds as follows: Sec. 2 describes the SCPPM code and a version of the BCJR algorithm by Bahl et al. [3] used by the standard SCPPM decoder. Sec. 3 introduces our method to transform the BCJR algorithm to allow higher parallelism and lower the decoding time. Sec. 4 describes the implementation of the algorithms using a GPU and CUDA. Sec. 5 shows simulation results in terms of runtime performance and frame error rate vs. average photon rates, and Sec. 6 concludes the paper.
Let the input of a code be a sequence with entries in a finite field U and outputs from a finite field F , which could be the as U . A K-symbol input sequence U = U 0 , U 1 , . . . , U K-1 is mapped to an N output sequence X where N = K n k . In this work, the value of k will be k = 1 and the output symbols are n-symbol words. Given an initial state S 0 , a convolutional code maps a K-symbol input sequence U to a unique output sequence X and a unique sequence of state S 1 , S 2 , . . . , S K , via the output and state transition functions. If the code is terminated, additional symbols are appended to the input symbol sequence U, to ensure that the last state S end , generally same as S 0 , is reached. In the next two sections, we will call the input to the inner code a K-bit sequence U and the output a N -bit sequence X, which is mapped to M -ary PPM symbols C 1 , C
2 , . . . with M ∈ {4, 8, 32, 64, 128, 256} each encoding m bits, where m ∈ {2, 3, . . . , 8}. The encoder's modulated output is a vector of symbols C 1 , . . . , C N/m , where symbol C i is a Time Symbol Times Symbol Symbol Symbol Pulsed Slot Non-pulsed Slots Slot Times Codeword Time Each symbol consists of 3-bits, that represent a position in 0, 1, . . . , 7. The 1 st modulated symbol is 011 = 3, the 3 rd modulated symbol is 010 = 2, the i-th modulated symbol is 100 = 4, and the last symbol is 101 = 5. A pulse is transmitted for each symbol, at positions 3, 2, 4, 5 for the 1 st , 3 rd , i th and last symbols respectively. No pulse is sent at any of the other slot positions on each symbol.
and C i has a single 1 entry. For convenience, each m-bits symbol constructed from the output of the inner code will be called an edge E, where E ∈ {0, 1, . . . , 2 m -1} and the input to the inner code will also be grouped into symbols
Since the inner code is rate 1, the value of k for the inner code is will be k = m.
The channel model is a Poisson channel defined by a background noise rate of K b average photons per non-pulsed PPM slot and K s + K b average photons per pulsed PPM slot.
The received signal for each transmitted symbol i is vector Y i,j and for each PPM slot j ∈ {0, . . . , 2 m -1} Y i,j is modeled by the probability mass functions (p.m.f.) Pr(Y i,j = y | C i,j = 0) and Pr(Y i,j = y | C i,j = 1) given by:
From now on, denote by p i,j (y) ≜ Pr(Y i,j = y | C i,j = 1) the channel's p.m.f. in Eq. ( 2).
The BCJR algorithm of Bahl et al. [3] is a standard component of the decoder of concatenated codes because it produces a-posteriori probabilities of its input and output sequences.
The inputs to the BCJR algorithm are the a-priori probabilities of the input and output sequences, where one of the two might be missing. The channel sequence symbol probabilities in (1) and ( 2) provide the output symbol probabilities for the inner-most decoder. In serially concatenated codes, the information sequence symbol a-posteriori probabilities are the output of only the outer-most decoder, thus the outermost decoder does not get updated a-priori probabilities of the information sequence from other decoders. The outputs of the BCJR algorithm are the vector of a-posteriori probabilities of the input sequence U 1 , U 2 , . . . , U K given by Pr(U i = u | Y = y) for each u ∈ U and the vector of a-posteriori probabilities of the output sequence X 1 , X 2 , . . . , X nK given by Pr(X i = x | Y = y) for each x ∈ F . Given code functions f i and the channel p.m.f., we can directly compute the probabilities Pr(Y i = y i | C i,j = 1) via:
For the rest of this section we explain how to compute the output probabilities Pr(
) and a-priori probabilities Pr(U i = u; I) produced by the outer decoder, where I stands for intrinsic input probabilities. Note that for the inner decoder the input sequence U is the output sequence of the outer decoder. In the next two sections we will refer to the input of the inner decoder by
and the output a sequence of edges
that would be produced by the inner are not used and thus are omitted. These probabilities are straightforward to obtain, see [12]. Next we express the channel probabilities
as functions of the input and state S i , that is, as
and each state s ∈ {1, 2, . . . , M } as follows:
We are interested in the conditional probabilities: Pr(U i = u | Y = y), which may be obtained from the joint probabilities Pr(U i = u, Y = y), i = 1, 2, . . . , N and the marginal Pr(Y = y) by:
) BCJR on Outer Decoder Convolutional De-interleave BCJR on Inner Code (Accumulator) Check CRC interleave Soft Extrinsic From Channel Soft Intrinsic Soft Intrinsic Soft Extrinsic
The inner decoder produces "extrinsic" log-likelihood ratios of its input, which is de-interleaved and used as input by the outer decoder. The outer decoder produces an estimate of the input bits and extrinsic log-likelihood ratios of its output sequence. If the estimate sequence is not a codeword of the CRC code, the extrinsic output of the outer code is interleaved and combined as intrinsic input with the channel probability logs and used as input in next iterations by the inner decoder. This process continues until a codeword of the CRC code is found or the max iterations are exceeded.
The bottom Pr(Y = y) is a constant given Y that can be ignored in many cases, or it can be computed via u∈U Pr(U i = u, Y = y). Then, is suffices to compute the probabilities Pr(U i = u, Y = y) as follows.
= s e:fs(u)=e
Since every next output is E i = f s (u) and every next state is S i+1 = h si (u), then the event {U i = u, S i = s} guarantees that Pr(Y = y, U i = u, S i = s, S i+1 ̸ = h s (u)) = 0, and:
The rest of this section focuses on computing the probabilities:
Let the beginning state of an edge e be b(e) and the end state be t(e). Also denote by Y j i = y i , Y i+1 , . . . , Y j , then:
1 ) as functions of the channel signal Y, the channel p.m.f and the a-priori probabilities
) may be computed by:
• Pr(
Since
. . , N/m are obtained from the channel via equations ( 1), ( 2) and (7).
Then:
= Pr(S i-1 = b(e), Y i-1
• Pr(
Computing the probabilities (25), ( 26) and (27)
Let the input symbol U i consists of m bits, that is:
, with b j ∈ {0, 1}. Denote each bit of U i by b i,j . We wish to compute the probabilities of each bit: by:
The probabilities
via Eq. ( 9), with:
The values of σ i,j (b), for b = 0, 1 and j = 1, 2, . . . , n, i = 1, 2, . . . , N/n may be obtained from α i-1 (b(e)), γ i (e), β i (t(e)) as follows:
We have considered a set of edges that necessarily includes repeated edge patterns. That is, if the number of inputs is 2 k and the number of states is 2 ν , then there are 2 ν 2 k edges to consider. If the length of the edges is n < ν +k then there are 2 ν+k-n edges with the same bit sequence. A different form that unambiguously accounts for all edges.
The values of σ i,j (b) may be computed using λ i,j (s ′ , s, b)
instead of γ i (e) as follows:
• Pr(
• Pr(
Let the input that produces edge e as output be u(e) ≜ f -1 b(e) (e). Then, for each edge e and each i = 1, 2, . . . , N/n the values of γ i (e) are computed from the channel probabilities Pr(Y i = y i | E i = e) and the a-priori probabilities Pr(U i = u(e)) as follows:
Note that two edges e and e ′ may have the same value but have different source state b(e) ̸ = b(e ′ ) or destination state t(e) ̸ = t(e ′ ). In such case, we will still treat them as different so that the source and destination states of each edge are unique and the event {S i-1 = b(e)} is unambiguous. However, the channel probabilities will be shared by these edges.
Let γ i (s ′ , u) ≜ γ i (f s ′ (u)), then the probabilities λ i,j (s, s ′ , b) may be computed from the values of γ i (e) or γ i (s ′ , u) as Figure 6: Binary tree representation of the process to compute the values of α i (s) for each i = 1, 2, . . . , N/m in a logarithmic number of steps, represented by the tree depth. First the values of δ i2 j+1 ,i2 j+1 +2 j (s ′ , s) are obtained for j = 1, 2, . . . , L, and each i = 1, 2, . . . , 2 L-j and each s, s ′ = 1, 2, . . . , M and then used to compute all α i (s).
follows:
The values of α i (s) and β i (s) are generally computed recursively, where α i (s), s = 1, . . . , M are used to obtain α i+1 (s) and β i+1 (s), s = 1, . . . , M are used to obtain β i (s), via the Forward and Backwards algorithm [3]. In this paper, we develop an algorithm to break the recursive dependencies, which is explained in the next section. The standard method to compute α i (s) and β i (s) for all i = 1, . . . , N and all s = 1, . . . , M follows:
1 , then:
The values of β i (s) are computed starting from β N/n (s) backwards as follows:
and the event {Y N i+1 = y N i+1 } is conditionally independent of {Y i = y i , S i-1 = s} given {S i = s ′ }, then:
Equation ( 61) concludes the review of the standard method to obtain the values of α i (s), β i (s) and γ i (e) needed to compute each σ i,j (b). However, with the recursive equations ( 54) and (61) the values of α i (s) need to be computed before the values of α i+1 (s) and the values of β i (s) need to be computed before the values of β i-1 (s), for each i = 1, 2, . . . , N and each s = 1, 2, . . . , M . In the next section, we break this dependency to enable computing the values of α i (s) and β i (s) for many values of i at the same time. The same values of δ i2 j+1 ,i2 j+1 +2 j (s ′ , s) used to compute each α i (s) are used to compute each β i (s).
In this section, we describe a method to compute the values of α i (s), β i (s) for many stages i simultaneously. The goal is to break the dependency of α i+1 (s) on α i (s) and of β i-i (s) on β i (s) and enable most of the speedup that parallel processing hardware would provide if these values could be computed independently.
The proposed method consists of a sequence of trellis reduction steps. These steps are shown in Fig. 5. At each step values of δ i,i+1 (s ′ , s) of a "reduced" trellis, that skips every other stage, are computed from pairs of these values of an original trellis. First the values of δ i,i+1 (s ′ , s) are computed as follows:
Pr(
where j can be any single value from 0, 1, . . . , k -1. Given the values of δ i,i+1 (s ′ , s) for all i and all s of an original trellis, the values of δ j,j+2 (s ′ , s) of a new trellis that joins every two stages of the original trellis may be obtained by the following recursion:
Equation ( 68) provides a method to compute the δ i,i+2 (s, s ′ ) given the values of deltas δ i,i+1 (s, s ′ ) and δ i+1,i+2 (s, s ′ ) of an original trellis. The initial values of δ i,i+1 (s, s ′ ) are obtained from Eq. ( 62) and Eq. ( 68). Then the values of δ i,i+2 (s, s ′ ) are obtained using Eq. ( 68), which become the values of δ i,i+1 (s, s ′ ) on a new trellis with half the stages.
To avoid describing the trellis we are operating on, we will define by δ i,i+2 (s ′ , s) the values of δ i,i+1 (s ′ , s) on the new trellis obtained after the first trellis reduction step, then define δ i,i+4 (s, s ′ ) the values of the new trellis after the second reduction step, and so on. Note that the values of i must be a multiple of number stages of the original trellis that are combined in the most recent "new trellis," and must be a power of 2.
For ease of notation suppose that n = 1 and the number of trellis stages is a power of 2 given by N = 2 L+1 for some L.
The method to construct all values of δ i,j (s ′ , s) of the form δ 2 l i,2 l i+2 l (s ′ , s) is illustrated in the trees at the top of Fig. 6 and Fig. 7, and is summarized in the following equations: The steps to construct such δ i,j (r, s ′ ) where i = 2 l and j = 2 l-1 (i + 1) are as follows:
. . . . . .
The α i (s) and β i (s), i = 1, 2, . . . , 2 L+1 , s = 1, 2, . . . , M may also be computed are also computed in a logarithmic number of steps, using the values of δ 2 l i,2 l i+2 l (s ′ , s), l = 0, 1, . . . , L, once the values of δ 0,2 L(s ′ , s) and δ 2 L ,2 L+1(s ′ , s) are obtained. This process is illustrated in the bottom trees of Fig. 6 and Fig. 7. To start set
s, 0), then, for each j = 0, 1, . . . , L -1 and each i = 0, 1, . . . , 2 L-j compute:
With Eq. (72) and Eq. ( 73) the values of α 2 j+1 -2 j (s) and β i2 j+1 -2 j (s) for fixed j can be computed in parallel, starting from j = L, L -1, . . . , 1, where L = log 2 (N ). This method allows to compute the β i (s) and β i (s) from which future values in a logarithmic number of given
3.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4.5 K s (dB) 2 4 6 8 10 12 14 16 20 22 24 26 28 30 32 34 Average iterations GPU with 15120 stages CPU with 15120 stages GPU with 2 14 = 16384 stages CPU with 2 14 = 16384 stages by L = log 2 (N ).
The same holds for the δ i2 j ,(i+1)2 j (s ′ , s) values.
In the case where the code length is not a power of 2, like the case of the SCPPM code with standard length N = 15120, two sets of δ 2 l i,2 l i+2 l (s ′ , s) may need to be computed, one to compute all α i (s) and one to compute all β i (s). The set to compute the values of α i (s) starts at the first stage, with i = 0, 1, . . . , 2 L-j for fixed j, and the set to compute the values of β i (s) starts at the last stage N . That is, computing instead δ N -2 l i,N -2 l i+2 l (s ′ , s) for i = 0, 1, . . . , 2 L-j for fixed j. However, if N is multiple of a power of 2, say 2 l with l < L, then majority of the values of δ 2 j i,2 j i+2 j (s ′ , s) may be re-used, that is, all for which j ≤ l.
The Graphics Processing Unit is a massively parallel computing device with thread numbers rising to thousands. Thread is a single unit of execution that completes instructions sequentially. While threads are not completely independent, separate cores use the same instructions making development easier.
We choose CUDA as it allows control of individual threads in contrast to other programming platforms supporting GPU such as PyTorch. Synchronization may be performed even on a subset of threads and memory is manually managed. Moreover, using CUDA's finer memory management interface we utilized data locality in trellis reduction operations.
GPU groups threads into Streaming Multiprocessors (SMs) that can be used to process relevant tasks using the same data. Additionally, synchronization within an SM is faster, compared to synchronization of all threads on a GPU. While "global" memory is accessible from all threads, each SM contains an additional faster "shared" memory accessible only from threads within an SM. Another memory level, "local", is only available within a thread and is the fastest. The whole workload is divided among available SMs while only the most accessed data is stored at local and shared
A naive implementation of the trellis reduction algorithm, as described in Algorithm 1, stores all data in global memory and assigns each thread a stage to process. As global memory is significantly larger than the size of all required variables, this approach is easily scaled up for larger number of stages, but runtime is limited by the memory bandwidth.
The memory hierarchy suggests further optimizations. Large arrays with rare accesses including the α i (s) and β i (s) values are located in global memory, while the γ i (e) values could remain local to the thread as the values are unused by other threads. On the other hand, the values of δ i,j (s ′ , s) are both compute-intensive and used by multiple threads. We place it in shared memory by exploiting the tree structure of trellis.
The upper part of trellis as seen in Fig. 6 is divided into blocks equal to the desired amount of SMs. Each block uses only its own shared memory. This approach not only improves memory access time but also removes the need for synchronization between SMs during trellis processing, as separate branches of the tree can be computed independently.
Shared memory usage acts as the main constraint for the implementation on modern Nvidia GPU architectures. The data structure containing the values of δ i,j (s ′ , s) is proportional in size to the code length and is divided between SMs. Given a GPU with enough SMs, any code length can be decoded using our implementation.
To verify the functionality of the proposed GPU implementation we compare GPU and baseline CPU implementations in terms of codeword error rate (CWER) performance. For the GPU implementation, we used 32 SMs with 256 threads each. CWER is measured by running each simulation until 100 codeword errors are observed (except for the last point in the GPU with N = 2 14 stages and N = 15120, where the number of errors is about 50). We run the original BCJR algorithm on the CPU and the trellis reduction algorithm on the GPU using CUDA as described in Section Performance results in of codeword error rate vs. average photons per pulsed slots K s (dB) are presented in Fig. 8. The data was obtained with K b = 0 dB and with 4 PPM. The curves are essentially identical, verifying the equivalent functioning of CPU and GPU decoders. Using this simulation framework, Fig. 9 shows the average number of iterations required to decode a single codeword. The maximum number of iterations permitted is 32. As the noise is decreases, the expected number of iterations also decreases. The curves further confirm that the CPU and GPU implementations require the same number of iterations to achieve decoding. Thus, the runtime of a given number of iterations can be used to estimate runtime for a given K s value using either Fig. 10 or Fig. 11.
To evaluate the runtime benefits of the GPU implementation over the CPU implementation, we decode 1000 codewords with 10 iterations required to decode a single codeword on average. In total 10000 iterations are performed to obtain a single data point. Results are gathered for GPU and CPU implementations and shown in Fig. 10 and Fig. 11. These graphs show the same data; one uses a linear scale for the Y-axis and the other uses a log scale. We used powers of 2 for the values of N to demonstrate how runtime increases as a function of N for these two implementations. Both figures show the significant decrease in runtime achieved by the GPU implementation as compared with the CPU implementation. More stages allow for more parallelism so that the GPU runtime remains relatively constant even as N increases as long as there are sufficient hardware resources. The exact difference can be seen on the log scale graph in Fig. 11. Even for 2 8 = 256 stages, GPU provides 38% speed-up. The benefit rises to 15 times runtime reduction with 2 14 = 16384 stages.
This work shows an efficient method to execute in parallel many of the computations of the BCJR algorithm that are generally computed sequentially. Provided sufficient hardware resource, a GPU can effectively use the benefits of the trellis reduction algorithm to significantly reduce the decoding time of each codeword. The GPU implementation decode also offers greater portability than other hardware implementations like an FPGA. The runtime of the GPU implementation is improved up to an order of magnitude compared to CPU when the blocklength reaches 2 14 -bits.
The current version of the algorithm relies on normalization of the values to maintain numerical stability. The log domain version of the algorithm has better stability properties, but the latest implementation does not attain the run-time performance of the original version. In the future, we plan to explore optimizations of the log domain implementation and better memory management to improve the run-time performance. Further improvement in the run-time performance may be attained by an integer version of the log domain algorithm, aided by look-up tables.
Dariush Divsalar received the Ph.D. degree in electrical engineering from UCLA, in 1978. Since then, he has been with the Jet Propulsion Laboratory (JPL), California Institute of Technology (Caltech), Pasadena, where he is a Fellow. At JPL, he has been involved with developing state-of-the-art technology for advanced deep-space communications systems and future NASA space exploration. Since 1986, he has taught graduate courses in communications and coding at UCLA and Caltech. He