Viterbi译码程序代码

合集下载

Viterbi译码的MATLAB仿真研究

Viterbi译码的MATLAB仿真研究

BUPT卷积码编码及Viterbi译码班级:07114学号:070422姓名:吴希龙指导老师:彭岳星邮箱:FusionBupt@1. 序言卷积码最早于1955年由Elias 提出,稍后,1957年Wozencraft 提出了一种有效地译码方法即序列译码。

1963年Massey 提出了一种性能稍差但是比较实用的门限译码方法,使得卷积码开始走向实用化。

而后1967年Viterbi 提出了最大似然译码算法,它对存储级数较小的卷积码很容易实现,被称作Viterbi 译码算法,广泛的应用于现代通信中。

2. 卷积码编码及译码原理2.1 卷积码编码原理卷积码是一种性能优越的信道编码,它的编码器和解码器都比较易于实现,同时还具有较强的纠错能力,这使得它的使用越来越广泛。

卷积码一般表示为(n,k,K)的形式,即将k 各信息比特编码为n 个比特的码组,K 为编码约束长度,说明编码过程中相互约束的码段个数。

卷积码编码后的n 各码元不经与当前组的k 个信息比特有关,还与前K-1个输入组的信息比特有关。

编码过程中相互关联的码元有K*n 个。

R=k/n 是编码效率。

编码效率和约束长度是衡量卷积码的两个重要参数。

典型的卷积码一般选n,k 较小,但K 值可取较大(>10),以获得简单而高性能的卷积码。

卷积码的编码描述方式有很多种:冲激响应描述法、生成矩阵描述法、多项式乘积描述法、状态图描述,树图描述,网格图描述等。

2.1.1 卷积码解析表示法卷积码的解析表示发大致可以分为离散卷积法,生成矩阵法,码多项式法。

下面以离散卷积为例进行说明。

卷积码的编码器一般比较简单,为一个具有k 个输入端,n 个输出端,m 级移位寄存器的有限状态有记忆系统。

下图所示为(2,1,7)卷积码的编码器。

若输入序列为u =(u 0u 1u 2u 3……),则对应两个码字序列c ①=(c 0①c 1①c 2①c 3①……)和c ②=(c 0②c 1②c 2②c 3②……)相应的编码方程可写为c ①=u ∗g ①,c ②=u ∗g ②,c=(c ①,c ②)。

FEC前向纠错,卷积编码之维特比译码

FEC前向纠错,卷积编码之维特比译码

FEC前向纠错,卷积编码之维特⽐译码因为要学习做WCDMA的流程解析,需要先提取卷积数据,⾸先就要做FEC卷积译码。

于是⽹上翻了好⼤⼀圈,特地学习了下viterbi译码算法,费很⼤⼒⽓才凑齐能够正确跑起来的代码,特记录⼀下。

说点题外话:viterbi是个⼈,全名Andrew J. Viterbi,⼀枚数学家,美国⾼通公司的创始⼈之⼀(没错,就是现在⼿机上那个⾼通芯⽚的⾼通),现代⽆线数字通讯⼯程的奠基⼈。

此外,viterbi算法不但可⽤来做卷积译码,另外⼀个⼴泛的应⽤是做股票证券量化预测计算(隐马模型预测)。

没错,就是⾼频量化交易炒股的算法基础。

美国量化对冲基⾦的传奇,⼤奖章基⾦背后的⽼板James Simons,⾸创将隐马模型预测⽤于炒股,也是⼀枚数学家。

两枚玩数学算法的巨富,谁说学数学算法没⽤的?FEC前向纠错译码⽤的viterbi算法代码摘抄⾃libfec,⾥⾯的实现很全⾯,还有mmx,sse加速版。

测试是做WCDMA,所以解交织和译码⽤的3GPP WCDMA的卷积编码参数(⼀帧编码长度262bits,卷积编码多项式是 r=1/2 k=9),我就只摘取了需要的viterbi29部分,有需要全部实现的同学可以⾃⾏去下载libfec。

下⾯是写来⽤于测试验证算法正确性的线程代码。

从IQ⽂件中⼀次读取2个float数,读取262+8(⼀帧viterbi译码长度)个数据后进⾏译码,若译码成功程序停在对应的测试断点上。

