• Matlab实现 维特比译码
千次阅读
2021-06-12 15:43:12

# 利用matlab实现维特比译码

## 话不多说，直接展示代码

1. 首先是利用卷积码进行编码
2. 然后加入awgn噪声
3. 最后利用维特比译码
4. 绘制误比特率图像

同时我也上传所有的代码文件
Matlab 实现维特比译码

在使用时，首先运行encode，对于随机产生的N=1000个比特流进行卷积编码
然后直接运行show errorbit，它会调用维特比译码的函数viterbi_hard()，同时绘制出在加入awgn白噪声后，不同信噪比下的误比特率。

卷积码进行编码
生成向量为 g1: [1,1,1], g2: [1,0,1]

encode

%%  初始化参数
N = 1000; %序列长度

code_in = randi(2,1,N)-1;
%code_in = [1,1,0,1,1,1,0,0,1];
N =  length(code_in);
g1 = [1,1,1];
g2 =  [1,0,1];

%% 计算 卷积码
%x_g1

x0 = g1(1) * code_in;
x1 = g1(2) * code_in;
x2 = g1(3)* code_in;
x = xor( x0(2:N), x1(1:N-1) );
x = [ x0(1), x, x1(N) ];
x3 = xor(x(3:N+1),x2(1:N-1));
x_g1 = [x(1:2) ,x3, x2(N) ];
%x_g2
x0 = g2(1) * code_in;
x1 = g2(2) * code_in;
x2 = g2(3)* code_in;
x = xor( x0(2:N), x1(1:N-1) );
x = [ x0(1), x, x1(N) ];
x3 = xor(x(3:N+1),x2(1:N-1));
x_g2 = [x(1:2) ,x3, x2(N) ];

%% 合并两路

x = zeros(1,size(x_g1, 2)+size(x_g2, 2));
x(1:2:end) = x_g1;
x(2:2:end) = x_g2;
x = x(1:length(x) - 4);



维特比译码
将状态转移图存在graph中
decode

function m = viterbi_hard(x)
%%
% a 00
% b 10
% c 01
% d 11
%x = [ 1 1 0 1 0 1 0 0 0 1];
sizex = size(x);
s = sizex(2)/2;
x = [0,0,0,0, x];
% to record the value
val_a = 0;
val_b = 0;
val_c = 0;
val_d = 0;

% graph       aa0    ab1    bc0   bd1    ca0     cb1   dc0   dd1
gra =          [ 0,0;    1,1;    1,0;    0,1;    1,1;    0,0;    0,1;    1,0  ];

% to record route
ma = zeros(1,s+2);
mb = zeros(1,s+2);
mc = zeros(1,s+2);
md = zeros(1,s+2);
%% stage 1 and 2

val_a = val_a + dis(gra(1,:), x(1:2));
ma(1)=0;
val_b = val_b + dis(gra(2,:), x(1:2));
mb(1)=1;

mc = mb;
md = mb;
val_c = val_b + dis(gra(3,:), x(3:4));
mc(2)=0;
val_d = val_b + dis(gra(4,:), x(3:4));
md(2)=1;

val_a = val_a + dis(gra(1,:), x(3:4));
ma(2)=0;
val_b = val_b + dis(gra(2,:), x(3:4));
mb(2)=1;

for i = 3:s+2
%     val_a_t =val_a;
%     val_b_t =val_b;
%     val_c_t =val_c;
%     val_d_t =val_d;
%     tempa = ma;
%     tempb = mb;
%     tempc = mc;
%     tempd = md;
% for val_a
if val_a + dis(gra(1,:), x(2*i-1:2*i)) >= val_c + dis(gra(5,:),x(2*i-1:2*i))
tempa = mc;
val_a_t = val_c + dis(gra(5,:),x(2*i-1:2*i));
tempa(i)=0;
else
val_a_t = val_a + dis(gra(1,:),x(2*i-1:2*i));
tempa = ma;
tempa(i)=0;
end
%for val_b
if val_a + dis(gra(2,:), x(2*i-1:2*i)) >= val_c + dis(gra(6,:),x(2*i-1:2*i))
tempb = mc;
val_b_t = val_c + dis(gra(6,:),x(2*i-1:2*i));
tempb(i)=1;
else
val_b_t = val_a + dis(gra(2,:),x(2*i-1:2*i));
tempb = ma;
tempb(i)=1;
end

%for val_c
if val_b + dis(gra(3,:), x(2*i-1:2*i)) >= val_d + dis(gra(7,:),x(2*i-1:2*i))
tempc = md;
val_c_t = val_d + dis(gra(7,:),x(2*i-1:2*i));
tempc(i)=0;
else
val_c_t = val_b + dis(gra(3,:),x(2*i-1:2*i));
tempc = mb;
tempc(i)=0;
end

%for val_c
if val_b + dis(gra(4,:), x(2*i-1:2*i)) >= val_d + dis(gra(8,:),x(2*i-1:2*i))
tempd = md;
val_d_t = val_d + dis(gra(8,:),x(2*i-1:2*i));
tempd(i)=1;
else
val_d_t = val_b + dis(gra(4,:),x(2*i-1:2*i));
tempd = mb;
tempd(i)=1;
end
val_a =val_a_t;
val_b =val_b_t;
val_c =val_c_t;
val_d =val_d_t;
ma = tempa;
mb = tempb;
mc = tempc;
md = tempd;
end

if val_a <= val_b
m = ma;
t = val_a;
else
m = mb;
t = val_a;
end

