卷积码编码和维特比译码
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;}。
卷积码及维特比译码 notes

卷积码在CDMA中使用,表示方式包含多项式表示法,状态转移图法和网格图法。
编码器:n (编码后输出的码字长度)(模二加法器个数)k (输入的比特信息长度)(每个寄存器的位数)L (约束长度/记忆深度)(寄存器个数)R表示为R = k/n。
莫二加法器以寄存器的数位为单位进行选取计算,计算方式自定义,形成网格图作为码本莫二加法运算,等同于“异或”运算。
两个序列按位相加,即两个序列中对应位相加,不进位。
效果是相同为0,不同为1。
例: 1+1 = 0+0 = 01+0 = 0+1 = 10 1 0 1+ 0 0 1 1──────0 1 1 0下图是输出码长n=2,输入比特k=1,记忆深度L=2的,(2,1,2)卷积码编码器。
如编码序列“0 1 1 0 0”在图中的序列如下:汉明距离两个二进制数之间进行逐位对比,得到不同的个数如1000,与1100为1,与1110为2,与1111为3维特比算法综合状态之间的转移概率和前一层各状态的概率情况计算出概率最大的状态转换路径,从而推断出隐含状态的序列的情况。
的分支度量(汉明距离)。
其中有两条路径的分支量度为0。
3.寻找最大似然路径 - 译码过程维特比算法的关键点在于,接收机可以使用分支度量和先前计算的路径度量递推地计算当前状态的路径度量。
初始时,状态00代价为0,其它状态代价为正无穷(∞)。
算法的主循环由两个主要步骤组成:首先计算下一时刻监督比特序列的分支度量,然后计算该时刻各状态的路径度量。
路径度量的计算可以认为是一个“加-比-选”的过程1)将分支度量与上一时刻状态的路径度量相加。
2)每一状态比较来自前一时刻状态可达到的所有路径(只有两条这样的路径进行比较)3)每一状态删除其余到达路径,选择最小度量的路径保留(称为幸存路径/存活路径)若进入某个状态的部分路径中,有两条路径的度量值相等,则可以任选其一作为幸存路径。
下图显示了维特比译码的过程。
此例接收到的位序列为11 10 11 00 01 10(偷偷告诉你:这是有误码的信息)此时,产生了具有相同路径度量的四个不同路径,通向这四个状态的任一路径都是可能发送的比特序列(它们都具有度量为2的汉明距离)。
卷积码的维特比译码

卷积码编码原理 维特比译码原理 Matlab实验 DSP仿真实验
1.卷积码编码
卷积编码则把k比特信息段编成n比特的 码组,但所编的n长码组不仅同当前的k 比特信息段有关联,而且还同前面的(N1)个信息段有关联,人们常称这N为该卷 积码的约束长度。
一般来说,对于卷积码,k和n是较小的 整数, 常把卷积码记作(n,k,N)卷积码,它 的编码效率为R=k/n。
据
主程序完成 对输入比特的
译码
探针将已译 码的数据写入
文件
译码实验结果对比
MATLAB译码 结果
DSP译码结 果
图2 :(2,1,2)卷积码的树状图
2.
格图也称网络图或篱笆图,它由状态图在时间上展开而得到,如图3所示。图 中画出了所有可能的数据输入时, 状态转移的全部可能轨迹,实线表示数据为 0, 虚线表示数据为 1, 线旁数字为输出码字,节点表示状态。
图3:(2,1,2)卷积码的网格图
2.卷积码的维特比译码
4.卷积码编码的C54x Simulator仿 真试验
初始化
探针从PC文件 中读取比特数
据
主程序完成 对输入比特的
卷积编码
探针将已编 码的数据写入
文件
编码结果对比
未编码数据
DSP编码后的数 据
MATLAB编码 后的数据
Viterbi译码的C54x Simulator仿真 试验
初始化
探针从PC文件 中读取比特数
维特比译码是一种最大似然译码算法。 最大似然译码算法的基本思路是: 把接收码字与所有可能的码字比较,选择一种码距最小的码字作为解码输出。 由于接收序列通常很长,所以维特比译码时最大似然译码做了简化, 即它把接 收码字分段累接处理,每接收一段码字,计算、 比较一次, 保留码距最小的路 径,直至译完整个序列。
卷积码编码和维特比译码的原理、性能与仿真分析

