OFDM_matlab源程序.doc

1、%main_OFDM.m%一个相对完整的 OFDM 通信系统的仿真设计，包括编码，调制，IFFT，%上下变频，高斯信道建模，FFT，PAPR 抑制，各种同步，解调和解码等模%块，并统括系统性能的仿真验证了系统设计的可靠性。clear allclose allclc%+全局变量+% seq_num 表示当前帧是第几帧% count_dds_up 上变频处的控制字的累加% count_dds_down 下变频处的控制字的累加（整整）% count_dds_down_tmp 下变频处的控制字的累加（小数）% dingshi 定时同步的定位% m_syn 记录定时同步中的自相关平台global seq

2、_numglobal count_dds_upglobal count_dds_downglobal count_dds_down_tmpglobal dingshi global m_syn%+% SNR_Pre 设定用于仿真的信噪比的初值% interval_SNR 设定用于仿真的信噪比间隔% frame_num 每一个信噪比下仿真的数据帧数% err_int_final 用于计算每帧出现的误比特数% fwc_down 设定的接收机初始载波频率控制字% fre_offset 设定接收机初始载波频率偏移调整量（单位为 Hz）% k0 每次进入卷积编码器的信息比特数% G 卷积编码的生成矩阵S

3、NR_Pre=-5;interval_SNR=1;for SNR_System=SNR_Pre:interval_SNR:5frame_num=152;dingshi=250;err_int_final=0;fwc_down=16.050;fre_offset=0;k0=1;G=1 0 1 1 0 1 1;1 1 1 1 0 0 1 ;disp(-start-);for seq_num=1:frame_num, %frame_num 帧数%+以下为输入数据部分+datain=randint(1,90); %+%+以下为信道卷积编码部分+encodeDATA=cnv_encd(G,k0,data

4、in);%+%+信道交织编码+interlacedata=interlacecode(encodeDATA,8,24);%+%+以下为 QPSK 调制部分 +QPSKdata=qpsk(interlacedata);%+%+生成训练序列+if seq_num2seq_num-2fftw=32+dingshi;rxFIRdata_syn=rxFIRdatai(fftw:fftw+255)+j*rxFIRdataq(fftw:fftw+255);FFTdata=fft_my(rxFIRdata_syn);%+%+降 PAPR 逆矩阵变换+fftdata_reg=FFTdata(193:256),F

