%%%%%%%%%%%%%%    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')