ofdmsystemsframeworkprogramme120408.m

OFDM 框架程序 修改日期2008-12-04 仿真环境 Matlab6.5

%%%%%%%%%%%%%%    OFDM 框架程序   修改日期2008-12-04         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 仿真环境 Matlab6.5
% Author 。
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%barthannwin hamming   bartlett hann blackman kaiser blackmanharris nuttallwin    各种窗函数 
% bohmanwin   parzenwin chebwin  rectwin  flattopwin  triang  gausswin  tukeywin  
%n=64;wvtool(hamming(n),hann(n),nuttallwin(n),blackmanharris(n),barthannwin(n));%%%显示各种窗函数
% function();
% starttime=['starttime' date '.fig'];
% hgsave(['E:\BUPT paper and matlab programme\OFDM systems framework programme\OFDM systems framework programme result\' starttime ]);
% 确定仿真开始的时间。在设定的文件夹下查看被保存文件“starttime05-May-2008”的属性，例如可以看到：创建时间: 2008年5月5日, 11:16:17。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%if exist('1time_iffttest_030-Oct-2003.mat');
   %load ('1time_iffttest_030-Oct-2003.mat');A=1; 
   %elseif
    %A=2;
    %else
    %A=3;
    %end
%可以将已经存在的文件调入，利用调入的文件中的已有数据
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;% clear workspace      清空MATLAB的工作区
close all;% close all the figure 关闭所有的图形
clc;      % clear command window 清空命令窗口
'If programme is pausing, knock any key to activate'
%pause(10) 暂停10秒
%pause      %暂停直到有键盘输入
%-------对MATLB的工作环境进行前期处理--------------------------------------------------
tic % tic starts a stopwatch timer. toc prints the elapsed time since tic was used. t = toc returns the elapsed time in t.
global N; %define the global variable.N is the number point of ifft or fft.
N=1024;
M=2;     % M-QAM.
FrameNumber=100;
% ------定义了fft和ifft变换的点数 N ; 仿真的帧数 FrameNumber ；多进制（M）数字调制.-------------
finaldata=zeros(N,FrameNumber); %Create an array of all zeros，default value is 'double'。用于存放接收端逆映射（解调）后的信息BITS.
Original_X=randint(N,FrameNumber,M);%产生串并变换后的信息BITS. generate original  data of N line*FrameNumber row random int matrix.
%The entries are uniformly distributed and independently chosen from the range [0,M-1]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Xk1=modmap(Original_X,1,1,'qam',M);
[Original_X_Real,Original_X_Imaginary]=qaskenco(Original_X,M);
Xk= Original_X_Real + Original_X_Imaginary*j;%映射（调制），产生频域信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
modmap('qam',M,'*b')%%% 显示MQAM的星座图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)        %%% 画出MQAM星座图
title('QAM constellation ')
text(0.1,0.2,'QPSK constellation')
plot(sqrt(M),sqrt(M),'*')
hold on
plot(-sqrt(M),-sqrt(M),'*')
hold on
plot(Original_X_Real,Original_X_Imaginary,'*r')
xlabel('Re')
ylabel('Im')
grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear Original_X_Real;clear Original_X_Imaginary;%注意及时清除不在使用的矩阵，节省内存，提高运算速度。
%-------产生了一个N*FrameNuMber的复数矩阵 X(K)，其中每一个元素都是由M-QAM映射来的.--------------
xn=ifft(Xk,N);% 对 X(k) 进行N点的离散傅立叶反变换得到 x(n).在这里与xn=ifft(Xk)等价，因为Xk是一个N行的矩阵。产生IFFT后的时域信号
                    %//=fft(xn);
CP=6;
xn_add_cp=[xn(N-CP+1:N,:);xn(:,:)];% add CP to x(n)。CP使得OFDM传输在一定条件下可以完全消除ISI和ICI的影响。
%-------将x(n)的后 N-CP+1 到 N 的元素平移到 x(n)的前面.与xnaddcp=[xn(N-CP+1:N,:);xn]是一样的结果。
%//h=randint(6,1,10001)/10000;%产生6×1的多径信道各径的衰落因子矩阵h(t)，每一个元素是0－1的随机数。
                    %在数字通信第四版有这方面的介绍。何波比较熟悉。
h = [1,0.6,0.4,0.16,0.04,0.01]'; % 'denote 共轭转置.或者h＝zeros[6,1];h[1,1]=1....先产生6×1的零矩阵，再给每个元素赋值。
h = [1,0.6,0.4,0.16,0.04,0.01].';% .' denotes 转置.
% h = [0.6+0.6i;zeros(6,1);0.6;0;0;0.5i;zeros(10,1);0.2-0.1i;zeros(8,1);0.3;0].'; %深衰落
% h = [1+i;zeros(6,1);0.5-0.3i;0;0.3+0.2i;0;zeros(10,1);0.2-0.1i;zeros(8,1);0.1;0].';
% 产生多径瑞利衰落信道（Rayleigh fading channel）。即，信道响应的能量或包络服从瑞利分布，均值为0，且相位服从0 到2π 的均匀分布。
% 通常将信道增益以等效基带信号表示，即用一复数表示信道的幅度和相位特性。瑞利衰落可由这一复数表示，它的实部和虚部服从于零均值的独立同分布高斯过程。
h_shield = [1,1,1,1,1,1];% 对应的O位置，说明该径为0.
h_power  = [1,0.6,0.4,0.16,0.04,0.01]; % randn() pseudorandom,scalar value drawn from a normal distribution with mean 0 and standard deviation 1.
h = randn(1,length(h)).*h_shield*sqrt( diag( h_power ) )/sqrt(2)+randn(1,length(h)).*h_shield*sqrt( diag( h_power ) )/sqrt(2)*j
h = h/sqrt(h*h')% 使得每径方差的和为1，保证发射功率=接收功率，能量守恒，即h*h.'＝１;%　
h = h.' % 转成列向量。
% 产生多径瑞利衰落信道完毕。
h_Estimation = h;
H_Estimation=fft( h_Estimation,N);%产生系统函数H(w),此处N不能省略。因为h只有6行，而不是N=64行。
H=fft(h,N);%产生系统函数H(w),此处N不能省略。因为h只有6行，而不是N=64行。
%-------产生调制信号，加上CP，并进行傅立叶变换。另外，产生多径信道h(t)和系统函数H(w)-----------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fre_Offset= 0.0;   % 频偏，在时域能加小数频偏
Time_Offset=0.0;% 时偏，在时域能加小数频偏
Offset_factor_matrix = exp(i*2*pi*Fre_Offset*[0:N+CP-1]'/N)*exp(i*2*pi*Fre_Offset*Time_Offset/N);%产生在CP情况下的频偏和时偏系数矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for runningtime=1:1 %设定仿真次数。
MIN_SNR=10;   %设定最小的仿真信噪比。系统正常工作时所需的SNR与调制方式和纠检错码有关系，通常QPSK系统需要10dB, 64QAM需要20dB以上。
MAX_SNR=20;   %设定最大的仿真信噪比
BEN=zeros(1,MAX_SNR-MIN_SNR+1);%定义1×MAX_SNR的矩阵[0,0,0,....]，初始值为0，每个元素用于存放每一个SNR时的误码数
BER=zeros(1,MAX_SNR-MIN_SNR+1);%定义1×MAX_SNR的矩阵[0,0,0,....]，初始值为0，每个元素用于存放每一个SNR时的误码率
for SNR=MIN_SNR:1:MAX_SNR;
    runningtime,SNR
    for oneframe=1:FrameNumber
        inchannel=filter(h,1,xn_add_cp(:,oneframe));%让每一帧数据，依次通过多径信道h(n)。
        % f=conv(f1,f2); 求卷积，f1的长度是L1＝n1－n2,f2的长度是L2=m1－m2,f的长度是L1+L2-1
        inchannel=awgn(inchannel,SNR,'measured');%加上AWGN噪声。
        inchannel=inchannel.*Offset_factor_matrix;   % 加上频偏和时偏
        %--每一帧信号通过信道，即与信道h(n)卷积，并根据给定的SNR加上AWGN噪声，同时加上频偏.
        rec=inchannel(CP+1:N+CP);%去除CP,只留从CP+1行到最后一行
        Y=fft(rec,N);
        Y_H=Y./H;%频域均衡,去除信道系统函数H(w)的影响。”./”运算是将两个矩阵的对应元素相除
        %Y_H=Y./H_Estimation;%利用估计出的信道进行频域均衡,去除信道系统函数H(w)的影响。”./”运算是将两个矩阵的对应元素相除
        % C=Y_H-Xk(:,oneframe)  %%调试程序时使用的校验值，如果信道中没有噪声，同时没有频偏的影响。C应该等于 zeros(N,1)
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Y_H_REAL(:,oneframe)=real(Y_H); Y_H_IMAG(:,oneframe)=imag(Y_H);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        final_y=qaskdeco(real(Y_H),imag(Y_H),M);
        final_data(:,oneframe)=final_y;
        %[errorNumber,errorRatio]=biterr(Original_X(:,oneframe),final_data(:,oneframe));%Compute number of bit errors and bit error rate
        [errorNumber,errorRatio] =symerr(Original_X(:,oneframe),final_data(:,oneframe));%Compute number of symbol errors and symbol error rate
        BEN(:,SNR-MIN_SNR+1)=BEN(:,SNR-MIN_SNR+1)+errorNumber;%计算出特定SNR情况下所有帧总的误码数。
        BER(:,SNR-MIN_SNR+1)=BER(:,SNR-MIN_SNR+1)+errorRatio; %计算出特定SNR情况下所有帧总的误码率。
    end
    %[errorNumber1,errorRatio1]=biterr(Original_X,finaldata)
    %errorNumber1=errorNumber1/FrameNumber,errorRatio1=errorRatio1
end
BEN=BEN/FrameNumber;%计算出对应于每一个SNR时的每一帧的误码数。
BER=BER/FrameNumber;%计算出对应于每一个SNR时的每一帧的误码率。
% 以上程序计算出来结果。 以下绘出图形并保存结果－把结果写到动态文件里。每次动态文件名由计算结果而决定。
SNR=MIN_SNR:1:MAX_SNR;%把SNR做成1×MAX_SNR-MIN_SNR+1的矩阵。
figure(4) 
semilogy(SNR,BEN,'-*');%绘图，Plot the real-world values of fi objects on a linearly scaled x-axis and a logarithmically scaled y-axis
legend(['\epsilon=',num2str(M) ' FrameNumber=: ',num2str(FrameNumber)]);
xlabel('SNR')
ylabel('BER')
title('feedback MLSE detection over fading')
title('{\itAe}^{-\alpha\itt}sin\beta{\itt} \alpha<<\beta');%%常用的希腊符号在malab的 Text Properties上
% \alpha%\upsilon\sim~\beta\phi\leq\gamma\chi\infty\delta\psi\clubsuit\epsilon\omega\diamondsuit\zeta\Gamma\heartsuit\eta\Delta\spadesuit\theta
text(0.1,0.9, 'M=     FrameNumber=          ' , 'sc')
text(0.15,0.9, num2str(M) , 'sc')
text(0.4,0.9, int2str(FrameNumber) , 'sc')% num2str(M)
hold on   %在同一个窗口上绘下一个图
semilogy(SNR,BER,'-or');  %绘图
legend(['BEN MQAM= ',num2str(M) ' FrameNumber=: ',num2str(FrameNumber) ' Fre_Offset=  ',num2str(Fre_Offset) ' h_length=  ',num2str(length(h))],['BER Time_Offset = ',num2str(Time_Offset) ' FrameNumber=: ',num2str(FrameNumber)],['BER Time_Offset = ',num2str(Time_Offset) ' FrameNumber=: ',num2str(FrameNumber)],3);
%legend('BEN','BER',2)%图注，2 denote 2 curse in the figure.
grid on   %显示网格
if FrameNumber>=2000
    figurename = ['M_' int2str(M)  '_' 'FrameNumber_' int2str(FrameNumber) '-' int2str(floor(BEN(1))) date '.fig'];%把整型数字HSNR转化成字符串
    hgsave(['C:\MATLAB6p5\work\李军的matlab\OFDM systems framework programme\OFDM systems framework programme运行结果\' figurename]);%保存文件到制定的目录下
    matname = ['M_' int2str(M) '_' 'FrameNumber_' int2str(FrameNumber) '_' int2str(floor(BEN(1)))  date  '.mat'];%floor是地板函数
    save(['C:\MATLAB6p5\work\李军的matlab\OFDM systems framework programme\OFDM systems framework programme运行结果\' matname],'BEN','BER');%保存文件到制定的目录下
%     matname = ['M_' int2str(M) '_' 'FrameNumber_' int2str(FrameNumber) '_' int2str(floor(BEN(1)))  date  '.txt'];%floor是地板函数
%     save(['D:\' matname],'BEN','BER');%保存文件到制定的目录下 
%     save('MSE.txt','MSES');
    %close all; %把图形关闭
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5)%%% 画出星座图
modmap ('qam',M)%%% 显示MQAM的星座图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5)%%% 画出解调后的星座点
hold on
title('QPSK constellation  of demodulated signal')
plot(sqrt(M),sqrt(M),'*')
hold on
plot(-sqrt(M),-sqrt(M),'*')
hold on
% plot(Y_H_REAL,Y_H_IMAG,'.')
plot(Y_H_REAL,Y_H_IMAG,'.b')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% key=input('Want to adjust the displaying area(y/n)? [n]: ', 's' );
% if isempty(key)
%     key= 'n';
% end
% while key=='y'
%     a=111111
% end
% while key=='n'
%     a=222222
%     break
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    figure     %打开一个图形窗口
%    xlabel('SNR')
%    ylabel('BEN 或者 BER')
%    title('不同曲线分别表示特定SNR情况下所有帧总的 误码数BEN 和 误码率BNR')
%    x=0:0.05:pi;
%    y=sin(x);
%    plot(x,y,'r-- *')
%    hold on
%    plot(x,y*0.2,'b:  *')
%    text(0.8,0.55,'-* - * - *-BEN','sc')
%    gtext        ('-o-o- BER')