if val_c <= t
m = mc;
t =val_c;
end

if val_d <= t
m = md;
t = val_d;
end
%% 去掉最开始的00
m = m(3:s+2);


译码中间用到dis 函数，测量两段序列的海明距离

dis

%% function to compute distances between x,y in the form of [0,0]  [1,1]
function d = dis(x, y)
d = sum(xor(x,y));



加入噪声，绘制误比特率随snr变化的曲线
show error bit


errbit  = zeros(1,21);
for j = -5:15

y = awgn(x, j, 'measured');
for i = 1:2*N
if y(i) >=0.5
y(i)=1;
else
y(i)=0;
end
end

% 译码
m=viterbi_hard(y);
errbit(j+6) = sum(m~=code_in)/N ;

end

logerr = 10*log10(errbit);
plot( -5:15, logerr );

更多相关内容
• 维特比解码matlab代码使用维特比解码和路径度量的错误检测和纠正 该项目承担的任务是编写用于g1 = 110111和g2 = 111011的卷积（2、1、9）解码的MATLAB代码，并使用维特比算法和路径度量进行百分比错误检测和校正的...
• 维特比解码matlab代码Materl Viterbi解码器算法的实现 维特比算法作为卷积码的最大似然（ML）解码技术而闻名。 （n，k，m）维特比解码器中的路径存储单元负责跟踪与由路径度量单元指定的尚存路径相关联的信息位。 二...
• 维特比解码matlab代码维特比解码器 Matlab Viterbi解码器的实现此程序要求输入1位（nc），代码长度（l），约束长度（kc），要添加到编码后的代码字的错误数，然后程序生成给定长度的随机比特序列，使用随机生成...
• 维特比译码matlab代码
• 维特比解码matlab代码使用维特比算法解码卷积码 我可以说Python比MATLAB慢，并且比C语言慢得多
• 维特比解码matlab代码无线网络 Matlab项目：matlab中的Viterbi算法 目标：使用维特比算法解码卷积编码的数据。 使用的技术：MatLab及其编码和解码功能用于生成网格并比较原始序列和解码序列。 说明：最初以二进制...
• 卷积码是一种性能优越的信道编码，它的编码器和解码器都比较易于实现，同时还具有较强的纠错能力，这使得它的使用越来越广泛。
• 信息论课程的一套代码，附有详细的每一部分的代码说明文档。卷积编码采用自己写的代码，完成了217和319码。在MATLAB2016b环境下编译通过
• 维特比解码matlab代码CS5114维特比项目 2014年夏天，我在VT参加了CS5114，算法理论。 我们的第一个项目是分析动态编程算法，我选择研究用于解码卷积码的Viterbi算法。 该存储库包括该项目的结果。 内容 该存储库包含...
• 维特比解码matlab代码卷积通道编码和维特比解码器的实现 从头开始进行卷积通道编码和维特比解码器的MATLAB实现。 卷积编码器/解码器的实现可以使用任何首选的生成多项式。 除了信息速率r之外，所使用的生成多项式的...
• 维特比解码matlab代码完整的安全和无线通信系统项目 标题 用于已开发功能的应用程序之一（可以用于不同类型的网络协议） 卷积编码器，直接序列扩频，维特比解码器的安全加密和传输数据/消息的性能，非常强大的纠错...
• 维特比解码matlab代码YAMCRF-另一个Matlab条件随机场工具包 苏志良（c）2011版权所有。保留所有权利。 作者：苏志良 (zsu2 [at] buffalo [dot] edu) University at Buffalo, SUNY South China Agricultrual ...
• 维特比解码matlab代码施普林格细分代码 基于持续时间的HMM的心音分割代码 这是运行出版物中概述的心音分割算法的Matlab代码： D.Springer等人，“基于Logistic回归-基于HSMM的心音分割”，IEEE Trans。 生物医学。 ...
• 维特比解码matlab代码霍夫曼和卷积编码器 霍夫曼编码器和解码器以及通道卷积编码器和维特比解码器的Matlab实现。 main.m文件使用不同的SNR级别测试实现的性能，并将其与没有卷积编码的通道的性能进行比较。
• 维特比解码matlab代码卷积编码和维特比编码 问题陈述 用于g1 = 110111和g2 = 11101101的卷积码（ 2、1、10 ）解码的MATLAB代码，并使用维特比解码和路径度量分析错误检测和纠正的百分比。 项目贡献者 戈文德·杰文...
• 采用MATLAB自编函数对卷积码以及维特比译码进行仿真，且对其性能进行分析。由于卷积码有性能floor,编码增益随信噪比降低而体现不明显。 　1.引言 　卷积码的编码器是由一个有k位输入、n位输出，且具有m位移位...
• 完整的维特比译码的教学资料，并结合详细的实例进行解释，看了就能完全掌握维特比译码的知识
• matlab实现的维特比译码，已用数字通信第二版实例验证通过
• 利用matlab的communication toolbox实现AWGN信道下采用QPSK调制和卷积码编码，然后接收端采用维特比译码并且采用硬判决的系统最终得到的误码率曲线，并且采用BERtool工具将其与理论值进行比较。
• FPGA信号处理系列文章——卷积编码与维特比译码前言卷积编码的概念和定义卷积编码算法理解卷积编码FPGA实现 前言 之前在 北斗二号导航系统中的RDSS系统以及Galileo系统中接触到了其电文是有卷积编码的，这里我们...

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

