MATLAB OFDM卷积编码程序及代码

%bin22deci.m

function y=bin22deci(x)

%将二进制数转化为十进制数

t=size(x,2);

y=(t-1:-1:0);

y=2.^y;

y=x*y';

%************************end of file***********************************

%comb.m

%AWGN加噪声程序

function[iout,qout]=comb(idata,qdata,attn)

%******************variables*************************

%idata:输入I信道数据

%qdata:输入Q信道数据

%iout输出I信道数据

%qout输出Q信道数据

%attn:由信噪比导致的衰减系数

%******************************************************

iout=randn(1,length(idata)).*attn;

qout=randn(1,length(qdata)).*attn;

iout=iout+idata(1:length(idata));

qout=qout+qdata(1:length(qdata));

%************************end of file***********************************

%crdemapping.m

%数据逆映射载波程序

function[iout,qout]=crdemapping(idata,qdata,fftlen,nd);

%******************variables*************************

%idata:输入I信道的数据

%qdata:输入Q信道的数据

%iout:输出I信道的数据

%qout:输出Q信道的数据

%fftlen:FFT的长度

%nd:OFDM符号数

%*****************************************************

iout(1:26,:)=idata(2:27,:);

qout(1:26,:)=qdata(2:27,:);

iout(27:52,:)=idata(39:64,:);

qout(27:52,:)=qdata(39:64,:);

%********************end of file***************************

%crmapping.m

%数据映射载波程序

function[iout,qout]=crmapping(idata,qdata,fftlen,nd);

%******************variables*************************

%idata:输入I信道的数据

%qdata:输入Q信道的数据

%iout:输出I信道的数据

%qout:输出Q信道的数据

%fftlen:FFT的长度

%nd:OFDM符号数

%*****************************************************

iout=zeros(fftlen,nd);

qout=zeros(fftlen,nd);

iout(2:27,:)=idata(1:26,:);

qout(2:27,:)=qdata(1:26,:);

iout(39:64,:)=idata(27:52,:);

qout(39:64,:)=qdata(27:52,:);

%********************end of file***************************

%deci22bin.m

function y=deci22bin(x,t)

%十进制数x转化为二进制数，二进制数至少表示为t位

y=zeros(size(x,1),t);

for j=1:size(x,1)

i=1;

while x(j)>=0&i<=t

y(j,i)=rem(x(j),2);%x(j)为偶数时,y(j,i)为0;反之为1

x(j)=(x(j)-y(j,i))/2;

i=i+1;

end

y(j,:)=y(j,t:-1:1);%倒序排列

end

%************************end of file***********************************

%giins1.m

%插入保护间隔程序

function[iout,qout]=giins1(idata,qdata,fftlen,gilen,nd);

%******************变量*************************

%idata:输入I信道数据

%qdata:输入Q信道数据

%iout:输出I信道数据

%qout:输出Q信道数据

%fftlen:FFT长度(points)

%gilen:保护间隔长度(points)

%*****************************************************

idata1=reshape(idata,fftlen,nd);

qdata1=reshape(qdata,fftlen,nd);

idata2=[idata1(fftlen-gilen+1:fftlen,:);idata1];

qdata2=[qdata1(fftlen-gilen+1:fftlen,:);qdata1];

iout=reshape(idata2,1,(fftlen+gilen)*nd);

qout=reshape(qdata2,1,(fftlen+gilen)*nd);

%********************end of file***************************

%qpskdemod1.m

%QPSK解调程序

function[demodata]=qpskdemod1(idata,qdata,para,nd,ml)

%******************variables*************************

%idata:输入I信道数据

%qdata:数据Q信道数据

%demodata:解调后数据(para-by-nd matrix)

%para:并行信道数

%nd:符号数

%ml:调制数

%(QPSK->2 16QAM->4)

%*****************************************************

demodata=zeros(para,ml*nd);

demodata((1:para),(1:ml:ml*nd-1))=idata((1:para),(1:nd))>=0; demodata((1:para),(2:ml:ml*nd))=qdata((1:para),(1:nd))>=0;

%************************end of file***********************************

%qpskmod1.m

%QPSK调制程序

function[iout,qout]=qpskmod1(paradata,para,nd,ml)

%******************variables*************************

%paradata:输入数据

%iout:输出数据I

%qout:输出数据Q

%para:并行信道数

%nd:数据数

%ml:调制数

%(QPSK->2 16QAM->4)

%*****************************************************

m2=ml./2;

paradata2=paradata.*2-1;

count2=0;

for jj=1:nd

isi=zeros(para,1);

isq=zeros(para,1);

isi=isi+paradata2((1:para),1+count2);

isq=isq+paradata2((1:para),2+count2);

iout((1:para),jj)=isi;

qout((1:para),jj)=isq;

count2=count2+ml;

end

%********************end of file***************************

%viterbi.m

%viterbi解码程序

function[decoder_output,survivor_state,cumulated_metric]=viterbi(G,k,channel_output) %[decoder_output,survivor_state,cumulated_metric]=viterbi(G,k,channel_output)

%其中G是一个n行L*k列矩阵，它的每一行决定了从移位寄存器到输入码字的连接方式.

%survivor_state是一个矩阵，它显示了通过网格的最优路径，这个矩阵通过一个单独

%的函数metric(x,y)给出。

n=size(G,1);

%检验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('channle output not of the right size')

end

L=size(G,2)/k;

number_of_states=2^((L-1)*k);

%产生状态转移矩阵，输出矩阵和输入矩阵

for j=0:number_of_states-1

for t=0:2^k-1

[next_state,memory_contents]=nxt_stat(j,t,L,k);

input(j+1,next_state+1)=t;