卷积码编码和维特比译码的原理、性能与仿真分析1.引言卷积码的编码器是由一个有k位输入、n位输出,且具有m位移位寄存器构成的有限状态的有记忆系统,通常称它为时序网络。
编码器的整体约束长度为v,是所有k个移位寄存器的长度之和。
具有这样的编码器的卷积码称作[n,k,v]卷积码。
对于一个(n,1,v)编码器,约束长度v等于存储级数m.卷积码是由k个信息比特编码成n(n>k)比特的码组,编码出的n比特码组值不仅与当前码字中的k个信息比特值有关,而且与其前面v个码组中的v*k个信息比特值有关。
卷积码有三种译码方式:序列译码、门限译码和概率译码。
其中,概率译码根据最大似然译码原理在所有可能路径中求取与接收路径最相似的一条路径,具有最佳的纠错性能,维特比译码是概率译码中极重要的一种方式。
序列译码和门限译码则不一定能找出与接收路径最相似的一条路径。
不同于维特比译码,门限译码与序列译码所需的计算量是可变的且对于给定信息分组的最终判决仅仅基于(m+1)个接收分组,而不是基于整个接收序列。
与维特比译码所使用的对数似然量度不同,序列译码所使用的量度为Fano量度。
在接收序列受扰严重的情况下,序列译码的计算量大于维特比译码所需的固定计算量,虽然序列译码要求的平均计算次数通常小于维特比译码。
在采用并行处理的情况下,维特比译码的速度会优于序列译码。
在同样码率和存储级数的条件下,门限译码的性能比维特比译码低大约3dB.维特比译码的数据输出方式有硬判决及软判决两种方式,本文选取生成多项式为561,753的(2,1,8)卷积码对硬判决的性能进行分析,并依据维特比译码的原理以及卷积码的特性,对卷积码编码和维特比译码过程在加性高斯白噪声(AWGN)信道下进行仿真,并且根据仿真结果对维特比译码(硬判决)的结果进行分析。
由于卷积码的生成可以看做一个马尔科夫过程,因此,不同状态间的转移概率对描述这个过程有极关键的作用。
本文则基于MATLAB对不同状态间的转移概率进行求解,从而更准确地分析维特比译码的性能。
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 译码算法和BCJR 译码算法。
基于某种准那么,这两种算法都是最优的。
1967 年,Viterbi 提出了卷积码的Viterbi 译码算法,后来Omura 证明Viterbi 译码算法等效于在加权图中寻找最优途径问题的一个动态规划〔Dynamic Programming〕解决方案,随后,Forney 证明它实际上是最大似然〔ML,Maximum Likelihood〕译码算法,即译码器选择输出的码字通常使接收序列的条件概率最大化。
BCJR 算法是1974 年提出的,它实际上是最大后验概率〔MAP,Maximum A Posteriori probability〕译码算法。
这两种算法的最优化目的略有不同:在MAP 译码算法中,信息比特错误概率是最小的,而在ML 译码算法中,码字错误概率是最小的,但两种译码算法的性能在本质上是一样的。
由于Viterbi 算法实现更简单,因此在实际应用比较广泛,但在迭代译码应用中,例如逼近Shannon 限的Turbo 码,常使用BCJR 算法。
另外,在迭代译码应用中,还有一种Viterbi 算法的变种:软输出Viterbi 算法〔SOV A,Soft-Output Viterbi Algorithm〕,它是Hagenauer 和Hoeher 在1989 年提出的。
为了理解Viterbi 译码算法,我们需要将编码器状态图按时间展开〔因为状态图不能反映出时间变化情况〕,即在每个时间单元用一个分隔开的状态图来表示。
例如〔3,1,2〕非系统前馈编码器,其生成矩阵为:G(D)=[1+D1+D21+D+D2]〔1〕图1 〔a〕〔3,1,2〕编码器〔b〕网格图〔h=5〕假定信息序列长度为h=5,那么网格图包含有h+m+1=8 个时间单元,用0 到h+m=7 来标识,如图1〔b〕所示。
假设编码器总是从全0 态S0 开始,又回到全0 态,前m=2 个时间单元对应于编码器开始从S0“启程〞,最后m=2 个时间单元对应于向S0“返航〞。
动态规划:卷积码Viterbi译码算法