# 前言

之前在北斗导航系统中的RDSS系统以及Galileo系统中接触到了其电文是有卷积编码的，这里我们就来回顾一下卷积编码和维特比译码的相关内容

# 卷积编码的概念和定义

通常以 (n, k, L) 来描述卷积编码
• k 表示编码器的 输入码元数
• n 表示编码器的 输出码元数
• L表示编码器的 约束长度

由输入的 k 个信息比特，得到 n 个编码结果，所以 编码效率 = k/n。约束长度 L 的意思是，编码结果不仅取决于当前输入的 k 比特信息，还取决于前 (L-1) 段时间内的信息位。在 k = 1 的条件下，编码器需要的 移位寄存器级数 m = L - 1。

Galileo中采用的卷积码参数为：

因此：Galileo 卷积编码 (n, k, L) 描述为（2,1,7）

这里有个特殊的地方：
第二分支出来的数据要反向，这是Galileo特有的

还有一个特点是，（2,1,7）卷积码的卷尾是6个0。

# 卷积编码算法理解

ICD的附录中有一个这样的例子

编码前： 244个bit

11111111 11110000 11001100 10101010 00000000 00001111 00110011 01010101
11100011 11101100 11011111 10001010 00011100 00010011 00100000 01110101
01010101 01100001 01100010 01100011 10101010 10011110 10011101 10011100
00011100 00100011 00001101 00111001 11011100 11101100 0000


编码后：

10001100 00011010 10101010 01110011 00110001 01011010 01101111 01011001
01111000 10010101 01010101 10001100 11001110 10100101 10010000 10100110
10000100 00000010 00000010 00011011 10011001 11101011 01011100 00011000
10111011 11111101 11111101 11100100 01100110 00010100 10100011 11100111
10100110 01100110 01101011 00010000 00011100 00011101 01010001 10101110
00111001 11101001 10010100 11101111 11100011 11100010 10101110 01010001
11111101 00111101 11110000 10101001 11011000 00110010 00111111 10000100
10010010 10010110 00100100 10101011 01001110


正好，我可以这个数据作为算法验证的数据

卷积编码的matlab程序，注意标红部分为额外取反处理

a = '1111111111110000110011001010101000000000000011110011001101010101111000111110110011011111100010100001110000010011001000000111010101010101011000010110001001100011101010101001111010011101100111000001110000100011000011010011100111011100111011000000';

for i = 1 : 244
msg(i)=str2num(a(i));
end

trel = poly2trellis(7,[171,133]); % Define trellis
code = convenc(msg,trel); % Encode
code(2:2:end) = not(code(2:2:end))


得出的结果与ICD上的一致。

# 卷积编码FPGA实现

卷积编码的verilog实现也非常简单，也就不多说了，直接上代码，与上面Galileo 卷积编码对应

module Convolutional_Coding_2_1_7(
input clk,
input encode_rst,
input data_in,
input nd,
output rdy,
output [1:0]data_out_v,
output data_out_s,
output rdy_s
);

wire [6:0]PolyA;
wire [6:0]PolyB;
reg [6:0]wA;
reg [6:0]wB;
reg [6:0]ShReg;
reg nd_d1;
reg nd_d2;
reg nd_d3;

assign PolyA = 7'b1_011_011;
assign PolyB = 7'b1_111_001;

always @(posedge clk)
if(encode_rst)
begin
ShReg <= 0;
end
else if(nd)
begin
ShReg[6] <= data_in;
ShReg[5:0] <= ShReg[6:1];
end

always @(posedge clk)
if(encode_rst)
begin
wA <= 0;
wB <= 0;
end
else
begin
wA <= PolyA & ShReg;
wB <= PolyB & ShReg;
end

assign data_out_v[0] = ^wA;
assign data_out_v[1] = ^wB;

always @(posedge clk)
if(encode_rst)
begin
nd_d1 <= 0;
nd_d2 <= 0;
nd_d3 <= 0;
end
else
begin
nd_d1 <= nd;
nd_d2 <= nd_d1;
nd_d3 <= nd_d2;
end

assign rdy = nd_d2;

assign rdy_s = nd_d2 | nd_d3;
assign data_out_s = nd_d2 ? data_out_v[0] : ~data_out_v[1];


# 维特比译码的说明

卷积编码解码的过程就是维特比译码，编码容易解码难。难道什么程度，难道xilinx 的viterbi ip都是要收费的就知道这个有多复杂了。

但我们可以先利用matlab提供的函数搞清楚先。

维特比解码分为硬判决和软判决

根据之前的经验，软判决一般采用3,4,5bit输入即可，再增加输入性能增加不明显，而且资源显著增加。性能软判决相比较硬判决灵敏度可能有个3-5dB的提升。（精确的理论分析见一些论文，实际工程与设计好坏也有关系）

# 维特比译码算法理解

维特比的matlab算法如下，注意几个地方
1、解码后的前42个为无效，与回溯深度tblen=42有关。
2、如果译码不是以数据流的形式进行，再数据结尾需要再补齐tblen*2个随机数据才行，不然会少输出tblen个解码数据

硬判决：

clc;clear all;close all;

a = '1111111111110000110011001010101000000000000011110011001101010101111000111110110011011111100010100001110000010011001000000111010101010101011000010110001001100011101010101001111010011101100111000001110000100011000011010011100111011100111011000000';

for i = 1 : 244
msg(i)=str2num(a(i));
end