branch_output=rem(memory_contents*G',2);

nextstate(j+1,t+1)=next_state;

output(j+1,t+1)=bin2deci(branch_output);

end

end

input;

state_metric=zeros(number_of_states,2);

depth_of_trellis=length(channel_output)/n;

channel_output_matrix=reshape(channel_output,n,depth_of_trellis);

survivor_state=zeros(number_of_states,depth_of_trellis+1);

[row_survivor col_survivor]=size(survivor_state);

%开始非尾信道输出的解码

%i为段，j为每一阶段的状态，t为输入

for i=1:depth_of_trellis-L+1

flag=zeros(1,number_of_states);

if i<=L

step=2^((L-i)*k);

else

step=1;

end

for j=0:step:number_of_states-1

for t=0:2^k-1

branch_metric=0;

binary_output=deci2bin(output(j+1,t+1),n);

for tt=1:n

branch_metric=branch_metric+metric(channel_output_matrix(tt,i),binary_output(tt));

end

if

((state_metric(nextstate(j+1,t+1)+1,2)>state_metric(j+1,1)+branch_metric)|flag(nextstate(j+1 ,t+1)+1)==0)

state_metric(nextstate(j+1,t+1)+1,2)=state_metric(j+1,1)+branch_metric;

survivor_state(nextstate(j+1,t+1)+1,i+1)=j;

flag(nextstate(j+1,t+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 tt=1:n

branch_metric=branch_metric+metric(channel_output_matrix(tt,i),binary_output(tt));

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

%从最优路径产生解码输出

%由段得到状态序列，再由状序列从input矩阵中得到该段的输出

state_sequence=zeros(1,depth_of_trellis+1);

size(state_sequence);

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

state_sequence;

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));

cumulated_metric=state_metric(1,1);

%************************end of file***********************************

%ofdm1.m

%QPSK调制、AWGN信道下仿真程序

%**********************初始参数***************************

para=52;%并行信道数

fftlen=64;%FFT长度

noc=52;%载波数

nd=6;%每循环中OFDM符号数

ml=2;%调制水平:QPSK

gilen=16;%保护间隔长度(points)

ebn0=2;%信噪比

sr=250000;%OFDM symbol rate(250 ksyombol/s)

br=sr.*ml;%Bit rate per carrier

%**************************主循环部分**************************

nloop=200;%仿真的循环数

noe1=0;%信道解码前错误数据数

nod1=0;%信道解码前传输数据数

noe2=0;%信道解码后错误数据数

nod2=0;%信道解码后传输数据数

%**************************发射机*********************************

for iii=1:nloop

%**************************数据产生**************************** %信源编码

t=[0:pi/25:2*pi];

xx=sin(t);

init=[-1:.1:1];

partition=[-1:.1:.9];

predictor=[0 1];

encode=dpcmenco(xx,init,partition,predictor);

encode2=reshape(encode,51,1);

bin=deci22bin(encode2,6);

recode=reshape(bin,1,306);

%信道编码

k0=1;

G=[1 0 1 1 0 1 1;1 1 1 1 0 0 1];

channelencode=cnv_encd(G,k0,recode);

%******************串并转换***********************

paradata=reshape(channelencode,para,nd*ml);%reshape:内建功能

%**************************QPSK调制***************************** [ich,qch]=qpskmod1(paradata,para,nd,ml);

kmod=1/sqrt(2);%sqrt:内建功能

ich=ich.*kmod;

qch=qch.*kmod;

%数据映射

[ich1,qch1]=crmapping(ich,qch,fftlen,nd);

%*******************IFFT************************

x=ich1+qch1.*i;

y=ifft(x);%ifft:内建功能

ich2=real(y);%real:内建功能(实部)

qch2=imag(y);%imag:内建功能（虚部）

%*********插入保护间隔**********

[ich3,qch3]=giins1(ich2,qch2,fftlen,gilen,nd);

fftlen2=fftlen+gilen;

%------------------衰减计算-----------------------

spow=sum(ich3.^2+qch3.^2)/nd./para;%sum:内建功能

attn=0.5*spow*sr/br*10.^(-ebn0/10);

attn=sqrt(attn);

%---------------AWGN addition---------------

[ich4,qch4]=comb(ich3,qch3,attn);

%******************去除保护间隔*********

[ich5,qch5]=girem1(ich4,qch4,fftlen2,gilen,nd);

%******************FFT******************

rx=ich5+qch5.*i;

ry=fft(rx);

ich6=real(ry);

qch6=imag(ry);

%载波逆映射

[ich7,qch7]=crdemapping(ich6,qch6,fftlen,nd);

%*****************解调*******************

ich7=ich7./kmod;

qch7=qch7./kmod;

[demodata]=qpskdemod1(ich7,qch7,para,nd,ml);

%**************并串转换*****************

demodata1=reshape(demodata,1,para*nd*ml);

%信道解码

[channeldecode]=viterbi(G,k0,demodata1);

%信源解码

reshapechanneldecode=reshape(channeldecode,51,6);

deci=bin22deci(reshapechanneldecode);

redeci=reshape(deci,1,51);

codebook=[-1:.1:5.3];

decode=dpcmdeco(redeci,codebook,predictor);

%**************************比特误码率(BER)**************************** noe10=sum(abs(demodata1-channelencode));

nod10=length(channelencode);

noe20=sum(abs(channeldecode-recode));

nod20=length(recode);

noe1=noe10+noe1;

nod1=nod10+nod1;

noe2=noe20+noe2;

nod2=nod20+nod2;

fprintf('%d\t%e\t%e\n',iii,noe10/nod10,noe20/nod20);

end

ber1=noe1/nod1;

ber2=noe2/nod2;

%**********************输出结果***************************

fprintf('%f\t%e\t%e\t%d\t\n',ebn0,ber1,ber2,nloop);

%************************end of file***********************************