OFDM matlab源程序

1、%main_OFDM.m%一個相對完整的OFDM通信系統(tǒng)的仿真設(shè)計，包括編碼，調(diào)制，IFFT，%上下變頻，高斯信道建模，F(xiàn)FT，PAPR抑制，各種同步，解調(diào)和解碼等模%塊，并統(tǒng)括系統(tǒng)性能的仿真驗證了系統(tǒng)設(shè)計的可靠性。clear allclose allclc%+全局變量+% seq_num 表示當(dāng)前幀是第幾幀% count_dds_up 上變頻處的控制字的累加% count_dds_down 下變頻處的控制字的累加（整整）% count_dds_down_tmp 下變頻處的控制字的累加（小數(shù)）% dingshi 定時同步的定位% m_syn 記錄定時同步中的自相關(guān)平臺global seq_nu

2、mglobal count_dds_upglobal count_dds_downglobal count_dds_down_tmpglobal dingshi global m_syn%+% SNR_Pre 設(shè)定用于仿真的信噪比的初值% interval_SNR 設(shè)定用于仿真的信噪比間隔% frame_num 每一個信噪比下仿真的數(shù)據(jù)幀數(shù)% err_int_final 用于計算每幀出現(xiàn)的誤比特數(shù)% fwc_down 設(shè)定的接收機初始載波頻率控制字% fre_offset 設(shè)定接收機初始載波頻率偏移調(diào)整量（單位為Hz）% k0 每次進入卷積編碼器的信息比特數(shù)% G 卷積編碼的生成矩陣SNR_P

3、re=-5;interval_SNR=1;for SNR_System=SNR_Pre:interval_SNR:5 frame_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 幀數(shù) %+以下為輸入數(shù)據(jù)部分+datain=randint(1,90); %+%+以下為信道卷積編碼部分+encodeDATA=cnv_encd(G

4、,k0,datain);%+%+信道交織編碼+interlacedata=interlacecode(encodeDATA,8,24);%+%+以下為QPSK調(diào)制部分+QPSKdata=qpsk(interlacedata);%+%+生成訓(xùn)練序列+if seq_num<3trainsp_temp=seq_train(); end%+%+插入導(dǎo)頻+PILOT=(1+j);m_QPSKdata=QPSKdata;data2fft_temp=m_QPSKdata(1:8),PILOT,m_QPSKdata(9:16),PILOT,m_QPSKdata(17:24),PILOT,m_QPSKda

5、ta(25:32),PILOT,m_QPSKdata(33:40),PILOT,m_QPSKdata(41:48),m_QPSKdata(49:56),PILOT,m_QPSKdata(57:64),PILOT,m_QPSKdata(65:72),PILOT,m_QPSKdata(73:80),PILOT,m_QPSKdata(81:88),PILOT,m_QPSKdata(89:end);%+trainsp_temp2=trainsp_temp,zeros(1,128);trainsp=trainsp_temp2(65:256),trainsp_temp2(1:64);%+降PAPR矩陣變換