动态规划:卷积码的Viterbi译码算法学院:网研院姓名:xxx 学号:xxx 一、动态规划原理动态规划(dynamic programming)是运筹学的一个分支,是求解决策过程(decision process)最优化的数学方法。
动态规划算法通常用于求解具有某种最优性质的问题。
在这类问题中,可能会有许多可行解,每一个解都对应于一个值,我们希望找到具有最优值的解。
动态规划算法与分治法类似,其基本思想也是将待求解问题分解成若干个子问题,先求解子问题,然后从这些子问题的解得到原问题的解。
与分治法不同的是,适合于用动态规划求解的问题,经分解得到子问题往往不是互相独立的。
若用分治法来解这类问题,则分解得到的子问题数目太多,有些子问题被重复计算了很多次。
如果我们能够保存已解决的子问题的答案,而在需要时再找出已求得的答案,这样就可以避免大量的重复计算,节省时间。
动态规划程序设计是对解最优化问题的一种途径、一种方法,而不是一种特殊算法。
不象搜索或数值计算那样,具有一个标准的数学表达式和明确清晰的解题方法。
动态规划程序设计往往是针对一种最优化问题,由于各种问题的性质不同,确定最优解的条件也互不相同,因而动态规划的设计方法对不同的问题,有各具特色的解题方法,而不存在一种万能的动态规划算法,可以解决各类最优化问题。
二、卷积码的Viterbi译码算法简介在介绍维特比译码算法之前,首先了解一下卷积码编码,它常常与维特比译码结合使用。
(2,1,3)卷积码编码器是最常见的卷积码编码器,在本次实验中也使用了(2,1,3)卷积码编码器,下面介绍它的原理。
(2,1,3)卷积码是把信源输出的信息序列,以1个码元为一段,通过编码器输出长为2的一段码段。
该码段的值不仅与当前输入码元有关,而且也与其之前的2个输入码元有关。
如下图所示,输出out1是输入、第一个编码器存储的值和第二个编码器存储的值逻辑加操作的结果,输出out2是输入和第二个编码器存储的值逻辑加操作的结果。
卷积码编码和维特比译码