%编码过程
trel = poly2trellis(7,[171,133]); % Define trellis
code = convenc(msg,trel); % Encode
code(488+1:488+42*2) = randi([0 1],42*2,1)';
code(2:2:end) = not(code(2:2:end))

%译码过程

tblen = 42;
code(2:2:end) = not(code(2:2:end))
decoded = vitdec(code,trel,tblen,'cont','hard');

plot(msg == decoded(43:end))

%decoded的前42个为无效，与tblen=42有关，所以输入需要多输入tblen*2个编码数据才行


软判决：

clc;clear all;close all;

a = '1111111111110000110011001010101000000000000011110011001101010101111000111110110011011111100010100001110000010011001000000111010101010101011000010110001001100011101010101001111010011101100111000001110000100011000011010011100111011100111011000000';

for i = 1 : 244
msg(i)=str2num(a(i));
end

tblen = 48;

trel = poly2trellis(7,[171,133]); % Define trellis
code = convenc(msg,trel); % Encode
code(488+1:488+tblen*2) = randi([0 1],tblen*2,1)';

%code(1:5) = 0;
%code(2:2:end) = not(code(2:2:end))

decoded1 = vitdec(code,trel,tblen,'cont','hard');
%decoded的前42个为无效，与tblen有关

% Map "0" bit to 1.0 and "1" bit to -1.0.  Also add AWGN.
ucode = real(awgn(1-2*code, 2, 'measured'));

% Soft decision decoding with quantized inputs
[x,qcode] = quantiz(ucode,[1/(2^(8-1))-1:1/(2^(8-1)):1-1/(2^(8-1))],...
2^8-1:-1:0); % Values in qcode are between 0 and 2^8-1.
decoded2 = vitdec(qcode,trel,tblen,'cont','soft',8);

[n2,r2] = biterr(decoded2(tblen+1:end),msg(1:end));
disp(['The bit error rates are:   ',num2str([r2])])


软判决，软判决对编码后的数据加噪声，来体现它的性能

clc;clear all;close all;

a = '1111111111110000110011001010101000000000000011110011001101010101111000111110110011011111100010100001110000010011001000000111010101010101011000010110001001100011101010101001111010011101100111000001110000100011000011010011100111011100111011000000';

for i = 1 : 244
msg(i)=str2num(a(i));
end

%编码过程
tblen = 48;

trel = poly2trellis(7,[171,133]); % Define trellis
code = convenc(msg,trel); % Encode
code(2:2:end) = not(code(2:2:end)) %galileo 定义独有
code(488+1:488+tblen*2) = randi([0 1],tblen*2,1)';

code(2:2:end) = not(code(2:2:end))

%硬判决
decoded1 = vitdec(code,trel,tblen,'cont','hard');
%decoded的前42个为无效，与tblen有关

% Map "0" bit to 1.0 and "1" bit to -1.0.  Also add AWGN.
ucode = real(awgn(1-2*code, 3, 'measured'));

% Soft decision decoding with quantized inputs
[x,qcode] = quantiz(ucode,[1/(2^(8-1))-1:1/(2^(8-1)):1-1/(2^(8-1))],...
2^8-1:-1:0); % Values in qcode are between 0 and 2^8-1.

%软判决
decoded2 = vitdec(qcode,trel,tblen,'cont','soft',8);

[n2,r2] = biterr(decoded2(tblen+1:end),msg(1:end));
disp(['The bit error rates are:   ',num2str([r2])])


# 维特比译码的FPGA实现

1、xilinx ip的维特比译码器
不交钱就不给你玩，不过幸好我以前玩过。我们可以看看里面的一些设置

可以看到，也就这些参数的设置，所以如果能用ip的话就挺简单的

2、opencores网站上有一个开源项目

这个我也玩过，经过之前的验证，发现如果使用这个开源代码的话需要软判决位宽为6才能勉强和xilinx ip核的性能匹配上。有兴趣的小伙伴可以直接下载去验证下。

3、自己写

后面准备有空的时候，重新再去理解一下这个译码过程。当年虽然也花了些时间去理解，但始终云里雾里的感觉，忘的也差不多了，也许现在也不一定能看明白，到时候希望能搞明白点吧。有些东西可能过了几年回过头来看又有新的感悟呢

展开全文
• % 将经过信道的编码序列转换为以需要译码的个数为行数，以n为列数的矩阵 survivor_state=zeros(number_of_states,depth_of_trellis+1); % 4行，7列 (初始化) % 开始无尾信道输出解码 for i=1:depth_of_trellis-L+1 ...

### 213

function [mybit,decoder_output_table,BER,right_channel_output,channel_output_table,cumulated_metric_table,survivor_state]=viterbi_2_1_3
changdu=input('Please input the length of the code sequence\n\n');
mybit=fix(rand(1,changdu)+0.5);

%  (2,1,3)卷积编码
G=[1 1 1;1 0 1];                                                                          %  (2,1,3)生成矩阵G=[1 1 1;1 0 1];
k=1;
right_channel_output=cnv_encd(mybit);
lgth=length(right_channel_output);
BER=zeros(1,lgth);
cumulated_metric_table=zeros(1,lgth);
decoder_output_table=zeros(lgth,lgth/2-2);
channel_output_table=zeros(lgth,lgth);

%  模拟信道误码
for z=1:lgth                                                                             % 产生随机误码

%    channel=zeros(1,lgth);                                                              % 通过改变lgth的值 ,调整误码个数
%    for i=1:z:lgth
%        channel(i)=1;
%    end

