turbo系统 matlab 仿真代码 参考Dr.Wuyufei 有详细注释

clear all; close all; clc; % Copyright Nov. 1998 Yufei Wu % MPRG lab, Virginia Tech. % for academic use only % 代码修改：李亚坤 % 云南大学 % liyakun1990@gmail.com % 版本日期：11月2日 % 版本日期：11月3日 % 版本日期：11月5日 % 版本日期：11月6日 % 版本日期：11月7日 % 版本日期：11月8日 % 张敏，要做Turbo，感觉是绕不开Dr.wuyufei的数学架构的 % **************Dr. wuyufei的典型Turbo系统主函数******************************** % This script simulates the classical turbo encoding-decoding system. % It simulates parallel concatenated convolutional codes. % Two component rate 1/2 RSC (Recursive Systematic Convolutional) component encoders are assumed. % % First encoder is terminated with tails bits. (Info + tail) bits are scrambled and passed to % the second encoder, while second encoder is left open without tail bits of itself. % % Random information bits are modulated into +1/-1, and transmitted through a AWGN channel. % Interleavers are randomly generated for each frame. % While it's UNECESSARY in our first version! % % Log-MAP algorithm without quantization or approximation is used. % By making use of ln(e^x+e^y) = max(x,y) + ln(1+e^(-abs(x-y))), % the Log-MAP can be simplified with a look-up table for the correction function. % If use approximation ln(e^x+e^y) = max(x,y), it becomes MAX-Log-MAP. % 将结果写入txt diary turbo_logmap.txt % 选择编码算法,0为Log-MAP,1为SOVA，默认为Log-MAP dec_alg = input(' Please enter the decoding algorithm. (0:Log-MAP, 1:SOVA) default 0(We choose 0 indeed) '); if isempty(dec_alg) dec_alg = 0; end % 设定帧长，默认帧长400 L_total = input(' Please enter the frame size (= info + tail, default: 400) '); if isempty(L_total) L_total = 400; % 包括信息和附加码 end % 码元生成矩阵 g = input(' Please enter code generator: ( default: g = [1 1 1; 1 0 1 ] ) '); if isempty(g) g = [ 1 1 1; 1 0 1 ]; end %g = [1 1 0 1; 1 1 1 1]; %g = [1 1 1 1 1; 1 0 0 0 1]; [n,K] = size(g); %每一时刻送至编码器的输入信息元为k个，相应的编码输出码元为n个 m = K - 1;%m为移位寄存器的个数 nstates = 2^m;%卷积码的编码寄存器的状态 %puncture = 0,删余，码元速率1/2; %puncture = 1,不删余，码元速率不变 puncture = input(' Please choose punctured / unpunctured (0/1): default 0 '); if isempty(puncture) puncture = 0; end % 码率 rate = 1/(2+puncture); % 信道衰落幅度; a=1 为在高斯白噪声信道 a = 1; % 每帧迭代次数，默认为5次 niter = input(' Please enter number of iterations for each frame: default 5 '); if isempty(niter) niter = 5; end % 定义多少误帧数作为停止准则，默认为15 ferrlim = input(' Please enter number of frame errors to terminate: default 15 '); if isempty(ferrlim) ferrlim = 15; end %输入信噪比，默认2db EbN0db = input(' Please enter Eb/N0 in dB : default [2.0] '); if isempty(EbN0db) EbN0db = [2.0]; end %对TXT文档的记录输入 fprintf('\n\n----------------------------------------------------\n'); if dec_alg == 0 fprintf(' === Log-MAP decoder === \n'); else fprintf(' === SOVA decoder === \n'); end fprintf(' Frame size = %6d\n',L_total); fprintf(' code generator: \n'); for i = 1:n for j = 1:K fprintf( '%6d', g(i,j)); end fprintf('\n'); end if puncture==0 fprintf(' Punctured, code rate = 1/2 \n'); else fprintf(' Unpunctured, code rate = 1/3 \n'); end fprintf(' iteration number = %6d\n', niter); fprintf(' terminate frame errors = %6d\n', ferrlim); fprintf(' Eb / N0 (dB) = '); for i = 1:length(EbN0db) fprintf('%10.2f',EbN0db(i)); end fprintf('\n----------------------------------------------------\n\n'); fprintf('+ + + + Please be patient. Wait a while to get the result. + + + +\n'); for nEN = 1:length(EbN0db) en = 10^(EbN0db(nEN)/10); % 将噪声单位从分贝转换为十进制数 L_c = 4*a*en*rate; % 信道置信度Lc % 对于噪声服从分布N（0,No/2）的AWGN信道，信道置信度定义为... sigma = 1/sqrt(2*rate*en); % AWGN噪声标准差 % 误帧nferrr计数器和误比特errs计数器清零 errs(nEN,1:niter) = zeros(1,niter); nferr(nEN,1:niter) = zeros(1,niter); nframe = 0; % 发送帧计数清零 while nferr(nEN, niter)<ferrlim % 如果当该信噪比nEN下，误帧数小于最大误帧数 nframe = nframe + 1; % 发送帧数+1 x = round(rand(1, L_total-m)); % 对产生随机数进行就近取整，产生0-1码序列 %非零码组的持续长度(L_total)=非零信息码组的持续长度（L_info）+编码储存长度(m) [temp, alpha] = sort(rand(1,L_total));% 随机交织器的映射 %temp是原来的矩阵x按照行，每行从小到大重新排列得到的新矩阵。 %alpha告诉你重排的详细信息，也就是做了什么样的变动。 %例如b的第一行显示3 4 5 1 2，那么将原矩阵X的第一行的第3 4 5 1 2个元素取出来，顺次排列，就变成a矩阵的第一行。 en_output = encoderm( x, g, alpha, puncture ) ; % 编码器输出双极性码(+1/-1) r = en_output+sigma*randn(1,L_total*(2+puncture)); % 接收比特=编码器输出双极性码+标准差*随机正态分布噪声变量 % 此处就是模拟了高斯白噪声的信道环境 yk = demultiplex(r,alpha,puncture); % 分为多路，给译码器1和译码器2 % 扫描已接收比特 % 这行语句主要是接受比特收到信道衰减的影响 rec_s = 0.5*L_c*yk; % 初始化外部信息 L_e(1:L_total) = zeros(1,L_total); for iter = 1:niter % 译码器1 L_a(alpha) = L_e; % 一个先验信息 if dec_alg == 0 % log_map L_all = logmapo(rec_s(1,:), g, L_a, 1); % 完整信息 else L_all = sova0(rec_s(1,:), g, L_a, 1); % 完整信息 end % end if L_e = L_all - 2*rec_s(1,1:2:2*L_total)-L_a; % 外部信息 % 输出对数似然比信息=先验+系统+外部 % 求出外部信息用于另一个分量译码器，作为其先验信息 % 译码器2 L_a = L_e(alpha); % 一个先验信息 %L_a作为先验信息是另一个分量译码器生成的外部信息L_e经过交织后的对数似然比值 if dec_alg == 0 L_all = logmapo(rec_s(2,:), g, L_a, 2); % 完整信息 else L_all = sova0(rec_s(2,:), g, L_a, 2); % 完整信息 end % end if L_e = L_all - 2*rec_s(2,1:2:2*L_total) - L_a; % 外部信息 % 判决信息比特,解交织后存在xhat中，并由+1/-1转成0/1 % sign(x):if x>0,sign(x)=1; % if x=0,sign(x)=0; % if x<0,sign(x)=-1; xhat(alpha) = (sign(L_all)+1)/2; % 计算在当前循环中的误比特率 % 对比xhat和x err(iter) = length(find(xhat(1:L_total-m)~=x)); % 计算在当前循环中的误帧数 % 只要这一帧数据存在误比特,误帧数+1 if err(iter)>0 nferr(nEN,iter) = nferr(nEN,iter)+1; end end %iter,迭代循环+1 % 每次迭代循环，误比特数列举 errs(nEN,1:niter) = errs(nEN,1:niter) + err(1:niter); if rem(nframe,1)==0 | nferr(nEN, niter)==ferrlim % 误比特率 ber(nEN,1:niter) = errs(nEN,1:niter)/nframe/(L_total-m); % 误帧率 fer(nEN,1:niter) = nferr(nEN,1:niter)/nframe; % 在txt中列出最终的结果 fprintf('************** Eb/N0 = %5.2f db **************\n', EbN0db(nEN)); fprintf('Frame size = %d, rate 1/%d. \n', L_total, 2+puncture); fprintf('%d frames transmitted, %d frames in error.\n', nframe, nferr(nEN, niter)); fprintf('Bit Error Rate (from iteration 1 to iteration %d):\n', niter); for i=1:niter fprintf('%8.4e ', ber(nEN,i)); end fprintf('\n'); fprintf('Frame Error Rate (from iteration 1 to iteration %d):\n', niter); for i=1:niter fprintf('%8.4e ', fer(nEN,i)); end fprintf('\n'); fprintf('***********************************************\n\n'); % 存储中间结果 save turbo_sys_demo EbN0db ber fer end %iter end %while end %nEN diary off