UINT CALLBACK ProcessIQDataThread(void *p_param){Convolution2Viterbi_t *p_sys = (Convolution2Viterbi_t *)p_param;vector<double> *p_data = NULL;int framebits = 262;struct v29 *vp;vector<float> uninter1;vector<float> uninter2;vector<float> radioFrame;vector<float> frameL;vector<float> frameR;vector<uint8_t> viterbi_in;uint8_t viterbi_out[33];uninter2.resize((framebits + 8));uninter1.resize((framebits + 8) * 2);viterbi_in.resize((framebits + 8) * 2);if ((vp = (struct v29 *)create_viterbi29_port(framebits)) == NULL) {return -1;}while (!p_sys->exit_thread){if (!p_data){WaitForSingleObject(p_sys->iq_data_event, 200);p_sys->iq_data_list.Lock();p_data = p_sys->iq_data_list.GetUse();p_sys->iq_data_list.UnLock();}if (!p_data)continue;// 读取2个浮点数for (int i = 0; i < p_data->size() / 2; i++){frameL.push_back(p_data->at(0));frameR.push_back(p_data->at(1));}// ⼀帧译码长度if (frameL.size() >= (framebits + 8)){// 解交织inter2deInterleavingNP(30, inter2Perm, frameL, uninter2);radioFrame.insert(radioFrame.end(), uninter2.begin(), uninter2.end());deInterleavingNP(30, inter2Perm, frameR, uninter2);radioFrame.insert(radioFrame.end(), uninter2.begin(), uninter2.end());if (radioFrame.size() >= (framebits + 8) * 2){// 解交织inter1deInterleavingNP(inter1Columns[1], inter1Perm[1], radioFrame, uninter1);radioFrame.clear();// 使⽤Viterbi算法做FEC译码for (int i = 0; i < (framebits + 8) * 2; i++){viterbi_in.push_back( ((int)uninter1[i] >> 24) );}init_viterbi29_port(vp, 0);update_viterbi29_blk_port(vp, viterbi_in.data(), framebits + 8);chainback_viterbi29_port(vp, viterbi_out, framebits, 0);#ifdef _DEBUG// 数据验证for (int i = 0; i < 33 - 4; i++) {if (viterbi_out[i] == 0x98 && viterbi_out[i + 1] == 0x54 && viterbi_out[i + 2] == 0xA0 && viterbi_out[i + 3] == 0x00) int bbb = 0;if (viterbi_out[i] == 0x05 && viterbi_out[i + 1] == 0x33 && viterbi_out[i + 2] == 0x00 && viterbi_out[i + 3] == 0x0c) int bbb = 0;if (viterbi_out[i] == 0x00 && viterbi_out[i + 1] == 0x03 && viterbi_out[i + 2] == 0xf4 && viterbi_out[i + 3] == 0x40) int bbb = 0;if (viterbi_out[i] == 0xc5 && viterbi_out[i + 1] == 0x80 && viterbi_out[i + 2] == 0x05 && viterbi_out[i + 3] == 0x5a) int bbb = 0;if (viterbi_out[i] == 0xe4 && viterbi_out[i + 1] == 0x00 && viterbi_out[i + 2] == 0x33 && viterbi_out[i + 3] == 0xe6) int bbb = 0;if (viterbi_out[i] == 0x72 && viterbi_out[i + 1] == 0x08 && viterbi_out[i + 2] == 0x38)int bbb = 0;if (viterbi_out[i] == 0xe2 && viterbi_out[i + 1] == 0x00 && viterbi_out[i + 2] == 0x38)int bbb = 0;}#endif}frameL.clear();frameR.clear();}p_sys->iq_data_list.Lock();p_sys->iq_data_list.PushEmpty(p_data);p_data = p_sys->iq_data_list.GetUse();p_sys->iq_data_list.UnLock();}return 0;}解交织 Interleaving.h#pragma once#include <vector>#include <assert.h>const int inter1Columns[4] = {1, // TTI = 10ms2, // TTI = 20ms4, // TTI = 40ms8 // TTI = 80ms};// 25.212 4.2.5.2 table 4: Inter-Column permutation pattern for 1st interleaving:const char inter1Perm[4][8] = {{0}, // TTI = 10ms{0, 1}, // TTI = 20ms{0, 2, 1, 3}, // TTI = 40ms{0, 4, 2, 6, 1, 5, 3, 7} // TTI = 80ms};// 25.212 4.2.11 table 7: Inter-Column permutation pattern for 2nd interleaving:const char inter2Perm[30] = { 0,20,10,5,15,25,3,13,23,8,18,28,1,11,21,6,16,26,4,14,24,19,9,29,12,2,7,22,27,17 };/* FIRST steps two pointers through a mapping, one pointer into the interleaved* data and the other through the uninterleaved data. The fifth argument, COPY,* determines whether the copy is from interleaved to uninterleaved, or back.* FIRST assumes no padding is necessary.* The reason for the define is to minimize the cost of parameterization and* function calls, as this is meant for L1 code, while also minimizing the* duplication of code.*/#define FIRST(UNINTERLEAVED,UNINTERLEAVEDP,INTERLEAVED,INTERLEAVEDP,COPY) \assert(UNINTERLEAVED.size() == INTERLEAVED.size()); \unsigned int rows = UNINTERLEAVED.size() / columns; \assert(rows * columns == UNINTERLEAVED.size()); \const char *colp = permutation; \float *INTERLEAVEDP = &INTERLEAVED[0]; \for (unsigned i = 0; i < columns; i++) { \float *UNINTERLEAVEDP = &UNINTERLEAVED[*colp++]; \for (unsigned int j = 0; j < rows; j++) { \COPY; \UNINTERLEAVEDP += columns; \} \}/** interleaving with No Padding */void interleavingNP(const unsigned int columns, const char *permutation, vector<float> &in, vector<float> &out) {FIRST(in, inp, out, outp, *outp++ = *inp)}/** de-interleaving with No Padding */void deInterleavingNP(const unsigned int columns, const char *permutation, vector<float> &in, vector<float> &out) {FIRST(out, outp, in, inp, *outp = *inp++)}/* SECOND steps two pointers through a mapping, one pointer into the interleaved* data and the other through the uninterleaved data. The fifth argument, COPY,* determines whether the copy is from interleaved to uninterleaved, or back.* SECOND pads if necessary.* The reason for the define is to minimize the cost of parameterization and* function calls, as this is meant for L1 code, while also minimizing the* duplication of code.*/#define SECOND(UNINTERLEAVED,UNINTERLEAVEDP,INTERLEAVED,INTERLEAVEDP,COPY) \assert(UNINTERLEAVED.size() == INTERLEAVED.size()); \int R2 = (UNINTERLEAVED.size() + columns - 1) / columns; \int padding = columns * R2 - UNINTERLEAVED.size(); \int rows = R2; \int firstPaddedColumn = columns - padding; \const char *colp = permutation; \float *UNINTERLEAVEDP = &UNINTERLEAVED[0]; \for (int i = 0; i < columns; i++) { \int trows = rows - (*colp >= firstPaddedColumn); \float *INTERLEAVEDP = &INTERLEAVED[*colp++]; \for (int j = 0; j < trows; j++) { \COPY; \INTERLEAVEDP += columns; \} \}/** interleaving With Padding */void interleavingWP(const int columns, const char *permutation, vector<float> &in, vector<float> &out){SECOND(in, inp, out, outp, *outp = *inp++)}/** de-interleaving With Padding */void deInterleavingWP(const int columns, const char *permutation, vector<float> &in, vector<float> &out){SECOND(out, outp, in, inp, *outp++ = *inp)}/*Determining the constants.From the standard we know:* All frame sizes for the BCH.* transport block is 246 bits* there are two radio frames, 270 bits each* TTI is 20 ms* SF is 256* parity word Li is 16 bits* For all downlink TrCH except BCH, the radio frame size is 2*38400/SF = 76800/SF.* For SF=256 that's 300.* For SF=4 that's 19200.* The maximum code block size for convulutional coding is 540 bits (25.212 4.2.2.2).* That corresponds to a radio frame size of 1080 bits, or a spreading factor of 71,meaning that the smallest spreading factor that can be used is 128.* 76800/128 = 600 coded bits -> roughly 300 data bits.* That corresponds to an input rate of roughly 30 kb/s at.* The maximum code block size for turbo coding is 5114 bits (25.212 4.2.2.2).* That corresponds to a radio frame size of 15342 bits, or a spreading factor of 5,meaning that the smallest spreading factor that can be used is 8.* 76800/8 = 9600 coded bits -> roughly 3200 data bits.* That corresponds to an input rate of roughly 320 kb/s.OK - SO HOW DO YOU GET HIGHER RATES?? HOW CAN YOU USE SF=4??A: Use the full 5114-but code block and then expand it with rate-matching.You still can't get the full ~640 kb/s implied by SF=4, but you get to ~500 kb/s.(pat) A: They considered this problem. See 25.212 4.2.2.2 Code block segmentation.In Layer1, after transport block concatenation, you then simply chop the result upinto the largest pieces that can go through the encoder, then put them back together after. From "higher layers" we are given:* SF: 4, 8, 16, 32, 64, 128, 256.* P: 24, 16, 8, 0 bits* TTI: 10, 20, 40 ms.To simplify things, we set:* TTI 10 ms always on DCH and FACH, 20 ms on PCH and BCH* BCH and PCH are always rate-1/2 convolutional code* DCH and FACH are always rate-1/3 turbo code* no rate-matching, no puncturing* parity word is always 16 bits* So the only parameter than changes is spreading factor.* We will only support 15-slot (non-compressed) slot formats.From our simplifications we also know:* For non-BCH/PCH TrCH there is one radio frame,76800/SF channel (coded) bits, per transport block.* DCH and FACH always use rate-1/3 turbo code,which has 12 terminating bits in the output.* For DCH and FACH, the transport block size is((76800/SF - 12)/3) - P = (25600/SF) - 4 - P data bits,where P is the parity word size.* Fix P=16 for simplicity and transport block size is (25600/SF) - 20.* for SF=256, that's 80 bits.* for SF=16, that's 1580 bits.* for SF=8, that's 3180 bits.* For PCH there is one radio frame,76800/SF channel (coded) bits, per transport block.* SF=64, for that's 1200 channel bits.* It's a rate-1/2 conv code, so that's 1200/2 - 8 - P data bits.* P=16 so that's 1200/2 - 24 = 576 transport bits. Really?*/#if 0const int inter1Columns[] = { 1, 2, 4, 8 };const char inter1Perm[4][8] = {{0},{0, 1},{0, 2, 1, 3},{0, 4, 2, 6, 1, 5, 3, 7}};const char inter2Perm[] = {0, 20, 10, 5, 15, 25, 3, 13, 23, 8, 18, 28, 1, 11, 21,6, 16, 26, 4, 14, 24, 19, 9, 29, 12, 2, 7, 22, 27, 17};vector<char> randomBitVector(int n){vector<char> t(n);for (int i = 0; i < n; i++) t[i] = random() % 2;return t;}void testInterleavings(){int lth1 = 48;int C2 = 30;for (int i = 0; i < 4; i++) {vector<char> v1 = randomBitVector(lth1);vector<char> v2(lth1);vector<char> v3(lth1);v1.interleavingNP(inter1Columns[i], inter1Perm[i], v2);v2.deInterleavingNP(inter1Columns[i], inter1Perm[i], v3);cout << "first " << i << " " << (veq(v1, v3) ? "ok" : "fail") << endl;}for (int lth2 = 90; lth2 < 120; lth2++) {vector<char> v1 = randomBitVector(lth2);vector<char> v2(lth2);vector<char> v3(lth2);v1.interleavingWP(C2, inter2Perm, v2);v2.deInterleavingWP(C2, inter2Perm, v3);cout << "second " << lth2 << " " << (veq(v1, v3) ? "ok" : "fail") << endl;}for (int lth = 48; lth <= 4800; lth *= 10) {TurboInterleaver er(lth);cout << "Turbo Interleaver permutation(" << lth << ") " << (permutationCheck(er.permutation()) ? "ok" : "fail") << endl; vector<char> er1 = randomBitVector(lth);vector<char> er2(lth);er.interleave(er1, er2);vector<char> er3(lth);er.unInterleave(er2, er3);cout << "Turbo Interleaver(" << lth << ") " << (veq(er1.sliced(), er3) ? "ok" : "fail") << endl;}}#endifviterbi译码partab实现 partab.c/* Utility routines for FEC support* Copyright 2004, Phil Karn, KA9Q*/#include <stdio.h>#include "fec.h"unsigned char Partab[256] = {0};int P_init = 0;/* Create 256-entry odd-parity lookup table* Needed only on non-ia32 machines*/void partab_init(void) {int i, cnt, ti;/* Initialize parity lookup table */for (i = 0; i < 256; i++) {cnt = 0;ti = i;while (ti) {if (ti & 1)cnt++;ti >>= 1;}Partab[i] = cnt & 1;}P_init = 1;}/* Lookup table giving count of 1 bits for integers 0-255 */int Bitcnt[] = {0, 1, 1, 2, 1, 2, 2, 3,1, 2, 2, 3, 2, 3, 3, 4,1, 2, 2, 3, 2, 3, 3, 4,2, 3, 3, 4, 3, 4, 4, 5,1, 2, 2, 3, 2, 3, 3, 4,2, 3, 3, 4, 3, 4, 4, 5,2, 3, 3, 4, 3, 4, 4, 5,3, 4, 4, 5, 4, 5, 5, 6,1, 2, 2, 3, 2, 3, 3, 4,2, 3, 3, 4, 3, 4, 4, 5,2, 3, 3, 4, 3, 4, 4, 5,3, 4, 4, 5, 4, 5, 5, 6,2, 3, 3, 4, 3, 4, 4, 5,3, 4, 4, 5, 4, 5, 5, 6,3, 4, 4, 5, 4, 5, 5, 6,4, 5, 5, 6, 5, 6, 6, 7,1, 2, 2, 3, 2, 3, 3, 4,2, 3, 3, 4, 3, 4, 4, 5,2, 3, 3, 4, 3, 4, 4, 5,3, 4, 4, 5, 4, 5, 5, 6,2, 3, 3, 4, 3, 4, 4, 5,3, 4, 4, 5, 4, 5, 5, 6,3, 4, 4, 5, 4, 5, 5, 6,4, 5, 5, 6, 5, 6, 6, 7,2, 3, 3, 4, 3, 4, 4, 5,3, 4, 4, 5, 4, 5, 5, 6,3, 4, 4, 5, 4, 5, 5, 6,4, 5, 5, 6, 5, 6, 6, 7,3, 4, 4, 5, 4, 5, 5, 6,4, 5, 5, 6, 5, 6, 6, 7,4, 5, 5, 6, 5, 6, 6, 7,5, 6, 6, 7, 6, 7, 7, 8,};fec.h/* User include file for libfec* Copyright 2004, Phil Karn, KA9Q* May be used under the terms of the GNU Lesser General Public License (LGPL)*/#ifndef _FEC_H_#define _FEC_H_#ifdef __cplusplusextern "C" {#endif/* r=1/2 k=9 convolutional encoder polynomials */#define V29POLYA 0x1af#define V29POLYB 0x11dvoid *create_viterbi29_port(int len);void set_viterbi29_polynomial_port(int polys[2]);int init_viterbi29_port(void *p, int starting_state);int chainback_viterbi29_port(void *p, unsigned char *data, unsigned int nbits, unsigned int endstate); void delete_viterbi29_port(void *p);int update_viterbi29_blk_port(void *p, unsigned char *syms, int nbits);void partab_init();static inline int parityb(unsigned char x) {extern unsigned char Partab[256];extern int P_init;if (!P_init) {partab_init();}return Partab[x];}static inline int parity(int x) {/* Fold down to one byte */x ^= (x >> 16);x ^= (x >> 8);return parityb(x);}#ifdef __cplusplus}#endif#endif /* _FEC_H_ */viterbi29_port.c/* K=9 r=1/2 Viterbi decoder in portable C* Copyright Feb 2004, Phil Karn, KA9Q* May be used under the terms of the GNU Lesser General Public License (LGPL)*/#include <stdio.h>#include <stdlib.h>#include <memory.h>#include "fec.h"typedef union { unsigned int w[256]; } metric_t;typedef union { unsigned long w[8]; } decision_t;static union { unsigned char c[128]; } Branchtab29[2];static int Init = 0;/* State info for instance of Viterbi decoder */struct v29 {metric_t metrics1; /* path metric buffer 1 */metric_t metrics2; /* path metric buffer 2 */decision_t* dp; /* Pointer to current decision */metric_t* old_metrics, * new_metrics; /* Pointers to path metrics, swapped on every bit */decision_t* decisions; /* Beginning of decisions for block */};/* Initialize Viterbi decoder for start of new frame */int init_viterbi29_port(void* p, int starting_state) {struct v29* vp = p;int i;if (p == NULL)return -1;for (i = 0; i < 256; i++)vp->metrics1.w[i] = 63;vp->old_metrics = &vp->metrics1;vp->new_metrics = &vp->metrics2;vp->dp = vp->decisions;vp->old_metrics->w[starting_state & 255] = 0; /* Bias known start state */return 0;}void set_viterbi29_polynomial_port(int polys[2]) {int state;for (state = 0; state < 128; state++) {Branchtab29[0].c[state] = (polys[0] < 0) ^ parity((2 * state) & abs(polys[0])) ? 255 : 0; Branchtab29[1].c[state] = (polys[1] < 0) ^ parity((2 * state) & abs(polys[1])) ? 255 : 0; }Init++;}/* Create a new instance of a Viterbi decoder */void* create_viterbi29_port(int len) {struct v29* vp;if (!Init) {int polys[2] = { V29POLYA,V29POLYB };set_viterbi29_polynomial_port(polys);}if ((vp = (struct v29*)malloc(sizeof(struct v29))) == NULL)return NULL;if ((vp->decisions = (decision_t*)malloc((len + 8) * sizeof(decision_t))) == NULL) {free(vp);return NULL;}init_viterbi29_port(vp, 0);return vp;}/* Viterbi chainback */int chainback_viterbi29_port(void* p,unsigned char* data, /* Decoded output data */unsigned int nbits, /* Number of data bits */unsigned int endstate) /* Terminal encoder state */{struct v29* vp = p;decision_t* d;if (p == NULL)return -1;d = vp->decisions;/* Make room beyond the end of the encoder register so we can* accumulate a full byte of decoded data*/endstate %= 256;/* The store into data[] only needs to be done every 8 bits.* But this avoids a conditional branch, and the writes will* combine in the cache anyway*/d += 8; /* Look past tail */while (nbits-- != 0) {int k;k = (d[nbits].w[(endstate) / 32] >> (endstate % 32)) & 1;data[nbits >> 3] = endstate = (endstate >> 1) | (k << 7);}return 0;}/* Delete instance of a Viterbi decoder */void delete_viterbi29_port(void* p) {struct v29* vp = p;if (vp != NULL) {free(vp->decisions);free(vp);}}/* C-language butterfly */#define BFLY(i) {\unsigned int metric,m0,m1,decision;\metric = (Branchtab29[0].c[i] ^ sym0) + (Branchtab29[1].c[i] ^ sym1);\ m0 = vp->old_metrics->w[i] + metric;\m1 = vp->old_metrics->w[i+128] + (510 - metric);\decision = (signed int)(m0-m1) > 0;\vp->new_metrics->w[2*i] = decision ? m1 : m0;\d->w[i/16] |= decision << ((2*i)&31);\m0 -= (metric+metric-510);\m1 += (metric+metric-510);\decision = (signed int)(m0-m1) > 0;\vp->new_metrics->w[2*i+1] = decision ? m1 : m0;\d->w[i/16] |= decision << ((2*i+1)&31);\}/* Update decoder with a block of demodulated symbols* Note that nbits is the number of decoded data bits, not the number* of symbols!*/int update_viterbi29_blk_port(void* p, unsigned char* syms, int nbits) { struct v29* vp = p;decision_t* d;if (p == NULL)return -1;d = (decision_t*)vp->dp;while (nbits--) {void* tmp;unsigned char sym0, sym1;int i;for (i = 0; i < 8; i++)d->w[i] = 0;sym0 = *syms++;sym1 = *syms++;for (i = 0; i < 128; i++)BFLY(i);d++;tmp = vp->old_metrics;vp->old_metrics = vp->new_metrics;vp->new_metrics = tmp;}vp->dp = d;return 0;}。

简单的维特比译码

简单的维特比译码
M1=length(r);%译码码字长度
L1=M1/N1;%译码输出码长
%w=zeros(1,M1); %w表示判决输出;
%for mm=1:M1
% if abs(r(mm))>1
% w(mm)=1 ;
%else w(mm)=0;
% end
for k=1:L1 %输出序列码长
for st=1:numStates %状态
for i=0:numInputs-1 %当前输入
dout=dec2bin(trellis.outputs(st,i+1),N1)-48;
dist=sum((r((k-1)*N1+1:k*N1)-dout).^2); %计算汉明码距
next_metric(nextState)=x;
sur_path(nextState,1:k)=[path(st,1:k-1) i];
else
if x<next_metric(nextState) % 择优路径
cur_metric=zeros(1,numStates)+Inf;%汉明码距
cur_metric(1)=0;
next_metric=-ones(1,numStates);%下一状态码距
path=zeros(numStates,L1);
sur_path=zeros(numStates,L1);
next_metric=-ones(1,numStates);
path=sur_path; %最佳路径
end
out1=path(1,:); %输出
next_metric(nextState)=x;
sur_path(nextState,1:k)=[path(st,1:k-1) i];

维特比译码程序

维特比译码程序

(n,k,N)卷积码的维特比译码算法实现#include<iostream>#define t_src 0#define t_des 1#define t_len 2#define t_flag 3#define t_in 4using namespace std;int myn=0;int stalen=0;int myg1[10]={0};int myg2[10]={0};int stan0[256][2]={0};//输入0 时个状态的输出int stan1[256][2]={0};//输入1 时各状态的输出int stachn[256][2]={0};//状态装换表int path[256][100]={0};//存储路径int calpath[256]={0};//存储路径长度int myin[24]; //一次处理12 次int myout[200]; //int myoutsym=0;int pthsym;int outfull=0; //决定是否输出int table1[8]={1,2,4,8,16,32,64,128};void chartobits(char ch,int *bits);char bitstochar(int *bits);int calluj(int a1,int a2,int b1,int b2);void initpath(void);void selpath(int a1,int a2);void wridata(void);void viterbit(void);void writdataedn(void);void creatsta(void);void myinput(void);int main(){myinput();creatsta();viterbit();}void myinput(void){int i,j;输入编码的约束长度cin>>myn;stalen=int(pow(2.0,myn-1));选择默认的编码矢量则输入1,输入 2 则可输入其他的编码矢量cin>>i;if(i==1){switch(myn){case 3: myg1[0]=1,myg1[1]=1,myg1[2]=1;myg2[0]=1,myg2[1]=0,myg2[2]=1;break;case 4: myg1[0]=1,myg1[1]=1,myg1[2]=1,myg1[3]=1;myg2[0]=1,myg2[1]=0,myg2[2]=1,myg2[3]=1;break;case 5: myg1[0]=1,myg1[1]=0,myg1[2]=1,myg1[3]=1,myg1[4]=1;myg2[0]=1,myg2[1]=1,myg2[2]=0,myg2[3]=1,myg2[4]=1;break;case 6: myg1[0]=1,myg1[1]=0,myg1[2]=1,myg1[3]=1,myg1[4]=1,myg1[5]=1;myg2[0]=1,myg2[1]=1,myg2[2]=0,myg2[3]=1,myg2[4]=0,myg2[5]=1;break;case 7: myg1[0]=1,myg1[1]=0,myg1[2]=0,myg1[3]=1,myg1[4]=1,myg1[5]=1,myg1[6]=1;myg2[0]=1,myg2[1]=1,myg2[2]=0,myg2[3]=1,myg2[4]=1,myg2[5]=0,myg2[6]=1;break;case 8: myg1[0]=1,myg1[1]=0,myg1[2]=0,myg1[3]=1,myg1[4]=1,myg1[5]=1,myg1[6]=1,myg1[7]=1;myg2[0]=1,myg2[1]=1,myg2[2]=1,myg2[3]=0,myg2[4]=0,myg2[5]=1,myg2[6]=0,myg2[7]=1;break;case 9: myg1[0]=1,myg1[1]=1,myg1[2]=0,myg1[3]=1,myg1[4]=0,myg1[5]=1,myg1[6]=1,myg1[7]=1,m yg1[8]=1;myg2[0]=1,myg2[1]=0,myg2[2]=0,myg2[3]=0,myg2[4]=1,myg2[5]=1,myg2[6]=1,myg2[7]=0,m yg2[8]=1;break;}}else{输入for(j=0;j<myn;j++)cin>>myg1[j];输入for(j=0;j<myn;j++)cin>>myg2[j];}连接矢量1 为for(j=0;j<myn;j++)cout<<endl;连接矢量2 为for(j=0;j<myn;j++)cout<<endl;cout<<endl;}void creatsta(void){int i,j,k,myi;int tembits[10];for(i=0;i<stalen;i++){stan1[i][0]=0;stan1[i][1]=0;stan0[i][0]=0;stan0[i][1]=0;stachn[i][0]=i/2;myi=i;for(j=0;j<myn;j++){if(myi>=pow(2.0,myn-1-j)){tembits[j]=1;myi=myi-pow(2.0,myn-1-j);}else{tembits[j]=0;}}for(k=0;k<myn;k++){stan0[i][0]+=myg1[k]*tembits[k];stan0[i][1]+=myg2[k]*tembits[k];}stan0[i][0]=stan0[i][0]%2;stan0[i][1]=stan0[i][1]%2;myi=i+int(pow(2.0,myn-1));stachn[i][1]=myi/2;for(j=0;j<myn;j++){if(myi>=pow(2.0,myn-1-j)){tembits[j]=1;myi=myi-pow(2.0,myn-1-j);}else{tembits[j]=0;}}for(k=0;k<myn;k++){stan1[i][0]+=myg1[k]*tembits[k];stan1[i][1]+=myg2[k]*tembits[k];}stan1[i][0]=stan1[i][0]%2;stan1[i][1]=stan1[i][1]%2;}状态转移出for(i=0;i<stalen;i++)cout<<endl;输入0 状态转移后的输出for(i=0;i<stalen;i++)cout<<endl;输入1 状态转移后的输出for(i=0;i<stalen;i++)cout<<endl;}void chartobits(char ch,int *bits){int i;for(i=0;i<8;i++){if(ch<0)bits[i]=1;elsebits[i]=0;ch=ch<<1;}}char bitstochar(int *bits){char temp=0;int i;for(i=0;i<8;i++){if(bits[i]==1)temp+=table1[7-i];}return temp;}int calluj(int a1,int a2,int b1,int b2){int y=0;if(a1!=b1)y++;if(a2!=b2)y++;return(y);}void initpath(){int tem;int t_tem[256][5]={0};int i,j,k,l;int ljtem[256][100];int pttem[256]={0};int staflag[256]={0};staflag[0]=1;int a1,a2;for(l=0;l<myn-1;l++){for(i=0;i<stalen;i++)for(j=0;j<pthsym;j++)ljtem[i][j]=path[i][j];i=0;a1=myin[2*l];a2=myin[2*l+1];for(j=0;j<stalen;j++){if(staflag[j]==1){tem=calluj(a1,a2,stan0[j][0],stan0[j][1]);t_tem[i][t_src]=j;t_tem[i][t_des]=stachn[j][0];t_tem[i][t_len]=calpath[j]+tem;t_tem[i][t_in]=0;tem=calluj(a1,a2,stan1[j][0],stan1[j][1]);t_tem[i+1][t_src]=j;t_tem[i+1][t_des]=stachn[j][1];t_tem[i+1][t_len]=calpath[j]+tem;t_tem[i+1][t_in]=1;i=i+2;}}for(k=0;k<stalen;k++)staflag[k]=0;for(k=0;k<i;k++){staflag[t_tem[k][t_des]]=1;calpath[t_tem[k][t_des]]=t_tem[k][t_len];for(j=0;j<pthsym;j++)path[t_tem[k][t_des]][j]=ljtem[t_tem[k][t_src]][j];path[t_tem[k][t_des]][pthsym]=t_tem[k][t_in];}pthsym++;}初始化后的路径长度for(int i=0;i<8;i++)cout<<endl;*/}void selpath(int a1,int a2){//16 选8int t_tem[512][5]={0};int i,j,tem;int ljtem[256][100];j=0;for(i=0;i<2*stalen;i=i+2){tem=calluj(a1,a2,stan0[j][0],stan0[j][1]);t_tem[i][t_src]=j;t_tem[i][t_des]=stachn[j][0];t_tem[i][t_len]=calpath[j]+tem;t_tem[i][t_flag]=0;t_tem[i][t_in]=0;//t_tem[i][t_rep]=0;tem=calluj(a1,a2,stan1[j][0],stan1[j][1]);t_tem[i+1][t_src]=j;t_tem[i+1][t_des]=stachn[j][1];t_tem[i+1][t_len]=calpath[j]+tem;t_tem[i+1][t_flag]=0;t_tem[i+1][t_in]=1;//t_tem[i][t_rep]=0;j++;}for(i=0;i<2*stalen;i++)for(j=i+1;j<2*stalen;j++)if(t_tem[i][t_des]==t_tem[j][t_des]){if(t_tem[i][t_len]<=t_tem[j][t_len]){t_tem[i][t_flag]=1;}else{t_tem[j][t_flag]=1;}}for(i=0;i<stalen;i++)for(j=0;j<pthsym;j++)ljtem[i][j]=path[i][j];for(i=0;i<2*stalen;i++)if(t_tem[i][t_flag]==1){calpath[t_tem[i][t_des]]=t_tem[i][t_len];for(j=0;j<pthsym;j++)path[t_tem[i][t_des]][j]=ljtem[t_tem[i][t_src]][j];path[t_tem[i][t_des]][pthsym]=t_tem[i][t_in];}if(pthsym>16)outfull=1;pthsym++;}void wridata(){int i,j,icout,equcout;icout=0;equcout=0;for(i=0;i<pthsym;i++){for(j=0;j<stalen-1;j++)if(path[j][i]==path[j+1][i])equcout++;if(equcout==stalen-1){myout[myoutsym++]=path[0][i];icout++;equcout=0;}elsebreak;}if(icout!=0){for(i=0;i<pthsym-icout;i++){for(j=0;j<stalen;j++)path[j][i]=path[j][i+icout];}}pthsym=pthsym-icout;outfull=0;}void writdataedn(void){int i,j;i=0;for(j=1;j<stalen;j++)if(calpath[i]>calpath[j])i=j;for(j=0;j<pthsym;j++)myout[myoutsym++]=path[i][j]; }void viterbit(){FILE *fp_input,*fp_output;exit(0);}elseexit(0);}elsechar ch;int count=0;int i,j;char wch;int wcout=0;int mybit[8];ch=fgetc(fp_input);chartobits(ch,mybit);for(i=0;i<8;i++)myin[i]=mybit[i];while(feof(fp_input)==0){ch=fgetc(fp_input);输入输入数据1 为for(temi=0;temi<8;temi++)cout<<endl;*/if(count==0){chartobits(ch,mybit);for(i=0;i<8;i++)myin[i+8]=mybit[i];initpath();for(j=myn-1;j<8;j=j++)selpath(myin[2*j],myin[2*j+1]);}else{chartobits(ch,myin);for(j=0;j<4;j++)selpath(myin[2*j],myin[2*j+1]);}count++;if(count==0)count=1;//if(outfull==1)wridata();if(myoutsym>=8){wcout=int(myoutsym/8);for(i=0;i<wcout;i++){for(j=0;j<8;j++)mybit[j]=myout[8*i+j];wch=bitstochar(mybit);输出为fputc(wch,fp_output);}for(i=0;i<myoutsym-wcout*8;i++)myout[i]=myout[wcout*8+i];myoutsym=myoutsym-wcout*8;}}writdataedn();if(myoutsym>=3){for(i=0;i<8-myoutsym;i++)myout[myoutsym++]=0;wch=bitstochar(myout);fputc(wch,fp_output);}fclose(fp_input);fclose(fp_output);cin>>i;}。

卷积码的Viterbi译码设计设计

卷积码的Viterbi译码设计设计

摘要在数字通信系统中,通常采用差错控制编码来提高系统的可靠性。

自P.Elias 首次提出卷积码编码以来,这一编码技术至今仍显示出强大的生命力。

目前,卷积码已广泛应用在无线通信标准中,如GSM,CDMA2000和IS-95等无线通信标准中。

针对N-CDMA数据传输过程中的误码问题,本文论述了旨在提高数据传输质量的维特比译码器的设计。

虽然Viterbi译码复杂度较大,实现较为困难,但效率高,速度快。

因此本文着重分析和讨论了1/2速率的(2,1,9)卷积码编码和其Viterbi译码算法。

深入研究卷积码编码原理和Viterbi算法原理后,提出了(2,1,9)卷积码编码以及Viterbi算法的初始化、加—比—选和回溯设计方案,运用查表的方法,避免了大量繁琐计算,使得译码简洁迅速,译码器的实时性能良好。

并充分利用TMS320C54X系列DSP芯片,用汇编语言完成了(2,1,9)卷积码编码和Viterbi 译码的程序。

关键词:差错控制编码、卷积码、Viterbi译码、TMS320C54X、DSPAbstractIn digital communication systems, error control coding is usually used to improve system reliability. Since P.Elias put forward the convolutional coding the first time, the coding is still showing strong vitality.,has become widely used in satellite communications, wireless communications and many other communication systemsas a kind of channel coding method. such as GSM, CDMA2000 and has been a wireless communication standards of IS-95.In view of the error problem in the process of N-CDMA data transmission, this paper discusses the aims to improve the quality of data transmission of victor design than the decoder.Although Viterbi decoding complexity is bigger, more difficult to achieve, but high efficiency and fast speed. So this article emphatically analyzed and discussed the 1/2 rate (2,1,9) convolution code coding and its Viterbi decoding algorithm. In-depth study on principle of convolution code coding and Viterbi algorithm, proposed the convolution code coding and Viterbi algorithm (2,1,9) initialization, add - than - choose and back design, using look-up table method, to avoid a large amount of tedious calculation, the decoding and quick, good real-time performance of the decoder. Make full use of the series of TMS320C54X DSP chip, using assembly language to complete the(2,1,9)convolution code coding and Viterbi decoding process.Keywords: error control coding, convolutional code, Viterbi decoding, TMS320C54X目录摘要 (1)Abstract (2)目录 (3)1.绪论 (1)1.1 移动通信及N-CDMA背景 (1)1.2 数字通信概述 (1)1.3 卷积编码与译码的发展 (3)1.4 主要研究工作 (3)2.DSP与CCS简介 (5)2.1 DSP概述 (5)2.1.1 DSP的主要特点 (5)2.1.2 CSSU单元概述 (7)2.2 CCS概述 (8)2.3 本章小结 (8)3.卷积码的理论基础 (9)3.1 卷积码的概述 (9)3.1.1 卷积码基本原理 (9)3.1.2 卷积码的纠错能力 (9)3.1.3 卷积码的表示方法 (10)3.2 Viterbi译码的概述 (11)3.3 本章小结 (14)4.卷积编码的实现 (15)4.1 (2,1,9)卷积码编码 (15)4.1.1 (2,1,9)卷积码编码设计方案 (15)4.1.2 (2,1,9)卷积码编码流程图 (16)4.1.3 (2,1,9)卷积编码程序实现 (16)4.1.4 (2,1,9)的程序仿真 (17)4.2 (2,1,9)卷积码状态转换表 (17)4.2.1 (2,1,9)卷积码状态转换表的设计算法 (18)4.2.2 (2,1,9)卷积码状态转换表的流程图 (18)4.2.3 (2,1,9)卷积码状态表 (18)4.2.4 (2,1,9)卷积码状态表的蝶形结构 (21)4.3 本章小结 (22)5. Viterbi译码的实现 (23)5.1 Viterbi译码基础 (23)5.2 Viterbi译码算法 (23)5.3 变量定义情况 (25)5.4 初始化 (26)5.4.1 初始化流程图 (27)5.4.2 初始化程序仿真 (27)5.5 加-比-选 (28)5.5.1加-比-选流程图 (29)5.5.2加-比-选程序仿真 (30)5.6 回溯 (31)5.6.1 回溯流程图 (32)5.6.2 回溯仿真图 (33)5.7 Viterbi纠错测试 (34)5.8 本章小结 (34)总结 (36)致谢 ............................................................................ 错误!未定义书签。

Viterbi译码的Matlab实现

Viterbi译码的Matlab实现

2010年12月(上)Viterbi 译码的Matlab 实现张慧(盐城卫生职业技术学院,江苏盐城224006)[摘要]本文主要介绍了Viterbi 译码是一种最大似然译码算法,是卷积编码的最佳译码算法。

本文主要是以(2,1,2)卷积码为例,介绍了Viterbi 译码的原理和过程,并用Matlab 进行仿真。

[关键词]卷积码;Viterbi 译码1卷积码的类型卷积码的译码基本上可划分为两大类型:代数译码和概率译码,其中概率译码是实际中最常采用的卷积码译码方法。

2Viterbi 译码Viterbi 译码是由Viterbi 在1967年提出的一种概率译码,其实质是最大似然译码,是卷积码的最佳译码算法。

它利用编码网格图的特殊结构,降低了计算的复杂性。

该算法考虑的是,去除不可能成为最大似然选择对象的网格图上的路径,即,如果有两条路径到达同一状态,则具有最佳量度的路径被选中,称为幸存路径(surviving path )。

对所有状态都将进行这样的选路操作,译码器不断在网格图上深入,通过去除可能性最小的路径实现判决。

较早地抛弃不可能的路径降低了译码器的复杂性。

为了更具体的理解Viterbi 译码的过程,我们以(2,1,2)卷积码为例,为简化讨论,假设信道为BSC 信道。

译码过程的前几步如下:假定输入数据序列m ,码字U ,接收序列Z ,如图1所示,并假设译码器确知网格图的初始状态。

图1时刻t 1接收到的码元是11,从状态00出发只有两种状态转移方向,00和10,如图a 所示。

状态转换的分支量度是2;状态转换的分支径量度是0。

时刻t 2从每个状态出发都有两种可能的分支,如图b 所示。

这些分支的累积量度标识为状态量度┎a ,┎b ,┎c ,┎d ,与各自的结束状态相对应。

同样地,图c 中时刻t 3从每个状态出发都有两个分支,因此,时刻时到达每个状态的路径都有两条,这两条路径中,累积路径量度较大的将被舍弃。

如果这两条路径的路径量度恰好相等,则任意舍弃其中一条路径。

matlab卷积编码与viterbi译码的实现

matlab卷积编码与viterbi译码的实现

matlab卷积编码与viterbi译码的实现MATLAB中viterbi译码算法讨论⼤家可以再评论区交流!!!MATLAB中实现viterbi译码的函数为:convenc其中:code = convenc(msg,trellis)vitdec其中:vitdec(code,trellis,tblen,opmode,dectype)code卷积编码,trellis⽹格表,tblen回溯长度,opmode:cont、term、trunc,dectype:unquant、hard、soft;本⼈最近在做⼀个关于viterbi译码算法,最终在FPGA中实现,在FPGA中最终的实现⽅案为xillinx IP核实现。

在此之前⽤MATLAB进⾏仿真验证。

matlab程序:Tre = poly2trellis(7,[133 171]);通过poly2trellis⽣成逻辑关系图,如下图所⽰。

逻辑关系图%卷积编码:msg = [0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 0 0 1 0 1 0 0 1];code = convenc(msg,Tre);%code = [0,0,1,1,0,1,0,0,0,1,1,0,0,1,0,1,1,0,0,0,1,0,0,0,1,0,1,0,0,1,1,1,1,1,1,1,0,0,1,1,0,1,0,1,1,1,1,0];%这是通过convenc函数⽣成的卷积码%vitdec译码:%在vitdec译码过程中采⽤硬判决,通过不同的tblen和opmode来找出其中关系。

%(1) opmode = conttblen = 12;msg_dat = vitdec(code,Tre,tblen,'cont','hard');%msg_dat =[ 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 0 1 1 1 1 ];%通过了解到cont模式中,vitdec译码会有延迟,延迟的长度为tblen长度,所以在此对vitdec进⾏修改code_temp = [code,zeros(1,24)];msg_temp = vitdec(code_temp ,T,12,'cont','hard')msg_dat = msg_temp(13:end);%msg_dat = [ 0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 0 0 1 0 1 0 0 1];%此时vitdec译码出来的数据和信源⼀样tblen = 18;code_temp = [code,zeros(1,24)];msg_temp = vitdec(code_temp ,T,12,'cont','hard')msg_dat = msg_temp(13:end);%msg_dat = [ 0 0 0 0 0 0 0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 0 0];%此时vitdec译码出来的数据和信源在后⾯最后⼀位不⼀样%(2) opmode = termtblen = 12;msg_dat = vitdec(code,Tre,tblen,'term','hard');%msg_dat = [0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 1 0 0 0 0 0 0 0];%此时vitdec译码出来的数据和信源⼀样前16位和信源⼀样后⾯的就出错了tblen = 18;msg_dat = vitdec(code,Tre,tblen,'term','hard');%msg_dat = [0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 1 0 0 0 0 0 0 0];%此时vitdec译码出来的数据和信源⼀样前16位和信源⼀样后⾯的就出错了%(3)opmode = trunctblen = 12;msg_dat = vitdec(code,Tre,tblen,'trunc','hard');%msg_dat = [ 0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 0 0 1 0 1 0 0 1];%此时vitdec译码出来的数据和信源⼀样tblen = 18;msg_dat = vitdec(code,Tre,tblen,'trunc','hard');%msg_dat = [ 0 1 0 1 1 0 1 0 1 1 1 1 0 0 1 1 0 0 1 0 1 0 0 1];%此时vitdec译码出来的数据和信源⼀样总结:以上通过⽐较tblen和opmode模式的不同对产⽣的结果,其中cont和trunc的模式总结起来就是cont有tblen延迟,但是trunc没有。

Viterbi译码程序代码

Viterbi译码程序代码

译码主要部分#include"stdafx.h"//#define DEBUGvoid deci2bin(int d, int size, int *b);int bin2deci(int *b, int size);int nxt_stat(int current_state, int input, int *memory_contents);void init_quantizer(void);void init_adaptive_quant(float es_ovr_n0);int soft_quant(float channel_symbol);int soft_metric(int data, int guess);int quantizer_table[256];void sdvd(int g[2][K], float es_ovr_n0, long channel_length, float*channel_output_vector, int *decoder_output_matrix){int i, j, l, ll; //循环控制变量long t; //时间int memory_contents[K]; //记录输入内容int input[TWOTOTHEM][TWOTOTHEM]; //对当前状态以及下一个状态映射int output[TWOTOTHEM][2]; //卷积码编码输出矩阵int nextstate[TWOTOTHEM][2]; //下一个状态矩阵int accum_err_metric[TWOTOTHEM][2]; //误差累计矩阵int state_history[TWOTOTHEM][K * 5 + 1]; //历史状态表int state_sequence[K * 5 + 1]; //状态序列int *channel_output_matrix; //信道输出序列int binary_output[2];int branch_output[2]; //0或者1的输出分支int m, n, number_of_states, depth_of_trellis, step, branch_metric,sh_ptr, sh_col, x, xx, h, hh, next_state, last_stop;n = 2; //1/2为卷积码传输数据的码率m = K - 1;//寄存器个数number_of_states = (int)pow(2.0, m);//状态个数number of states = 2^(K - 1) = 2^mdepth_of_trellis = K * 5;for (i = 0; i < number_of_states; i++){for (j = 0; j < number_of_states; j++)input[i][j] = 0; //输入数组初始化for (j = 0; j < n; j++){nextstate[i][j] = 0;//下一个状态数组初始化output[i][j] = 0; //输出数组初始化}for (j = 0; j <= depth_of_trellis; j++){state_history[i][j] = 0;//历史状态数组初始化state_history[4][16] }accum_err_metric[i][0] = 0;//误差累计矩阵第一列初始化为0accum_err_metric[i][1] = MAXINT;//误差累计矩阵第二列初始化为一个很大的数}/*前向纠错简称FEC(Forward Error Correction),其原理是:发送方将要发送的数据附加上一定的冗余纠错码一并发送,接收方则根据纠错码对数据进行差错检测,如发现差错,由接收方进行纠正*//*产生状态转移矩阵、输出矩阵、输入矩阵*///输入矩阵表示的是FEC编码传输给下一个状态//下一个状态由输入和当前状态给出//输出矩阵for (j = 0; j < number_of_states; j++){for (l = 0; l < n; l++){next_state = nxt_stat(j, l, memory_contents);input[j][next_state] = l;/*计算给定的卷积编码器输出当前状态数和输入值*/branch_output[0] = 0;branch_output[1] = 0;for (i = 0; i < K; i++){branch_output[0] = branch_output[0] ^ memory_contents[i] & g[0][i];branch_output[1] = branch_output[1] ^ memory_contents[i] & g[1][i];}nextstate[j][l] = next_state;//下一个状态output[j][l] = bin2deci(branch_output, 2);//输出十进制}}#ifdef DEBUGprintf("\nInput:");for (j = 0; j < number_of_states; j++){printf("\n");for (l = 0; l < number_of_states; l++)printf("%2d ", input[j][l]);}printf("\nOutput:");for (j = 0; j < number_of_states; j++){printf("\n");for (l = 0; l < n; l++)printf("%2d ", output[j][l]);}printf("\nNext State:");for (j = 0; j < number_of_states; j++){printf("\n");for (l = 0; l < n; l++)printf("%2d ", nextstate[j][l]);}#endifchannel_output_matrix =(int *)malloc(channel_length * sizeof(int));if (channel_output_matrix == NULL){printf("allocation is failure!!\n");exit(1);}printf("\n");/*信道输出为n行,2列,每行对应于一个通道符号给定的位和每一列对应于一个已编码的位*/ channel_length = channel_length / n;init_adaptive_quant(es_ovr_n0);//进行优化,匹配过信噪比的量化矩阵//量化信道输出,将浮点型数据转化为整形for (t = 0; t < (channel_length * n); t = t + n){for (i = 0; i < n; i++){*(channel_output_matrix + (t / n) + (i * channel_length)) =soft_quant(*(channel_output_vector + (t + i)));//printf("%d ",*(channel_output_matrix + (t / n) + (i * channel_length)));}}/*结束设置:利用网格遍历开始译码通道,在编码完成后结束*/for (t = 0; t < channel_length - m; t++){if (t <= m)//假设从零,所以只是计算路径的所有零状态step = (int)pow(2.0, m - t * 1);//如果不写成2.0,会出现函数重载不明确的错误elsestep = 1;//利用state_history矩阵作为循环缓冲区sh_ptr = (int)((t + 1) % (depth_of_trellis + 1));//sh_ptr为state history矩阵的指针for (j = 0; j < number_of_states; j += step){//重复每个可能的卷积编码器的输出组for (l = 0; l < n; l++){branch_metric = 0;//计算每个通道符号的分支度量,以及所有的信道总和在卷积编码器的输出组信道符号binary_output[0] = (output[j][l] & 0x00000002) >> 1;binary_output[1] = output[j][l] & 0x00000001;branch_metric = branch_metric +abs(*(channel_output_matrix + (0 * channel_length + t)) - 7 *binary_output[0])+ abs(*(channel_output_matrix + (1 * channel_length + t)) - 7 * binary_output[1]);//选择累加误差最小的if (accum_err_metric[nextstate[j][l]][1] > accum_err_metric[j][0] + branch_metric){accum_err_metric[nextstate[j][l]][1] = accum_err_metric[j][0] + branch_metric;state_history[nextstate[j][l]][sh_ptr] = j;//printf("state_history[%d][%d]=%d\n", nextstate[j][l], sh_ptr,state_history[nextstate[j][l]][sh_ptr]);}} //循环l结束} //j结束,更新网格//accum_err_metric矩阵第二列移到第一列,第二列标志为一个很大的数for (j = 0; j < number_of_states; j++){accum_err_metric[j][0] = accum_err_metric[j][1];//printf("accum_err_metric[%d][0]=%d\n", j, accum_err_metric[j][0]);accum_err_metric[j][1] = MAXINT;}//如果网格填充完成,现在需要追踪{for (j = 0; j <= depth_of_trellis; j++)//初始化状态序列矩阵state_sequence[j] = 0;// 找到的最小累积state_history元素x = MAXINT;for (j = 0; j < (number_of_states / 2); j++){if(accum_err_metric[j][0] < accum_err_metric[number_of_states - 1 - j][0]){xx = accum_err_metric[j][0];hh = j;}else{xx = accum_err_metric[number_of_states - 1 - j][0];hh = number_of_states - 1 - j;}if (xx < x){x = xx;h = hh;}}state_sequence[depth_of_trellis] = h;for (j = depth_of_trellis; j > 0; j--){sh_col = j + (sh_ptr - depth_of_trellis);if (sh_col < 0)sh_col = sh_col + depth_of_trellis + 1;state_sequence[j - 1] = state_history[state_sequence[j]][sh_col];}//找出输入序列对应的状态序列在最佳路径*(decoder_output_matrix + t - depth_of_trellis + 1) =input[state_sequence[0]][state_sequence[1]];//printf("译码输出:%d\n", *(decoder_output_matrix + t - depth_of_trellis + 1));} //if状态} // 结束t循环//译码信道中的数据for (t = channel_length - m; t < channel_length; t++){sh_ptr = (int)((t + 1) % (depth_of_trellis + 1));last_stop = number_of_states / pow(2.0, t - channel_length + m);//不需要考虑输入的状态是1,所以确定最高可能的状态数是0for (j = 0; j < last_stop; j++){branch_metric = 0;deci2bin(output[j][0], n, binary_output);for (ll = 0; ll < n; ll++){branch_metric = branch_metric + soft_metric(*(channel_output_matrix + (ll * channel_length + t)), binary_output[ll]);}if ((accum_err_metric[nextstate[j][0]][1] > accum_err_metric[j][0] +branch_metric)){accum_err_metric[nextstate[j][0]][1] = accum_err_metric[j][0] +branch_metric;state_history[nextstate[j][0]][sh_ptr] = j;}}for (j = 0; j < number_of_states; j++){accum_err_metric[j][0] = accum_err_metric[j][1];accum_err_metric[j][1] = MAXINT;}//对所选路径进行选择{for (j = 0; j <= depth_of_trellis; j++)state_sequence[j] = 0;x = accum_err_metric[0][0];h = 0;for (j = 1; j < last_stop; j++){if (accum_err_metric[j][0] < x){x = accum_err_metric[j][0];h = j;}}state_sequence[depth_of_trellis] = h;for (j = depth_of_trellis; j > 0; j--){sh_col = j + (sh_ptr - depth_of_trellis);if (sh_col < 0)sh_col = sh_col + depth_of_trellis + 1;state_sequence[j - 1] = state_history[state_sequence[j]][sh_col];}*(decoder_output_matrix + t - depth_of_trellis + 1) =input[state_sequence[0]][state_sequence[1]];} //if条件状态} //结束t循环for (i = 1; i < depth_of_trellis - m; i++)*(decoder_output_matrix + channel_length - depth_of_trellis + i) =input[state_sequence[i]][state_sequence[i + 1]];free(channel_output_matrix);return;}//初始化三位软判决量化编码器//加入噪声后的量化void init_adaptive_quant(float es_ovr_n0){int i, d;float es, sn_ratio, sigma;es = 1;sn_ratio = (float)pow(10.0, (es_ovr_n0 / 10.0));sigma = (float)sqrt(es / (2.0 * sn_ratio));d = (int)(32 * 0.5 * sigma);for (i = -128; i < (-3 * d); i++)quantizer_table[i + 128] = 7;for (i = (-3 * d); i < (-2 * d); i++)quantizer_table[i + 128] = 6;for (i = (-2 * d); i < (-1 * d); i++)quantizer_table[i + 128] = 5;for (i = (-1 * d); i < 0; i++)quantizer_table[i + 128] = 4;for (i = 0; i < (1 * d); i++)quantizer_table[i + 128] = 3;for (i = (1 * d); i < (2 * d); i++)quantizer_table[i + 128] = 2;for (i = (2 * d); i < (3 * d); i++)quantizer_table[i + 128] = 1;for (i = (3 * d); i < 128; i++)quantizer_table[i + 128] = 0;}//channel_symbol信道中的值为-1或+1的加噪信号int soft_quant(float channel_symbol){int x;x = (int)(32.0 * channel_symbol); //则x的平均值为-32或+32if (x < -128) x = -128; //小于-128,输出128if (x > 127) x = 127; //大于127则输出127return(quantizer_table[x + 128]); //查找量化表}/* this metric is based on the algorithm given in Michelson and Levesque, page 323. */int soft_metric(int data, int guess){return(abs(data - (guess * 7)));//当给出当前状态以及输入时,计算下一个状态,并且计算卷积码编码内容int nxt_stat(int current_state, int input, int *memory_contents){int binary_state[K - 1]; //二进制当前状态int next_state_binary[K - 1]; //二进制下一个状态int next_state; //十进制当前状态int i;deci2bin(current_state, K - 1, binary_state);next_state_binary[0] = input;for (i = 1; i < K - 1; i++)next_state_binary[i] = binary_state[i - 1];next_state = bin2deci(next_state_binary, K - 1);memory_contents[0] = input;for (i = 1; i < K; i++)memory_contents[i] = binary_state[i - 1]; //编码内容return(next_state);//返回十进制的下一个状态}//将十进制数据转化为特殊的二进制,例如:十进制“100011”输出二进制数据100011 void deci2bin(int d, int size, int *b) {int i;for (i = 0; i < size; i++)b[i] = 0;b[size - 1] = d & 0x01;for (i = size - 2; i >= 0; i--) {d = d >> 1;b[i] = d & 0x01;}}//将2进制数据转化为特殊的十进制,例如:二进制“100011”输出十进制数据100011 int bin2deci(int *b, int size) {int i, d;d = 0;for (i = 0; i < size; i++)d += b[i] << (size - i - 1);return(d);}程序辅助模块// viterbi.cpp : 定义控制台应用程序的入口点。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

译码主要部分#include"stdafx.h"//#define DEBUGvoid deci2bin(int d, int size, int *b);int bin2deci(int *b, int size);int nxt_stat(int current_state, int input, int *memory_contents);void init_quantizer(void);void init_adaptive_quant(float es_ovr_n0);int soft_quant(float channel_symbol);int soft_metric(int data, int guess);int quantizer_table[256];void sdvd(int g[2][K], float es_ovr_n0, long channel_length, float*channel_output_vector, int *decoder_output_matrix){int i, j, l, ll; //循环控制变量long t; //时间int memory_contents[K]; //记录输入内容int input[TWOTOTHEM][TWOTOTHEM]; //对当前状态以及下一个状态映射int output[TWOTOTHEM][2]; //卷积码编码输出矩阵int nextstate[TWOTOTHEM][2]; //下一个状态矩阵int accum_err_metric[TWOTOTHEM][2]; //误差累计矩阵int state_history[TWOTOTHEM][K * 5 + 1]; //历史状态表int state_sequence[K * 5 + 1]; //状态序列int *channel_output_matrix; //信道输出序列int binary_output[2];int branch_output[2]; //0或者1的输出分支int m, n, number_of_states, depth_of_trellis, step, branch_metric,sh_ptr, sh_col, x, xx, h, hh, next_state, last_stop;n = 2; //1/2为卷积码传输数据的码率m = K - 1;//寄存器个数number_of_states = (int)pow(2.0, m);//状态个数number of states = 2^(K - 1) = 2^mdepth_of_trellis = K * 5;for (i = 0; i < number_of_states; i++){for (j = 0; j < number_of_states; j++)input[i][j] = 0; //输入数组初始化for (j = 0; j < n; j++){nextstate[i][j] = 0;//下一个状态数组初始化output[i][j] = 0; //输出数组初始化}for (j = 0; j <= depth_of_trellis; j++){state_history[i][j] = 0;//历史状态数组初始化state_history[4][16] }accum_err_metric[i][0] = 0;//误差累计矩阵第一列初始化为0accum_err_metric[i][1] = MAXINT;//误差累计矩阵第二列初始化为一个很大的数}/*前向纠错简称FEC(Forward Error Correction),其原理是:发送方将要发送的数据附加上一定的冗余纠错码一并发送,接收方则根据纠错码对数据进行差错检测,如发现差错,由接收方进行纠正*//*产生状态转移矩阵、输出矩阵、输入矩阵*///输入矩阵表示的是FEC编码传输给下一个状态//下一个状态由输入和当前状态给出//输出矩阵for (j = 0; j < number_of_states; j++){for (l = 0; l < n; l++){next_state = nxt_stat(j, l, memory_contents);input[j][next_state] = l;/*计算给定的卷积编码器输出当前状态数和输入值*/branch_output[0] = 0;branch_output[1] = 0;for (i = 0; i < K; i++){branch_output[0] = branch_output[0] ^ memory_contents[i] & g[0][i];branch_output[1] = branch_output[1] ^ memory_contents[i] & g[1][i];}nextstate[j][l] = next_state;//下一个状态output[j][l] = bin2deci(branch_output, 2);//输出十进制}}#ifdef DEBUGprintf("\nInput:");for (j = 0; j < number_of_states; j++){printf("\n");for (l = 0; l < number_of_states; l++)printf("%2d ", input[j][l]);}printf("\nOutput:");for (j = 0; j < number_of_states; j++){printf("\n");for (l = 0; l < n; l++)printf("%2d ", output[j][l]);}printf("\nNext State:");for (j = 0; j < number_of_states; j++){printf("\n");for (l = 0; l < n; l++)printf("%2d ", nextstate[j][l]);}#endifchannel_output_matrix =(int *)malloc(channel_length * sizeof(int));if (channel_output_matrix == NULL){printf("allocation is failure!!\n");exit(1);}printf("\n");/*信道输出为n行,2列,每行对应于一个通道符号给定的位和每一列对应于一个已编码的位*/ channel_length = channel_length / n;init_adaptive_quant(es_ovr_n0);//进行优化,匹配过信噪比的量化矩阵//量化信道输出,将浮点型数据转化为整形for (t = 0; t < (channel_length * n); t = t + n){for (i = 0; i < n; i++){*(channel_output_matrix + (t / n) + (i * channel_length)) =soft_quant(*(channel_output_vector + (t + i)));//printf("%d ",*(channel_output_matrix + (t / n) + (i * channel_length)));}}/*结束设置:利用网格遍历开始译码通道,在编码完成后结束*/for (t = 0; t < channel_length - m; t++){if (t <= m)//假设从零,所以只是计算路径的所有零状态step = (int)pow(2.0, m - t * 1);//如果不写成2.0,会出现函数重载不明确的错误elsestep = 1;//利用state_history矩阵作为循环缓冲区sh_ptr = (int)((t + 1) % (depth_of_trellis + 1));//sh_ptr为state history矩阵的指针for (j = 0; j < number_of_states; j += step){//重复每个可能的卷积编码器的输出组for (l = 0; l < n; l++){branch_metric = 0;//计算每个通道符号的分支度量,以及所有的信道总和在卷积编码器的输出组信道符号binary_output[0] = (output[j][l] & 0x00000002) >> 1;binary_output[1] = output[j][l] & 0x00000001;branch_metric = branch_metric +abs(*(channel_output_matrix + (0 * channel_length + t)) - 7 *binary_output[0])+ abs(*(channel_output_matrix + (1 * channel_length + t)) - 7 * binary_output[1]);//选择累加误差最小的if (accum_err_metric[nextstate[j][l]][1] > accum_err_metric[j][0] + branch_metric){accum_err_metric[nextstate[j][l]][1] = accum_err_metric[j][0] + branch_metric;state_history[nextstate[j][l]][sh_ptr] = j;//printf("state_history[%d][%d]=%d\n", nextstate[j][l], sh_ptr,state_history[nextstate[j][l]][sh_ptr]);}} //循环l结束} //j结束,更新网格//accum_err_metric矩阵第二列移到第一列,第二列标志为一个很大的数for (j = 0; j < number_of_states; j++){accum_err_metric[j][0] = accum_err_metric[j][1];//printf("accum_err_metric[%d][0]=%d\n", j, accum_err_metric[j][0]);accum_err_metric[j][1] = MAXINT;}//如果网格填充完成,现在需要追踪{for (j = 0; j <= depth_of_trellis; j++)//初始化状态序列矩阵state_sequence[j] = 0;// 找到的最小累积state_history元素x = MAXINT;for (j = 0; j < (number_of_states / 2); j++){if(accum_err_metric[j][0] < accum_err_metric[number_of_states - 1 - j][0]){xx = accum_err_metric[j][0];hh = j;}else{xx = accum_err_metric[number_of_states - 1 - j][0];hh = number_of_states - 1 - j;}if (xx < x){x = xx;h = hh;}}state_sequence[depth_of_trellis] = h;for (j = depth_of_trellis; j > 0; j--){sh_col = j + (sh_ptr - depth_of_trellis);if (sh_col < 0)sh_col = sh_col + depth_of_trellis + 1;state_sequence[j - 1] = state_history[state_sequence[j]][sh_col];}//找出输入序列对应的状态序列在最佳路径*(decoder_output_matrix + t - depth_of_trellis + 1) =input[state_sequence[0]][state_sequence[1]];//printf("译码输出:%d\n", *(decoder_output_matrix + t - depth_of_trellis + 1));} //if状态} // 结束t循环//译码信道中的数据for (t = channel_length - m; t < channel_length; t++){sh_ptr = (int)((t + 1) % (depth_of_trellis + 1));last_stop = number_of_states / pow(2.0, t - channel_length + m);//不需要考虑输入的状态是1,所以确定最高可能的状态数是0for (j = 0; j < last_stop; j++){branch_metric = 0;deci2bin(output[j][0], n, binary_output);for (ll = 0; ll < n; ll++){branch_metric = branch_metric + soft_metric(*(channel_output_matrix + (ll * channel_length + t)), binary_output[ll]);}if ((accum_err_metric[nextstate[j][0]][1] > accum_err_metric[j][0] +branch_metric)){accum_err_metric[nextstate[j][0]][1] = accum_err_metric[j][0] +branch_metric;state_history[nextstate[j][0]][sh_ptr] = j;}}for (j = 0; j < number_of_states; j++){accum_err_metric[j][0] = accum_err_metric[j][1];accum_err_metric[j][1] = MAXINT;}//对所选路径进行选择{for (j = 0; j <= depth_of_trellis; j++)state_sequence[j] = 0;x = accum_err_metric[0][0];h = 0;for (j = 1; j < last_stop; j++){if (accum_err_metric[j][0] < x){x = accum_err_metric[j][0];h = j;}}state_sequence[depth_of_trellis] = h;for (j = depth_of_trellis; j > 0; j--){sh_col = j + (sh_ptr - depth_of_trellis);if (sh_col < 0)sh_col = sh_col + depth_of_trellis + 1;state_sequence[j - 1] = state_history[state_sequence[j]][sh_col];}*(decoder_output_matrix + t - depth_of_trellis + 1) =input[state_sequence[0]][state_sequence[1]];} //if条件状态} //结束t循环for (i = 1; i < depth_of_trellis - m; i++)*(decoder_output_matrix + channel_length - depth_of_trellis + i) =input[state_sequence[i]][state_sequence[i + 1]];free(channel_output_matrix);return;}//初始化三位软判决量化编码器//加入噪声后的量化void init_adaptive_quant(float es_ovr_n0){int i, d;float es, sn_ratio, sigma;es = 1;sn_ratio = (float)pow(10.0, (es_ovr_n0 / 10.0));sigma = (float)sqrt(es / (2.0 * sn_ratio));d = (int)(32 * 0.5 * sigma);for (i = -128; i < (-3 * d); i++)quantizer_table[i + 128] = 7;for (i = (-3 * d); i < (-2 * d); i++)quantizer_table[i + 128] = 6;for (i = (-2 * d); i < (-1 * d); i++)quantizer_table[i + 128] = 5;for (i = (-1 * d); i < 0; i++)quantizer_table[i + 128] = 4;for (i = 0; i < (1 * d); i++)quantizer_table[i + 128] = 3;for (i = (1 * d); i < (2 * d); i++)quantizer_table[i + 128] = 2;for (i = (2 * d); i < (3 * d); i++)quantizer_table[i + 128] = 1;for (i = (3 * d); i < 128; i++)quantizer_table[i + 128] = 0;}//channel_symbol信道中的值为-1或+1的加噪信号int soft_quant(float channel_symbol){int x;x = (int)(32.0 * channel_symbol); //则x的平均值为-32或+32if (x < -128) x = -128; //小于-128,输出128if (x > 127) x = 127; //大于127则输出127return(quantizer_table[x + 128]); //查找量化表}/* this metric is based on the algorithm given in Michelson and Levesque, page 323. */int soft_metric(int data, int guess){return(abs(data - (guess * 7)));//当给出当前状态以及输入时,计算下一个状态,并且计算卷积码编码内容int nxt_stat(int current_state, int input, int *memory_contents){int binary_state[K - 1]; //二进制当前状态int next_state_binary[K - 1]; //二进制下一个状态int next_state; //十进制当前状态int i;deci2bin(current_state, K - 1, binary_state);next_state_binary[0] = input;for (i = 1; i < K - 1; i++)next_state_binary[i] = binary_state[i - 1];next_state = bin2deci(next_state_binary, K - 1);memory_contents[0] = input;for (i = 1; i < K; i++)memory_contents[i] = binary_state[i - 1]; //编码内容return(next_state);//返回十进制的下一个状态}//将十进制数据转化为特殊的二进制,例如:十进制“100011”输出二进制数据100011 void deci2bin(int d, int size, int *b) {int i;for (i = 0; i < size; i++)b[i] = 0;b[size - 1] = d & 0x01;for (i = size - 2; i >= 0; i--) {d = d >> 1;b[i] = d & 0x01;}}//将2进制数据转化为特殊的十进制,例如:二进制“100011”输出十进制数据100011 int bin2deci(int *b, int size) {int i, d;d = 0;for (i = 0; i < size; i++)d += b[i] << (size - i - 1);return(d);}程序辅助模块// viterbi.cpp : 定义控制台应用程序的入口点。

相关文档
最新文档