if z==1
channel=zeros(1,lgth);
else
channel=fix(rand(1,lgth)+0.2);                                                       % 0.5可改
end
channel_output=mod(right_channel_output+channel,2);
channel_output_table(z,:)=channel_output;                                                % 数据转移
n=size(G,1);                                                                             % n=G的行数

%  检查大小
if rem(size(G,2),k)~=0
error('Size of G and k do not agree')                                                % 存在合理的约束长度
end
if rem(size(channel_output,2),n)~=0
error('channel output not of the right size')                                        % 经过信道后，进入译码器前，卷积编码的个数合理
end
L=size(G,2)/k;                                                                           % 约束长度＝寄存器个数+1
number_of_states=2^((L-1)*k);                                                            % 状态数

%  产生状态转移矩阵、输出矩阵和输入矩阵
for j=0:number_of_states-1                                                               % 对(2,1,3) j=0,1,2,3
for l=0:2^k-1                                                                        % l=0,1
[next_state,memory_contents]=nxt_stat(j,l,L,k);
input(j+1,next_state+1)=l;
branch_output=rem(memory_contents*G',2);                                         % 状态图上根据寄存器状态和input作出的编码
nextstate(j+1,l+1)=next_state;                                                   % 生成矩阵，含有十进制表示的下一状态
output(j+1,l+1)=bin2deci(branch_output);                                         % 生成矩阵，含有十进制表示的卷积编码
end
end
state_metric=zeros(number_of_states,2);                                                  % 4行，2列 (初始化)
depth_of_trellis=length(channel_output)/n;                                               % 需要译码的个数   设为6
channel_output_matrix=reshape(channel_output,n,depth_of_trellis);                        % 将经过信道的编码序列转换为以需要译码的个数为行数，以n为列数的矩阵
survivor_state=zeros(number_of_states,depth_of_trellis+1);                               % 4行，7列 (初始化)

%  开始无尾信道输出解码
for i=1:depth_of_trellis-L+1                                                             % i=1,2,3,4
flag=zeros(1,number_of_states);                                                      % 1行，4列 (初始化)
if i<=L
step=2^((L-i)*k);                                                                % i=1,2,3  step= 4,2,1
else
step=1;                                                                          % i=4  step=1
end
for j=0:step:number_of_states-1                                                      % i=1 j=0;i=2 j=0, 2;i=3,4 j=0,1,2,3
for l=0:2^k-1                                                                    % l=0,1
branch_metric=0;                                                             % 初始化
binary_output=deci2bin(output(j+1,l+1),n);                                   % 将此时刻十进制的编码结果用二进制表示
for ll=1:n                                                                   % ll=1,2
branch_metric=branch_metric+metric(channel_output_matrix(ll,i),...
binary_output(ll));                                                  % 将此时刻经过信道的编码与理论编码比较，得到分支度量值
end
if ((state_metric(nextstate(j+1,l+1)+1,2)>state_metric(j+1,1)...
+branch_metric)|flag(nextstate(j+1,l+1)+1)==0)
state_metric(nextstate(j+1,l+1)+1,2)=state_metric(j+1,1)+branch_metric;
survivor_state(nextstate(j+1,l+1)+1,i+1)=j;
flag(nextstate(j+1,l+1)+1)=1;
end
end
end
state_metric=state_metric(:,2:-1:1);
end

%  开始有尾信道输出解码
for i=depth_of_trellis-L+2:depth_of_trellis
flag=zeros(1,number_of_states);
last_stop=number_of_states/(2^((i-depth_of_trellis+L-2)*k));
for j=0:last_stop-1
branch_metric=0;
binary_output=deci2bin(output(j+1,1),n);
for ll=1:n
branch_metric=branch_metric+metric(channel_output_matrix(ll,i),binary_output(ll));
end
if ((state_metric(nextstate(j+1,1)+1,2)>state_metric(j+1,1)...
+branch_metric)|flag(nextstate(j+1,1)+1)==0)
state_metric(nextstate(j+1,1)+1,2)=state_metric(j+1,1)+branch_metric;
survivor_state(nextstate(j+1,1)+1,i+1)=j;
flag(nextstate(j+1,1)+1)=1;
end
end
state_metric=state_metric(:,2:-1:1);
end

%  从最佳路径中产生解码
state_sequence=zeros(1,depth_of_trellis+1);
state_sequence(1,depth_of_trellis)=survivor_state(1,depth_of_trellis+1);
for i=1:depth_of_trellis
state_sequence(1,depth_of_trellis-i+1)=survivor_state((state_sequence(1,depth_of_trellis+2-i)...
+1),depth_of_trellis-i+2);
end
decoder_output_matrix=zeros(k,depth_of_trellis-L+1);
for i=1:depth_of_trellis-L+1
dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
dec_output_bin=deci2bin(dec_output_deci,k);
decoder_output_matrix(:,i)=dec_output_bin(k:-1:1)';
end
decoder_output=reshape(decoder_output_matrix,1,k*(depth_of_trellis-L+1));
decoder_output_table(z,:)=decoder_output;                                                % 数据转移
cumulated_metric=state_metric(1,1);
cumulated_metric_table(z)=cumulated_metric;                                              % 数据转移
BER(z)=Bit_Error_Rate(mybit,decoder_output);
%  BER1(z)=Bit_Error_Rate(right_channel_output,channel_output)
end

%----------------------------------------subfunction Bit_Error_Rate-------------------------------------------------

function y=Bit_Error_Rate(a,b)
y=nnz(mod(a+b,2))/length(a);

%-------------------------------------------subfunction cnv_encd----------------------------------------------------

function the_output1=cnv_encd(mybit)       %  the_output1为完全正确译码，编码输出结果多出2*移位寄存器个数；  the_output2为不完全译码，译码结果与编码之前的原始信息差移位寄存器的个数
g=[1 1 1 0 1 1];
n=length(mybit);
juzhen=zeros(n,(n-1)*2+6);                 %  多余个数＝2*约束长度-2=2*移位寄存器个数
for i=1:n
juzhen(i,(i-1)*2+1:(i-1)*2+6)=g;
end
the_output=mybit*juzhen;
the_output=mod(the_output,2);
the_output1=the_output;                    %  需要时可以把分号去掉，看完全编码和不完全编码的结果
the_output2=the_output(1,1:n*2);

% survivor_state的每一行对应每一个状态，每一列对应每一个时刻，它的数值为从前一时刻到这一时刻、进入到所在的行状态所保留的幸存路径的来源状态
% cumulated_metric为译码结果所对应的编码序列与信道输入含误码序列的汉明距离

### 217

clc;
clear;
close all

mybit=fix(rand(1,215)+0.5)
lgth = 215;
%  (2,1,7)卷积编码
G=[1 1 1 1 0 0 1;1 0 1 1 0 1 1];                                                       %  (2,1,7)生成矩阵G=[1 1 1 1 0 0 1;1 0 1 1 0 1 1];
k=1;
right_channel_output=cnv_encd(mybit);
BER=zeros(1,lgth);
cumulated_metric_table=zeros(1,lgth);
decoder_output_table=zeros(lgth,lgth/2-6);
channel_output_table=zeros(lgth,lgth);

%  模拟信道误码
for z=1:lgth                                                                             % 产生随机误码

%    channel=zeros(1,lgth);                                                              % 通过改变lgth的值 ,调整误码个数
%    for i=1:z:lgth
%        channel(i)=1;
%    end

if z==1
channel=zeros(1,lgth);
else
channel=fix(rand(1,lgth)+0.2);                                                       % 0.5可改
end
channel_output=mod(right_channel_output+channel,2);
channel_output_table(z,:)=channel_output;                                                % 数据转移
n=size(G,1);                                                                             % n=G的行数

%  检查大小
if rem(size(G,2),k)~=0
error('Size of G and k do not agree')                                                % 存在合理的约束长度
end
if rem(size(channel_output,2),n)~=0
error('channel output not of the right size')                                        % 经过信道后，进入译码器前，卷积编码的个数合理
end
L=size(G,2)/k;                                                                           % 约束长度＝寄存器个数+1
number_of_states=2^((L-1)*k);                                                            % 状态数

%  产生状态转移矩阵、输出矩阵和输入矩阵
for j=0:number_of_states-1                                                               % 对(2,1,3) j=0,1,2,3
for l=0:2^k-1                                                                        % l=0,1
[next_state,memory_contents]=nxt_stat(j,l,L,k);
input(j+1,next_state+1)=l;
branch_output=rem(memory_contents*G',2);                                         % 状态图上根据寄存器状态和input作出的编码
nextstate(j+1,l+1)=next_state;                                                   % 生成矩阵，含有十进制表示的下一状态
output(j+1,l+1)=bin2deci(branch_output);                                         % 生成矩阵，含有十进制表示的卷积编码
end
end
state_metric=zeros(number_of_states,2);                                                  % 4行，2列 (初始化)
depth_of_trellis=length(channel_output)/n;                                               % 需要译码的个数   设为6
channel_output_matrix=reshape(channel_output,n,depth_of_trellis);                        % 将经过信道的编码序列转换为以需要译码的个数为行数，以n为列数的矩阵
survivor_state=zeros(number_of_states,depth_of_trellis+1);                               % 4行，7列 (初始化)

%  开始无尾信道输出解码
for i=1:depth_of_trellis-L+1                                                             % i=1,2,3,4
flag=zeros(1,number_of_states);                                                      % 1行，4列 (初始化)
if i<=L
step=2^((L-i)*k);                                                                % i=1,2,3  step= 4,2,1
else
step=1;                                                                          % i=4  step=1
end
for j=0:step:number_of_states-1                                                      % i=1 j=0;i=2 j=0, 2;i=3,4 j=0,1,2,3
for l=0:2^k-1                                                                    % l=0,1
branch_metric=0;                                                             % 初始化
binary_output=deci2bin(output(j+1,l+1),n);                                   % 将此时刻十进制的编码结果用二进制表示
for ll=1:n                                                                   % ll=1,2
branch_metric=branch_metric+metric(channel_output_matrix(ll,i),...
binary_output(ll));                                                  % 将此时刻经过信道的编码与理论编码比较，得到分支度量值
end
if ((state_metric(nextstate(j+1,l+1)+1,2)>state_metric(j+1,1)...
+branch_metric)|flag(nextstate(j+1,l+1)+1)==0)
state_metric(nextstate(j+1,l+1)+1,2)=state_metric(j+1,1)+branch_metric;
survivor_state(nextstate(j+1,l+1)+1,i+1)=j;
flag(nextstate(j+1,l+1)+1)=1;
end
end
end
state_metric=state_metric(:,2:-1:1);
end

%  开始有尾信道输出解码
for i=depth_of_trellis-L+2:depth_of_trellis
flag=zeros(1,number_of_states);
last_stop=number_of_states/(2^((i-depth_of_trellis+L-2)*k));
for j=0:last_stop-1
branch_metric=0;
binary_output=deci2bin(output(j+1,1),n);
for ll=1:n
branch_metric=branch_metric+metric(channel_output_matrix(ll,i),binary_output(ll));
end
if ((state_metric(nextstate(j+1,1)+1,2)>state_metric(j+1,1)...
+branch_metric)|flag(nextstate(j+1,1)+1)==0)
state_metric(nextstate(j+1,1)+1,2)=state_metric(j+1,1)+branch_metric;
survivor_state(nextstate(j+1,1)+1,i+1)=j;
flag(nextstate(j+1,1)+1)=1;
end
end
state_metric=state_metric(:,2:-1:1);
end

%  从最佳路径中产生解码
state_sequence=zeros(1,depth_of_trellis+1);
state_sequence(1,depth_of_trellis)=survivor_state(1,depth_of_trellis+1);
for i=1:depth_of_trellis
state_sequence(1,depth_of_trellis-i+1)=survivor_state((state_sequence(1,depth_of_trellis+2-i)...
+1),depth_of_trellis-i+2);
end
decoder_output_matrix=zeros(k,depth_of_trellis-L+1);
for i=1:depth_of_trellis-L+1
dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);
dec_output_bin=deci2bin(dec_output_deci,k);
decoder_output_matrix(:,i)=dec_output_bin(k:-1:1)';
end
decoder_output=reshape(decoder_output_matrix,1,k*(depth_of_trellis-L+1));
decoder_output_table(z,:)=decoder_output;                                                % 数据转移
cumulated_metric=state_metric(1,1);
cumulated_metric_table(z)=cumulated_metric;                                              % 数据转移

end

%----------------------------------------subfunction Bit_Error_Rate-------------------------------------------------

%-------------------------------------------subfunction cnv_encd----------------------------------------------------

D30

展开全文
• 卷积编码及维特比译码-卷积编码及维特比译码.rar 卷积编码及维特比译码相关文献 1.bmp
• BCC卷积码的维特比译码算法_matlab_demo 自己参照原理写的维特比译码算法 对（2,1,7）217卷积码数据进行加高斯白噪声，对比不同信噪比下的维特比译码误比特率 使用：在matlab直接跑CC_SimTest.m可运行
• 首先是演示代码，首先输入输入信号，用m序列加扰(演示里只用16位长度的m序列)，再1/3卷积编码，接收端先维特比解码，再解扰，得到原始信号 Demo.m   function demo input=[1 0 1 1 0 1 0 0 0 1 1];%输入信号 ...

首先是演示代码，首先输入输入信号，用m序列加扰(演示里只用16位长度的m序列)，再1/3卷积编码，接收端先维特比解码，再解扰，得到原始信号

Demo.m

function demo

input=[1 0 1 1 0 1 0 0 0 1 1];%输入信号
subplot(2,3,1);
drawSig(input);
input_r=scramble(input);%加扰
subplot(2,3,2);
drawSig(input_r);
Y=coder(input_r);%进行卷积编码
subplot(2,3,3);
drawSig(Y);
O=decoder(Y);%维特比解码
subplot(2,3,4);
drawSig(O);
Res=scramble(O);%解扰
subplot(2,3,5);
drawSig(Res);

end

得到图片

mxulie.m

%产生16位m序列
function Y=mxulie()
F(16)=0;
a3=1;
a2=0;
a1=0;
a0=0;
for i=1:16
F(i)=a0;
a0=a1;
a1=a2;
a2=a3;
a3=mod((a3+a0),2);
end
Y=F;
end


coder.m

%1/3,3卷积编码
function Y = coder( X )
len=length(X);
T=zeros(1,len+3,'int8');
F=zeros(1,len*3,'int8');
for i=4:len+3
T(i)=X(i-3);
end
for i=1:len
F(3*(i-1)+1)=T(i+3);
F(3*(i-1)+2)=xor(T(i+3),T(i+1));
F(3*(i-1)+3)=xor(xor(T(i+3),T(i+2)),T(i+1));
end
Y=F;
end

decoder.m

%1/3,3维特比卷积解码
function Y = decoder( X )
len=length(X);
Path=zeros(4,len/3,'int8');%4条幸存路径
Path_t=zeros(4,len/3,'int8');
Da=0;Db=0;Dc=0;Dd=0;%抵达a,b,c,d总汉明距离
Pa=1;Pb=2;Pc=3;Pd=4;%a,b,c,d路径指针
T(9)=0;
for i=1:9
T(i)=X(i);
end
%a
tp1=dist(T,[0 0 0 0 0 0 0 0 0],9);
tp2=dist(T,[1 1 1 0 0 1 0 1 1],9);
if(tp1<tp2)
Da=tp1;
Path(1,1)=0;Path(1,2)=0;Path(1,3)=0;
else
Da=tp2;
Path(1,1)=1;Path(1,2)=0;Path(1,3)=0;
end
%b
tp1=dist(T,[0 0 0 0 0 0 1 1 1],9);
tp2=dist(T,[1 1 1 0 0 1 1 0 0],9);
if(tp1<tp2)
Db=tp1;
Path(2,1)=0;Path(2,2)=0;Path(2,3)=1;
else
Db=tp2;
Path(2,1)=1;Path(2,2)=0;Path(2,3)=1;
end
%c
tp1=dist(T,[0 0 0 1 1 1 0 0 1],9);
tp2=dist(T,[1 1 1 1 1 0 0 1 0],9);
if(tp1<tp2)
Dc=tp1;
Path(3,1)=0;Path(3,2)=1;Path(3,3)=0;
else
Dc=tp2;
Path(3,1)=1;Path(3,2)=1;Path(3,3)=0;
end
%d
tp1=dist(T,[0 0 0 1 1 1 1 1 0],9);
tp2=dist(T,[1 1 1 1 1 0 1 0 1],9);
if(tp1<tp2)
Dd=tp1;
Path(4,1)=0;Path(4,2)=1;Path(4,3)=1;
else
Dd=tp2;
Path(4,1)=1;Path(4,2)=1;Path(4,3)=1;
end
%迭代
Dat=0;Dbt=0;Dct=0;Ddt=0;%交换缓冲
fga=0;fgb=0;fgc=0;fgd=0;%路径标志
rmSz=int32((len-9)/3);
for i=1:rmSz
T(1)=X(9+(i-1)*3+1);
T(2)=X(9+(i-1)*3+2);
T(3)=X(9+(i-1)*3+3);
%a
tp1=dist(T,[0 0 0],3)+Da;
tp2=dist(T,[0 1 1],3)+Dc;
if(tp1<tp2)
Dat=tp1;
fga=0;
else
Dat=tp2;
fga=1;
end
%b
tp1=dist(T,[1 1 1],3)+Da;
tp2=dist(T,[1 0 0],3)+Dc;
if(tp1<tp2)
Dbt=tp1;
fgb=0;
else
Dbt=tp2;
fgb=1;
end
%c
tp1=dist(T,[0 0 1],3)+Db;
tp2=dist(T,[0 1 0],3)+Dd;
if(tp1<tp2)
Dct=tp1;
fgc=0;
else
Dct=tp2;
fgc=1;
end
%d
tp1=dist(T,[1 1 0],3)+Db;
tp2=dist(T,[1 0 1],3)+Dd;
if(tp1<tp2)
Ddt=tp1;
fgd=0;
else
Ddt=tp2;
fgd=1;
end
%更新幸存路径
%a
if(fga==0)
Path_t(Pa,:)=Path(Pa,:);
Path_t(Pa,3+i)=0;
Da=Dat;
else
Path_t(Pa,:)=Path(Pc,:);
Path_t(Pa,3+i)=0;
Da=Dat;
end
%b
if(fgb==0)
Path_t(Pb,:)=Path(Pa,:);
Path_t(Pb,3+i)=1;
Db=Dbt;
else
Path_t(Pb,:)=Path(Pc,:);
Path_t(Pb,3+i)=1;
Db=Dbt;
end
%c
if(fgc==0)
Path_t(Pc,:)=Path(Pb,:);
Path_t(Pc,3+i)=0;
Dc=Dct;
else
Path_t(Pc,:)=Path(Pd,:);
Path_t(Pc,3+i)=0;
Dc=Dct;
end
%d
if(fgd==0)
Path_t(Pd,:)=Path(Pb,:);
Path_t(Pd,3+i)=1;
Dd=Ddt;
else
Path_t(Pd,:)=Path(Pd,:);
Path_t(Pd,3+i)=1;
Dd=Ddt;
end
%
Path(Pa,:)=Path_t(Pa,:);
Path(Pb,:)=Path_t(Pb,:);
Path(Pc,:)=Path_t(Pc,:);
Path(Pd,:)=Path_t(Pd,:);
end
k=min([Da Db Dc Dd]);
if(k==Da)
Y=Path(Pa,:);
end
if(k==Db)
Y=Path(Pb,:);
end
if(k==Dc)
Y=Path(Pc,:);
end
if(k==Dd)
Y=Path(Pd,:);
end
end



dist.m

%汉明距离
function d = dist(X1,X2,n)
sum=0;
for i=1:n
if X1(i)~=X2(i)
sum=sum+1;
end
end
d=sum;
end

scramble.m

%加扰
function Y = scramble( X )
n=length(X);
M=mxulie();
F(n)=0;
for i=1:n;
F(i)=xor(X(i),M(i));
end
Y=F;
end

drawSig.m

%用于绘制信号图
function drawSig(X)
n=length(X);
G(n*100)=0;
for i=1:n
if X(i)==0
for j=1:100
G((i-1)*100+j)=0;
end
else
for j=1:100
G((i-1)*100+j)=1;
end
end
end
plot(G);
axis([0 n*100 -2 2]);
end

展开全文
• ## 维特比译码详解

千次阅读 2022-03-10 12:26:18
维特比译码算法详细介绍
• MATLAB代码实现了卷积编码，BPSK调制，以及维特比译码，其中主要实现维特比译码，附加了Word文档说明。
• 2 维特比译码原理Viterbi译码是卷积码的最大似然译码算法，是一种实用化的概率算法。它的基本思想是把已接收序列与所有可能的发送序列作比较，选择其中码距最小的一个序列作为发送序列。从图2的卷积码网格图可以看出...
• viterbi译码算法是一种卷积码的解码算法。优点不说了。缺点就是随着约束长度的增加算法的复杂度增加很快。约束长度N为7时要比较的路径就有64条，为8时路径变为128条。 (2<先说编码(举例约束长度为7)：编码器7个...

...

matlab 订阅