卷积码编码维特比译码实验设计报告SUN一、实验目的掌握卷积码编码和维特比译码的基本原理,利用了卷积码的特性, 运用网格图和回溯以得到译码输出。
二、实验原理1.卷积码是由连续输入的信息序列得到连续输出的已编码序列。
其编码器将k个信息码元编为n个码元时,这n个码元不仅与当前段的k个信息有关,而且与前面的(m-1)段信息有关(m为编码的约束长度)。
2.一般地,最小距离d表明了卷积码在连续m段以内的距离特性,该码可以在m个连续码流内纠正(d-1)/2个错误。
卷积码的纠错能力不仅与约束长度有关,还与采用的译码方式有关。
3. 维特比译码算法基本原理是将接收到的信号序列和所有可能的发送信号序列比较,选择其中汉明距离最小的序列认为是当前发送序列。
卷积码的Viterbi 译码是根据接收码字序列寻找编码时通过网格图最佳路径的过程,找到最佳路径即完成了译码过程,并可以纠正接收码字中的错误比特。
4.所谓“最佳”, 是指最大后验条件概率:P( C/ R) = max [ P ( Cj/ R) ] , 一般来说, 信道模型并不使用后验条件概率,因此利用Beyes 公式、根据信道特性出结论:max[ P ( Cj/ R) ]与max[ P ( R/ Cj) ]等价。
考虑到在系统实现中往往采用对数形式的运算,以求降低运算量,并且为求运算值为整数加入了修正因子a1 、a2 。
令M ( R/ Cj) = log[ P ( R/ Cj) ] =Σa1 (log[ P( Rm/ Cmj ) ] + a2) 。
其中, M 是组成序列的码字的个数。
因此寻找最佳路径, 就变成寻找最大M( R/ Cj) , M( R/ Cj) 称为Cj 的分支路径量度,含义为发送Cj 而接收码元为R的似然度。
5.卷积码的viterbi译码是根据接收码字序列寻找编码时通过网格图最佳路径的过程,找到最佳路径即完成了译码过程并可以纠正接收码字中的错误比特。
三、实验代码#include<stdio.h>#include "Conio.h"#define N 7#include "math.h"#include <stdlib.h>#include<time.h>#define randomize() srand((unsigned)time(NULL))encode(unsigned int *symbols, /*编码输出*/unsigned int *data, /*编码输入*/unsigned int nbytes, /*nbytes=n/16,n为实际输入码字的数目*/unsigned int startstate /*定义初始化状态*/)////////////////////////////////////////////////////////////////////////////卷积码编码///////////////////////////////////////////////////////////////////////////////{unsigned int j;unsigned int input,a1=0,a2=0,a3=0,a4=0,a5=0,a6=0;for(j=0;j<nbytes;j++){input=*data;data++;*symbols = input^a1^a2^a3^a6; //c1(171)symbols++;*symbols = input^a2^a3^a5^a6; //c2(133)symbols++;a2=a1;a1=input;}return 0;}int trandistance(int m, int state1, int state2)/*符号m与从state1到state2时输出符号的汉明距离,如果state1无法到state2则输出度量值为100*/{int c;int sym,sym1,sym2;sym1=((state2>>1)&1)^(state2&1)^(state1&1);sym2=((state2>>1)&1)^(state1&1);sym=(sym1<<1) | sym2;if ( ((state1&2)>>1)==(state2&1))c=((m&1)^(sym&1))+(((m>> 1)&1)^((sym >> 1)&1));elsec=10000;return(c);}int traninput(int a,int b) /*状态从a到b时输入卷积码的符号*/{int c;c=((b&2)>>1);return(c);}int tranoutput(int a,int b) /*状态从a到b时卷积码输出的符号*/{int c,s1,s2;s1=(a&1)^((a&2)>>1)^((b&2)>>1);s2=(a&1)^((b&2)>>1);c=(s1<<1)|s2;return(c);}////////////////////////////////////////////////////////////////////////////维特比译码///////////////////////////////////////////////////////////////////////////////void viterbi(int initialstate, /*定义解码器初始状态*/int *viterbiinput, /*解码器输入码字序列*/int *viterbioutput /*解码器输出码字序列*/){struct sta /*定义网格图中每一点为一个结构体,其元素包括*/ {int met; /*转移到此状态累计的度量值*/int value; /*输入符号*/struct sta *last; /*及指向前一个状态的指针*/};struct sta state[4][N];struct sta *g,*head;int i,j,p,q,t,r,u,l;for(i=0;i<4;i++) /* 初始化每个状态的度量值*/for(j=0;j<N;j++)state[i][j].met=0;for(l=0;l<4;l++){state[l][0].met=trandistance(*viterbiinput,initialstate,l);state[l][0].value=traninput(initialstate,l);state[l][0].last=NULL;}viterbiinput++; /*扩展第一步幸存路径*/for(t=1;t<N;t++){for(p=0;p<4;p++){state[p][t].met=state[0][t-1].met+trandistance(*viterbiinput,0,p);state[p][t].value=traninput(0,p);state[p][t].last=&state[0][t-1];for(q=0;q<4;q++){if(state[q][t-1].met+trandistance(*viterbiinput,q,p)<state[p][t].met){state[p][t].met=state[q][t-1].met+trandistance(*viterbiinput,q,p);state[p][t].value=traninput(q,p);state[p][t].last=&state[q][t-1];}}}viterbiinput++;} /*计算出剩余的幸存路径*/r=state[0][N-1].met; /*找出n步后度量值最小的状态准备回溯路由*/g=&state[0][N-1];for(u=N;u>0;u--) /*向前递归的找出最大似然路径*/{*(viterbioutput+(u-1))=g->value;g=g->last;}/* for(u=0;u<8;u++)*(viterbioutput+u)=state[u][2].met; */ /*此行程序可用于检测第n列的度量值*/}void decode(unsigned int *input, int *output,int n){int viterbiinput[100];int j;for(j=0;j<n+2;j++){viterbiinput[j]=(input[j*2]<<1)|input[j*2+1];}viterbi(0,viterbiinput,output);}void main(){unsigned intencodeinput[100],wrong[10]={0,0,0,0,0,0,0,0,0,0},encodeoutput[100];int n=5,i,m,j=0,decodeinput[100],decodeoutput[100];randomize();for(i=0; i<n; i++)encodeinput[i]=rand()%2;encodeinput[n]= encodeinput[n+1]=0;encode(encodeoutput,encodeinput,n+2,0);printf("the input of encoder is :\n"); //信息源输入的信息码(随机产生)for(i=0;i<n; i++)printf("%2d",encodeinput[i]);printf("\n");printf("the output of encoder is :\n"); //编码之后产生的卷积码for(i=0;i<(n+2)*2;i++){printf("%2d",encodeoutput[i]);if(i%20==19)printf("\n");}printf("\n");printf("please input the number of the wrong bit\n"); //信道传输收到干扰而产生的错误码scanf("%d",&m);printf("please input the positions of the wrong bit(0-9)\n");for(i=0;i<m;i++){scanf("%d",&wrong[m]);if(encodeoutput[wrong[m]]==0)encodeoutput[wrong[m]]=1;elseencodeoutput[wrong[m]]=0;}printf("the input of decoder is :\n");for(i=0;i<(n+2)*2;i++){printf("%2d",encodeoutput[i]);if(i%20==19)printf("\n");}printf("\n");decode(encodeoutput,decodeoutput,n+2);printf("the output of decoder is :\n");for(i=0;i<n;i++)printf("%2d",decodeoutput[i]);printf("\n");for(i=0;i<n;i++){if(encodeinput[i]!=decodeoutput[i])j++;}printf("the number of incorrect bit is:%d\n",j);}四、实验总结(1)了解实验原理,分析实验所占数组变量很重要,也是相对考虑较多的;(2)对于读写文件,通过本实验更加熟悉;(3)记录实验程序最佳路径是本实验的难点;。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
{
for(p=0;p<4;p++)
{
state[p][t].met=state[0][t-1].met+trandistance(*viterbiinput,0,p);
state[p][t].value=traninput(0,p);
state[p][t].last=&state[0][t-1];
度量值*/
}
void decode(unsigned int *input, int *output,int n)
{
int viterbiinput[100];
int j;
for(j=0;j<n+2;j++)
{
viterbiinput[j]=(input[j*2]<<1)|input[j*2Ʋl++)
{
state[l][0].met=trandistance(*viterbiinput,initialstate,l);
state[l][0].value=traninput(initialstate,l);
state[l][0].last=NULL;
}
viterbiinput++; /*扩展第一步幸存路径*/
state[p][t].value=traninput(q,p);
state[p][t].last=&state[q][t-1];
}
}
}
viterbiinput++;
} /*计算出剩余的幸存路径*/
r=state[0][N-1].met; /*找出n步后度量值最小的状态准备回溯路由*/
g=&state[0][N-1];
randomize();
for(i=0; i<n; i++)
encodeinput[i]=rand()%2;
encodeinput[n]= encodeinput[n+1]=0;
encode(encodeoutput,encodeinput,n+2,0);
printf("the input of encoder is :\n"); //信息源输入的信息码(随机产生)
卷积码编码维特比译码实验设计报告
SUN
一、实验目的
掌握卷积码编码和维特比译码的基本原理,利用了卷积码的特性,运用网格图和回溯以得到译码输出。
二、实验原理
1.卷积码是由连续输入的信息序列得到连续输出的已编码序列。其编码器将k个信息码元编为n个码元时,这n个码元不仅与当前段的k个信息有关,而且与前面的 (m-1)段信息有关(m为编码的约束长度)。
for(q=0;q<4;q++)
{
if(state[q][t-1].met+trandistance(*viterbiinput,q,p)<state[p][t].met)
{state[p][t].met=state[q][t-1].met+trandistance(*viterbiinput,q,p);
2.一般地,最小距离d表明了卷积码在连续m段以内的距离特性,该码可以在m个连续码流内纠正(d-1)/2个错误。卷积码的纠错能力不仅与约束长度有关,还与采用的译码方式有关。
3.维特比译码算法基本原理是将接收到的信号序列和所有可能的发送信号序列比较,选择其中汉明距离最小的序列认为是当前发送序列。卷积码的Viterbi译码是根据接收码字序列寻找编码时通过网格图最佳路径的过程,找到最佳路径即完成了译码过程,并可以纠正接收码字中的错误比特。
for(i=0;i<n; i++)
printf("%2d",encodeinput[i]);
printf("\n");
printf("the output of encoder is :\n"); //编码之后产生的卷积码
for(i=0;i<(n+2)*2;i++)
5.卷积码的viterbi译码是根据接收码字序列寻找编码时通过网格图最佳路径的过程,找到最佳路径即完成了译码过程 并可以纠正接收码字中的错误比特。
三、实验代码
#include<>
#include ""
#define N 7
#include ""
#include <>
#include<>
#define randomize() srand((unsigned)time(NULL))
for(u=N;u>0;u--) /*向前递归的找出最大似然路径*/
{
*(viterbioutput+(u-1))=g->value;
g=g->last;
}
/* for(u=0;u<8;u++)
*(viterbioutput+u)=state[u][2].met; */ /*此行程序可用于检测第n列的
}
viterbi(0,viterbiinput,output);
}
void main()
{
unsigned int
encodeinput[100],wrong[10]={0,0,0,0,0,0,0,0,0,0},encodeoutput[100];
int n=5,i,m,j=0,decodeinput[100],decodeoutput[100];
4.所谓“最佳”,是指最大后验条件概率:P( C/ R) = max [ P ( Cj/ R) ] ,一般来说,信道模型并不使用后验条件概率,因此利用Beyes公式、根据信道特性出结论:max[ P ( Cj/ R) ]与max[ P ( R/ Cj) ]等价。考虑到在系统实现中往往采用对数形式的运算,以求降低运算量,并且为求运算值为整数加入了修正因子a1、a2。令M ( R/ Cj) = log[ P ( R/ Cj) ] =Σa1 (log[ P( Rm/ Cmj ) ] + a2)。其中, M是组成序列的码字的个数。因此寻找最佳路径,就变成寻找最大M( R/ Cj) , M( R/ Cj)称为Cj的分支路径量度,含义为发送Cj而接收码元为R的似然度。
encode(
unsigned int *symbols, /*编码输出*/
unsigned int *data, /*编码输入*/
unsigned int nbytes, /*nbytes=n/16,n为实际输入码字的数目*/
unsigned int startstate /*定义初始化状态*/
)