5、FTdata(1:64);dematrix_data=fftdata_reg*pinv(matix_data);%+%+相位补偿+rx_qpsk_din_th=phase_comp(dematrix_data);%+%+QPSK 解调部分+% figure% plot(rx_qpsk_din_th,.)% xlabel(星座图)datatemp4=deqpsk(rx_qpsk_din_th);datatemp4=sign(datatemp4);for m=1:192if datatemp4(m)=-1datatemp4(m)=1;elseif datatemp4(m)=1datatemp4(m

6、)=0;endend%+%+信道解交织+interdout=interlacedecode(datatemp4,8,24);%+%+以下为 viterbi 译码部分+decodeDATA=viterbi(G,k0,interdout);%+%+误比特统计+err_final=sum(abs(decodeDATA-datain) err_int_final=err_int_final+err_final endenddisp(-);SNR_Systemerr_rate_final(SNR_System-SNR_Pre)./interval_SNR+1)=err_int_final/(90*(fr

7、ame_num-2)disp(-);enddisp(-end-);SNR_System=SNR_Pre:interval_SNR:5;figuresemilogy(SNR_System,err_rate_final,b-*);xlabel(信噪比/dB)ylabel(误码率)axis(-5,5,0,1)grid on%+%*beginning of file*%cnv_encd.m%卷积码编码程序function output=cnv_encd(G,k0,input)% cnv_encd(G,k0,input),k0 是每一时钟周期输入编码器的 bit 数，% G 是决定输入序列的生成矩阵，它

8、有 n0 行 L*k0 列 n0 是输出 bit 数，% 参数 n0 和 L 由生成矩阵 G 导出，L 是约束长度。 L 之所以叫约束长度% 是因为编码器在每一时刻里输出序列不但与当前输入序列有关，% 而且还与编码器的状态有关，这个状态是由编码器的前(L-1)k0。% 个输入决定的,通常卷积码表示为 (n0,k0,m)，m=(L-1)*k0 是编码% 器中的编码存贮个数，也就是分为 L-1 段，每段 k0 个% 有些人将 m=L*k0 定义为约束长度，有的人定义为 m=(L-1)*k0% 查看是否需要补 0，输入 input 必须是 k0 的整数部%+variables+% G 决定输入序列的

9、生成矩阵% k0 每一时钟周期输入编码器的 bit 数% input 输入数据% output 输入数据%+ if rem(length(input),k0)0input=input,zeros(size(1:k0-rem(length(input),k0);endn=length(input)/k0;% 检查生成矩阵 G 的维数是否和 k0 一致if rem(size(G,2),k0)0error(Error,G is not of the right size.)end% 得到约束长度 L 和输出比特数 n0L=size(G,2)/k0;n0=size(G,1);% 在信息前后加 0，使存

10、贮器归 0，加 0 个数为(L-1)*k0 个u=zeros(size(1:(L-1)*k0),input,zeros(size(1:(L-1)*k0);% 得到 uu 矩阵,它的各列是编码器各个存贮器在各时钟周期的内容u1=u(L*k0:-1:1);%将加 0 后的输入序列按每组 L*k0 个分组，分组是按 k0 比特增加%从 1 到 L*k0 比特为第一组，从 1+k0 到 L*k0+k0 为第二组， 。 。 。 。 ，%并将分组按倒序排列。for i=1:n+L-2u1=u1,u(i+L)*k0:-1:i*k0+1);enduu=reshape(u1,L*k0,n+L-1);% 得到输出

11、，输出由生成矩阵 G*uu 得到output=reshape(rem(G*uu,2),1,n0*(L+n-1);% *end of file*%*beginning of file*%interlacecode.mfunction dout=interlacecode(din,m,n)%实现信道的交织编码%din 为输入交织编码器的数据，m，n 分别为交织器的行列值%+variables+% din 输入数据% m 交织器的行值% n 交织器的列值% dout 输出数据%+ for j=1:mtemp(j,:)=din(j*n-(n-1):j*n);enddout_temp=reshape(t

12、emp,1,length(din);dout=dout_temp(1:end);%*end of file*%*beginning of file*%qpsk.m%QPSK 调制映射function dout=qpsk(din)%+variables+% din 输入数据% dout 输出数据%+ din2=1-2*din;din_temp=reshape(din2,2,length(din)/2);for i=1:length(din)/2,dout(i)=din_temp(1,i)+j*din_temp(2,i);end% *end of file*%*beginning of file*

13、%seq_train.m%生成用于同步的训练符号function dout=seq_train()%第一帧产生短训练序列，第二帧产生长训练序列%每个短训练符号由 16 个子载波组成，短训练序列%是由伪随机序列经过数字调制后插 0 后，再经过%IFFT 之后得到的。具体过程如下：首先采用抽头%系数为1 0 0 1 的 4 级移位寄存器产生长度为%15 的伪随机序列之后末尾补 0，经过 QPSK 调制之%后的伪随机序列只在 16 的整数倍位置上出现，其%余的位置补 0，产生长度为 128 的序列，此序列再%补 128 个 0 经过数据搬移后做 256 点的 IFFT 变换就%得到 16 个以 16

14、 为循环的训练序列，经过加循环前%后缀就会产生 20 个相同的短训练序列。长训练序%列的产生同短训练序列。global seq_num if seq_num=1fbconnection=1 0 0 1;QPSKdata_pn=m_sequence(fbconnection),0;QPSKdata_pn=qpsk(QPSKdata_pn);elseif seq_num=2fbconnection=1 0 0 0 0 0 1;QPSKdata_pn=m_sequence(fbconnection),0;QPSKdata_pn=qpsk(QPSKdata_pn);endcountmod=0;for

15、k=1:128 if seq_num=1 if mod(k-1,16)=0 %生成 16 位循环的短训练符号 countmod=countmod+1;trainsp_temp(k)=QPSKdata_pn(countmod);elsetrainsp_temp(k)=0;endelseif seq_num=2 if mod(k-1,2)=0countmod=countmod+1;trainsp_temp(k)=QPSKdata_pn(countmod);elsetrainsp_temp(k)=0;endendenddout=trainsp_temp;% *end of file*%*beginn

16、ing of file*%m_sequence.m%用线性移位寄存器产生 m 序列function mseq= m_sequence(fbconnection);%+variables+% fbconnection 线性移位寄存器的系数% mseq 生成的 m 序列%+ n = length(fbconnection);N = 2n-1;register = zeros(1,n - 1) 1;%定义移位寄存器的初始状态mseq(1)= register(n);for i = 2:Nnewregister(1)= mod(sum(fbconnection.*register),2);for j

17、= 2:n,newregister(j)= register(j-1);end;register = newregister;mseq(i) = register(n);end% *end of file*%*beginning of file*%nyquistimp_PS.m%使用改进的 Nyquist 脉冲实现 OFDM 信号的 PAPR 抑制function dout=nyquistimp_PS()%改进的 Nyquist 脉冲整形方法能够显著改善 OFDM 信%号的 PAPR 分布；该方法实现简单，和 PTS 和 SLM 相比%不需迭代计算多个 IFFT 操作，不需传送边带信息，%不会

18、引起信号的畸变；通用性强，可以调整滚降%系数以适应任何子载波数的通信系统。当然，%Nyquist 脉冲成形的方法由于扩展了频谱，一定程% 度上降低了频谱利用率。%creat a matrix to shape the subcarries.%the spectrum of the pulse is as follows:% if abs(f)Bw*(1-b) % end% end% end% end N=106;L=11;b=0.22;% N=84;% L=22;% b=0.5;% N=98;% L=15;% b=0.3;% N=116;% L=6;% b=0.1;T=0.004;Ts=T/N

19、;Bw=1/Ts;begin=-Bw*(1+b)+Bw*(1+b)/128;finish=Bw*(1+b)-Bw*(1+b)/128;distance=Bw*(1+b)/64;kk=0;aa=log(2)/(b.*Bw);for f=begin:distance:finishkk=kk+1;if abs(f)Bw*(1-b) endendendend C(kk)=spec;endfor m=0:N-1for k=0:(N+2*L-1)p(m+1,k+1)=C(k+1)*exp(-i*2*pi.*m.*(k-L)./N);endenddout=p;% *end of file*%*beginni

20、ng of file*%fft_my.m%实现 N 点 FFT 运算function dout=fft_my(din)%本程序对输入序列 din 实现 DITFFT 基 2 算法，点数取大于等于 din 长度的 2 的幂次%+variables+% din 输入数据% dout 输出数据%+ m=nextpow2(din) ;N=2m ;if length(din)=128count_dds_up=count_dds_up-128;endendup_sin,up_cos=ram_sin(count_dds_up);up_sin_t(mk)=up_sin;up_cos_t(mk)=up_cos;

21、endfor xl=1:length(dini)DUCdata(xl)=dini(xl)*up_cos_t(xl)-dinq(xl)*up_sin_t(xl);enddout=DUCdata;% *end of file*%*beginning of file*%ram_sin.m%构造用于 DDS 的查找表function ysin,ycos=ram_sin(adr)%dds 方式需要的 sin 表%ram_sin 为寄存器名称%adr 为输入地址,y 为读出数据%+variables+% adr 输入地址% ysin 输出的正弦值% ycos 输出的余弦值%+ anl_inc=2*pi/1

22、28;for n=1:128data_sin(n)=sin(n-1)*anl_inc);data_cos(n)=cos(n-1)*anl_inc);endysin=data_sin(adr+1);ycos=data_cos(adr+1);% *end of file*%*beginning of file*%guiyi_DUCdata.mfunction dataout,datamax=guiyi_DUCdata(datain)%实现数据的归一化%+variables+% datain 输入数据% dataout 输出数据% datamax 输入数据中的最大值%+ datamax=max(ab

23、s(datain);dataout=datain./datamax;% *end of file*%*beginning of file*%down_convert_ofdm.mfunction douti,doutq=down_convert_ofdm(fwc_down,din)%用 DDS 的方式实现下变频%+variables+% fwc_down 下变频处的频率控制字% din 输入数据% douti 输出数据的实部% doutq 输出数据的虚部%+ global seq_num global count_dds_downglobal count_dds_down_tmpfor mkd

24、=1:length(din)if (seq_num=1) count_dds_down_tmp=0;else count_dds_down=round(count_dds_down_tmp+fwc_down);count_dds_down_tmp=count_dds_down_tmp+fwc_down;if count_dds_down=128count_dds_down=count_dds_down-128;count_dds_down_tmp=count_dds_down_tmp-128;endendup_sin_d,up_cos_d=ram_sin(count_dds_down);up_

25、sin_td(mkd)=up_sin_d;up_cos_td(mkd)=up_cos_d;DDCdatai(mkd)=din(mkd)*up_cos_td(mkd);DDCdataq(mkd)=-din(mkd)*up_sin_td(mkd);enddouti=DDCdatai;doutq=DDCdataq;% *end of file*%*beginning of file*%cic_deci.m%接收机的 CIC 滤波器设计function dout=cic_deci(din,r,init)%抽取 CIC 滤波器通过降采样实现%+variables+% din 输入数据% r 降采样的抽取

26、因子% init 设定的初始值% dout 输出数据%+ for i=1:length(din),int1(i)=sum(din(1:i);int2(i)=sum(int1(1:i);enddata_diff=zeros(1,length(int2)/r);data_diff=int2(init:r:end);data_1d=data_diff,0;diff1=0,data_diff-data_1d;diff1_1d=diff1,0;diff2=0,diff1-diff1_1d;dout=diff2(1:end-2);% *end of file*%*beginning of file*%ti

27、me_syn.m%系统的定时同步function time_syn(datai,dataq)%通过前导结构的两个训练帧的延时自相关算法和本地%互相关检测可以实现精确度非常高的定时同步。global seq_num global dingshi global m_syn if seq_num=1for nc=1:length(datai)-64 %计算相关值 for m=1:32m1_syn(m)=(datai(nc+m-1)+j*dataq(nc+m-1)*conj(datai(nc+m-1+16)+j*dataq(nc+m-1+16);endm2_syn(nc)=sum(m1_syn);m_

28、syn(nc)=abs(m2_syn(nc); %自相关自相关判决函数end% figure% plot(m_syn)% xlabel(采样点索引号)% ylabel(自相关判决函数)elseif seq_num=2local_seq= -1.0000 + 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i 1.0000 -

29、 1.0000i 1.0000 - 1.0000i 1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.

30、0000i 1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i -1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000

31、i 1.0000 + 1.0000i 1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i 1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i -1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1

32、.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i 1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 - 1.0000i 1.00

33、00 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 +

34、 1.0000i -1.0000 + 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 + 1.0000i 1.0000 + 1.0000i 1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.00

35、00 - 1.0000i -1.0000 - 1.0000i 1.0000 - 1.0000i 1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i -1.0000 + 1.0000i ;for nn=1:length(datai)-128 %计算相关值 (输

36、入信号与本地信号互相关)for m=1:100t1_syn(m)=(datai(nn+m-1)+j*dataq(nn+m-1)*conj(local_seq(m);endlolol(nn)=sum(t1_syn);t_syn(nn)=abs(sum(t1_syn); %输入信号与本地信号互相关判决函数endfor ni=1:length(t_syn)if t_syn(ni)60dingshi=find(t_syn(ni:ni+6)=max(t_syn(ni:ni+6)+ni-1;breakendend% figure% plot(t_syn) % xlabel(采样点索引号)% ylabel(互相关判决函数)end% *end of file*%*beginning of file*