6、+matix_data=nyquistimp_PS();matrix_mult=data2fft_temp*matix_data;%+data2fft2=matrix_mult(65:128),zeros(1,128),matrix_mult(1:64);%+ifft運算+if seq_num=1 ifftin=trainsp;elseif seq_num=2 ifftin=trainsp;else ifftin=data2fft2;endIFFTdata=fft_my(conj(ifftin)/256);IFFTdata=conj(IFFTdata);% figure% plot(real(

7、IFFTdata)% xlabel('realIFFTdata')% figure% plot(imag(IFFTdata)% xlabel('imagIFFTdata') %+%+以下為插入循環(huán)前后綴，2倍升采樣+data2fir=add_CYC_upsample(IFFTdata,2);% +% +fir低通濾波+guiyi_a=0.00172160.0101620.0255120.028801-0.0059219-0.060115-0.04960.0914310.296360.39560.296360.091431-0.0496-0.060115-0.00

8、592190.0288010.0255120.0101620.0017216 ;%抽樣截止頻率為128kHZ,通帶截止頻率為20kHZ，阻帶截止頻率為40kHZ，帶內(nèi)紋波動小于1dB，帶外衰減100dBtxFIRdatai=filter(guiyi_a,1,real(data2fir);txFIRdataq=filter(guiyi_a,1,imag(data2fir);% +%+發(fā)射機cic濾波+CICidatai=cic_inter(txFIRdatai,20);CICidataq=cic_inter(txFIRdataq,20);%+%+上變頻+fwc_up=16; %控制字可以選擇DU

9、Cdata=up_convert_ofdm(fwc_up,CICidatai,CICidataq);%+%+高斯白噪聲信道+DUCdata,datamax=guiyi_DUCdata(DUCdata);awgn_data=awgn(DUCdata,SNR_System); %+%*接受機*%+下變頻+DUCdata_tmp=awgn_data;fwc_down=fwc_down+(fre_offset*128/2560000);r_fre_offset=2560000*(fwc_down-fwc_up)/128);DDCdatai,DDCdataq=down_convert_ofdm(fwc_

10、down,DUCdata_tmp);%+%+接收機cic濾波+CICddatai=cic_deci(DDCdatai,40,40);CICddataq=cic_deci(DDCdataq,40,40);%+%+fir低通濾波+guiyi_b= 0.019527-0.039340.049055-0.018102-0.10030.59440.5944-0.1003-0.0181020.049055-0.039340.019527;%抽樣截止頻率為64kHZ,通帶截止頻率為20kHZ，阻帶截止頻率為30kHZ，帶內(nèi)紋波動小于1dB，帶外衰減60dBrxFIRdatai=filter(guiyi_b,

11、1,CICddatai);rxFIRdataq=filter(guiyi_b,1,CICddataq);%+%+量化+q_rxFIRdatai=sign(rxFIRdatai); q_rxFIRdataq=sign(rxFIRdataq);%+%+定時同步檢測+if seq_num<3time_syn(q_rxFIRdatai,q_rxFIRdataq);end%+%+頻率同步+fre_offset=fre_syn(rxFIRdatai,rxFIRdataq);%+ %+fft運算+if seq_num>2 seq_num-2 fftw=32+dingshi; rxFIRdata_

12、syn=rxFIRdatai(fftw:fftw+255)+j*rxFIRdataq(fftw:fftw+255); FFTdata=fft_my(rxFIRdata_syn); %+%+降PAPR逆矩陣變換+fftdata_reg=FFTdata(193:256),FFTdata(1:64);dematrix_data=fftdata_reg*pinv(matix_data);%+%+相位補償+rx_qpsk_din_th=phase_comp(dematrix_data);%+%+QPSK解調(diào)部分+% figure% plot(rx_qpsk_din_th,'.')% xl

13、abel('星座圖')datatemp4=deqpsk(rx_qpsk_din_th);datatemp4=sign(datatemp4);for m=1:192 if datatemp4(m)=-1 datatemp4(m)=1; elseif datatemp4(m)=1 datatemp4(m)=0; endend%+%+信道解交織+interdout=interlacedecode(datatemp4,8,24);%+%+以下為viterbi譯碼部分+decodeDATA=viterbi(G,k0,interdout);%+%+誤比特統(tǒng)計+err_final=sum(a

14、bs(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*(frame_num-2)disp('-');enddisp('-end-');SNR_System=SNR_Pre:interval_SNR:5;figuresemilogy(SNR_System,err_rate_final,'

15、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 數(shù)，% G 是決定輸入序列的生成矩陣，它有 n0 行 L*k0 列 n0 是輸出 bit 數(shù)，% 參數(shù) n0 和 L 由生成矩陣 G 導(dǎo)出，L 是約束長度。L 之所以叫約束長度% 是因為編碼器在每一時刻里輸

16、出序列不但與當(dāng)前輸入序列有關(guān)，% 而且還與編碼器的狀態(tài)有關(guān)，這個狀態(tài)是由編碼器的前(L-1)k0。% 個輸入決定的,通常卷積碼表示為(n0,k0,m)，m=(L-1)*k0 是編碼% 器中的編碼存貯個數(shù)，也就是分為 L-1 段，每段 k0 個% 有些人將 m=L*k0 定義為約束長度，有的人定義為 m=(L-1)*k0% 查看是否需要補 0，輸入 input 必須是 k0 的整數(shù)部%+variables+% G 決定輸入序列的生成矩陣% k0 每一時鐘周期輸入編碼器的 bit 數(shù)% input 輸入數(shù)據(jù)% output 輸入數(shù)據(jù)%+ if rem(length(input),k0)>0i

17、nput=input,zeros(size(1:k0-rem(length(input),k0);endn=length(input)/k0;% 檢查生成矩陣 G 的維數(shù)是否和 k0 一致if rem(size(G,2),k0)>0error('Error,G is not of the right size.')end% 得到約束長度 L 和輸出比特數(shù) n0L=size(G,2)/k0;n0=size(G,1);% 在信息前后加 0，使存貯器歸 0，加 0 個數(shù)為(L-1)*k0 個u=zeros(size(1:(L-1)*k0),input,zeros(size(1:

18、(L-1)*k0);% 得到 uu 矩陣,它的各列是編碼器各個存貯器在各時鐘周期的內(nèi)容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);% 得到輸出，輸出由生成矩陣 G*uu 得到output=reshape(rem(G*uu,2),1,n0*(L+n-1);% *end of file*%*be

19、ginning of file*%interlacecode.mfunction dout=interlacecode(din,m,n)%實現(xiàn)信道的交織編碼%din為輸入交織編碼器的數(shù)據(jù)，m，n分別為交織器的行列值%+variables+% din 輸入數(shù)據(jù)% m 交織器的行值% n 交織器的列值% dout 輸出數(shù)據(jù)%+ for j=1:m temp(j,:)=din(j*n-(n-1):j*n);enddout_temp=reshape(temp,1,length(din);dout=dout_temp(1:end);%*end of file*%*beginning of file*%q

20、psk.m%QPSK調(diào)制映射function dout=qpsk(din)%+variables+% din 輸入數(shù)據(jù)% dout 輸出數(shù)據(jù)%+ 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*%seq_train.m%生成用于同步的訓(xùn)練符號function dout=seq_train()%第一幀產(chǎn)生短訓(xùn)練序列，第二幀產(chǎn)生長訓(xùn)練序列%每個短訓(xùn)

21、練符號由16個子載波組成，短訓(xùn)練序列%是由偽隨機序列經(jīng)過數(shù)字調(diào)制后插0后，再經(jīng)過%IFFT之后得到的。具體過程如下：首先采用抽頭%系數(shù)為1 0 0 1 的4級移位寄存器產(chǎn)生長度為%15的偽隨機序列之后末尾補0，經(jīng)過QPSK調(diào)制之%后的偽隨機序列只在16的整數(shù)倍位置上出現(xiàn)，其%余的位置補0，產(chǎn)生長度為128的序列，此序列再%補128個0經(jīng)過數(shù)據(jù)搬移后做256點的IFFT變換就%得到16個以16為循環(huán)的訓(xùn)練序列，經(jīng)過加循環(huán)前%后綴就會產(chǎn)生20個相同的短訓(xùn)練序列。長訓(xùn)練序%列的產(chǎn)生同短訓(xùn)練序列。global seq_num if seq_num=1 fbconnection=1 0 0 1; QPS

22、Kdata_pn=m_sequence(fbconnection),0; QPSKdata_pn=qpsk(QPSKdata_pn);elseif seq_num=2 fbconnection=1 0 0 0 0 0 1; QPSKdata_pn=m_sequence(fbconnection),0; QPSKdata_pn=qpsk(QPSKdata_pn);endcountmod=0;for k=1:128 if seq_num=1 if mod(k-1,16)=0 %生成16位循環(huán)的短訓(xùn)練符號 countmod=countmod+1; trainsp_temp(k)=QPSKdata_p

23、n(countmod); else trainsp_temp(k)=0; end elseif seq_num=2 if mod(k-1,2)=0 countmod=countmod+1; trainsp_temp(k)=QPSKdata_pn(countmod); else trainsp_temp(k)=0; end endenddout=trainsp_temp;% *end of file*%*beginning of file*%m_sequence.m%用線性移位寄存器產(chǎn)生m序列function mseq= m_sequence(fbconnection);%+variables+

24、% fbconnection 線性移位寄存器的系數(shù)% mseq 生成的m序列%+ n = length(fbconnection);N = 2n-1;register = zeros(1,n - 1) 1;%定義移位寄存器的初始狀態(tài)mseq(1)= register(n);for i = 2:N newregister(1)= mod(sum(fbconnection.*register),2); for j = 2:n, newregister(j)= register(j-1); end; register = newregister; mseq(i) = register(n);end%

25、 *end of file*%*beginning of file*%nyquistimp_PS.m%使用改進的Nyquist脈沖實現(xiàn)OFDM信號的PAPR抑制function dout=nyquistimp_PS()%改進的Nyquist脈沖整形方法能夠顯著改善OFDM信%號的PAPR分布；該方法實現(xiàn)簡單，和PTS和SLM相比%不需迭代計算多個IFFT操作，不需傳送邊帶信息，%不會引起信號的畸變；通用性強，可以調(diào)整滾降%系數(shù)以適應(yīng)任何子載波數(shù)的通信系統(tǒng)。當(dāng)然，%Nyquist脈沖成形的方法由于擴展了頻譜，一定程% 度上降低了頻譜利用率。%creat a matrix to shape the

26、 subcarries.%the spectrum of the pulse is as follows:% if abs(f)<=Bw*(1-b)% spec=1;% else if (abs(f)>Bw*(1-b)&(abs(f)<=Bw) % spec=exp(aa.*(Bw.*(1-b)-abs(f);% else if (abs(f)>Bw)&(abs(f)<Bw*(1+b) % spec=1-exp(aa.*(abs(f)-Bw.*(1+b);% else if abs(f)>=Bw*(1+b)% spec=0; % end% e

27、nd% 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;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:finish kk=kk+1; if abs(f)<=Bw*(1-b) spec=1; else if (a

28、bs(f)>Bw*(1-b)&(abs(f)<=Bw) spec=exp(aa.*(Bw.*(1-b)-abs(f); else if (abs(f)>Bw)&(abs(f)<Bw*(1+b) spec=1-exp(aa.*(abs(f)-Bw.*(1+b); else if abs(f)>=Bw*(1+b) spec=0; end end end end C(kk)=spec;endfor m=0:N-1 for k=0:(N+2*L-1) p(m+1,k+1)=C(k+1)*exp(-i*2*pi.*m.*(k-L)./N); endenddo

29、ut=p;% *end of file*%*beginning of file*%fft_my.m%實現(xiàn)N點FFT運算function dout=fft_my(din)%本程序?qū)斎胄蛄衐in實現(xiàn)DITFFT基2算法，點數(shù)取大于等于din長度的2的冪次%+variables+% din 輸入數(shù)據(jù)% dout 輸出數(shù)據(jù)%+ m=nextpow2(din) ;N=2m ;if length(din)<N din=din,zeros(1,N-length(din);endnxd=bin2dec(fliplr(dec2bin(1:N-1,m)+1;y=din(nxd);for mm=1:mNmr

30、=2mm;u=1;WN=exp(-i*2*pi/Nmr); for j=1:Nmr/2 for k=j:Nmr:N kp=k+Nmr/2; t=y(kp)*u; y(kp)=y(k)-t; y(k)=y(k)+t; end u=u*WN; endenddout=y;% *end of file*%*beginning of file*%add_GI_upsample.m%加循環(huán)前后綴和升采樣程序function dout=add_CYC_upsample(din,upsample)%插入循環(huán)前后綴是將每個OFDM符號的前32個數(shù)據(jù)放%到符號尾部，將每個OFDM符號的后32個數(shù)據(jù)放到符號頭部，%

31、升采樣是通過中間插零的方式實現(xiàn)%+variables+% din 輸入數(shù)據(jù)% upsample 升采樣倍數(shù) % dout 輸出數(shù)據(jù)%+ data_buf=din(225:256),din,din(1:32);temp(1,:)=data_buf;temp(2:upsample,:)=zeros(upsample-1,length(data_buf);dout_temp=reshape(temp,1,length(data_buf)*upsample);dout=dout_temp(1:end);% *end of file*%*beginning of file*%cic_inter.m%發(fā)射

32、機的CIC濾波器設(shè)計function dout=cic_inter(din,r)%插值CIC濾波器通過復(fù)制樣值實現(xiàn)升采樣%+variables+% din 輸入數(shù)據(jù)% r 升采樣的插值因子% dout 輸出數(shù)據(jù)%+ din_1d=0,din;diff1=din,0-din_1d;diff1_1d=0,diff1(1:end-1);diff2=diff1(1:end-1),0-diff1_1d;for i=1:r temp1(i,:)=diff2(1:end-1);endtemp2=reshape(temp1,1,length(diff2(1:end-1)*r);data_int=temp2;f

33、or i=1:length(data_int) int1(i)=sum(data_int(1:i); int2(i)=sum(int1(1:i);enddout=int2;%*end of file*%*beginning of file*%up_convert_ofdm.m%用DDS的方式實現(xiàn)上變頻function dout=up_convert_ofdm(fwc_up,dini,dinq)%fwc_up是上變頻處的頻率控制字，每個控制字對應(yīng)一個載波頻率%count_dds_up是用于查找存儲表的整數(shù)值，%+variables+% fwc_down 上變頻處的頻率控制字% dini 輸入數(shù)據(jù)

34、的實部% dinq 輸入數(shù)據(jù)的虛部% dout 輸出數(shù)據(jù)%+ global seq_numglobal count_dds_upfor mk=1:length(dini) if (seq_num=1) & (mk=1) count_dds_up=0; else count_dds_up=count_dds_up+fwc_up; if count_dds_up>=128 count_dds_up=count_dds_up-128; end end up_sin,up_cos=ram_sin(count_dds_up); up_sin_t(mk)=up_sin; up_cos_t(mk)=up_cos;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%構(gòu)造用于DDS的查找表function ysin,ycos=ram_sin(adr)%dds方式需要的sin表%ram_sin為寄存器名稱%adr為輸入地址,y為讀出數(shù)據(jù)%+variables+% adr 輸入地址% ysin 輸出的正弦值% ycos 輸出的余