function [ratio, errors] = FMT_AWGN_BSE_DFE_SER(SNRpBit,modulation,alphabet)
% This file is used to simulate performance of the FMT system
% with different modulation method in baseband.

echo off;
% Initialize (設定)sampling times（抽樣時間）(baud rate), over sampling rate,
% modulation,alphabet（字母）, subcarrier numbers, FFT point numbers, frame
% numbers, and we'll use fractional（微小的） sampling
Fd = 2e7; Fs = 2e7; N = Fs/Fd;
if (nargin<1)
    error('No arguments are input!');
elseif (nargin>3)
    error('Too many arguments');
elseif nargin==1
    modulation='qask';
    alphabet=16;%alpha:系數a
end

subCarrNum=64;    %number of ifft/fft point ,the subchannel numbers.
gama=16;          %length of polyphase filter
NoF=1000;   %number of frame
Nf=(2*gama-1);% Select FSE and DFE lengths
Nd=gama;


% Define alphabet (quaternary)四個一組. Signal to Noise ratio (SNR).
M = alphabet; kk = log2(M);SNR = SNRpBit + 10*log10(kk);
% Set number of symbols per iteration （重復）and number of iterations
% Expected number of symbol error count
symbPerIter = subCarrNum*NoF; iters = 3; expSymErrs = 60;
numSymbTot = symbPerIter * iters;

% Set random number seeds for uniform and Gaussian noise
rand('state', 123456789);  randn('state', 987654321);%rand:產生均勻分布的隨機數
%randn:產生正態分布的隨機數
%產生噪聲，在兩個狀態內存放

% Create Gray encoding and decoding arrays（排列）
grayencod	= bitxor([0:M-1],floor([0:M-1]/2));%bitxor:按位異或，floor:向負無窮方向取整
[dummy graydecod] = sort(grayencod);  graydecod = graydecod - 1;%sort:將元素按升序排列

% Create the Prototype filter   %prototype filter:原型濾波器
Prototype=ProtoFilter(subCarrNum,gama);
% Create the polyphase components of the prototype filter and the match
% filter（匹配濾波器）
tPolyphase=reshape(Prototype,subCarrNum,gama);%polyphase:多相，gama多相濾波器的長度，subCarrNum為子載波個數, reshape:改變數組的配置 ％％％％％％一個簡單的變換，得到多相成分
mPolyphase=SubMatchFilter(tPolyphase,subCarrNum);
function mPolyphase=SubMatchFilter(tPolyphase,Row)
%this function is used to generate the polyphase components of the match filter in fmt system
%Input argument:tPolyphase——the the polyphase components of the transmit filter
%               Column    ——the column number of the input matrix.（點陣式）
%Output argument:mPolyphase——a matrix,it contains the polyphase components of the match filter.

for jj=1:Row
    mPolyphase(jj,:)=tPolyphase(Row-jj+1,:);
end
for ii=1:subCarrNum
    subChannelMat(ii,:)=conv(tPolyphase(ii,:),mPolyphase(ii,:));%conv:求兩個多項式的乘積
end
for ii=1:subCarrNum
    C(ii,:)=subChannelMat(ii,:)/norm(subChannelMat(ii,:));%%%norm:矩陣或向量的范數
end
%%%此處C矩陣是什么作用 ？？？？
% because FIR is used here , so the delta1 is in the middle
for ii=1:subCarrNum
    [Y(ii),delta1(ii)]=max(C(ii,:));  
end

% to channel-FSE combination delay
% decision need 1 delay
delta2=(delta1+1);

% Generate random numbers from in the range [0,M-1]
msg_orig = randsrc(symbPerIter,1,[0:M-1]);

% Gray encode symbols
if (strcmp(modulation,'psk'))
    msg_gr_orig = grayencod(msg_orig+1)';
else
    msg_gr_orig=msg_orig;
end

% Map the digital signal to analog signal
msg_map_gry_orig=modmap(msg_gr_orig,Fd,Fs,modulation,M);      
msg_to_ifft=msg_map_gry_orig(:,1)+i.*msg_map_gry_orig(:,2);%將矩陣變為復數形式
msg_tx=ifft(reshape(msg_to_ifft,subCarrNum,N*NoF));

for k=1:length(SNR),  
    % Normalize channel for unit received power 
%     snr=10^(SNR(k)/10);           %這歸化法有到據，能源序用調方有，們里時掉
%     desirednorm=2*(1-1/snr);
%     h=C*sqrt(desirednorm); 
    h=C;

    % Find UMMSE-FSE+DFE coefficients 
    for ii=1:subCarrNum
        [Ecf(:,ii),Ecd(:,ii),mmse(k,ii)]=ummse_fse_dfe(h(ii,:),SNR(k),Nf,Nd,delta1(ii),delta2(ii)); %設計均衡器，Ecf，Ecd，誤差在mmse中,h為信道參數矩陣
    end
    ser(k)=qam_ser(M,mmse(k,1));  %theoretical  %該函數是跟據符號SNR計算理論誤碼性能的！
    
    % Get dd
    for(ii=1:subCarrNum)
        fh(:,ii)=conv(Ecf(:,ii),h(ii,:)).';   % the convolution of the feed forward filter and the equivalent channel
        dd(:,ii)=[zeros(delta2(ii)-1,1); fh(delta2(ii):delta2(ii)+Nd,ii); ...
                zeros(size(fh,1)-Nd-delta2(ii),1)];
    end
    
    % transmit the mapped signal over the equivalent channel(tPolyphase+mPolyphase+FFF+FBF)
    for (ii=1:subCarrNum)
        n(ii,:)=awgn(msg_tx(ii,:),SNR(k),'measured',[],'dB')-msg_tx(ii,:);%得到要加入的高斯白噪聲
        msg_rx1(ii,:)=filter(h(ii,:),1,msg_tx(ii,:))+n(ii,:);%  the received signals without equalization,在濾波之前加噪聲還是在濾波之后加？
        msg_rx(ii,:)=msg_rx1(ii,201:N*NoF);% remove the first 200 signals
        nn(ii,:)=filter(Ecf(:,ii),1,n(ii,:));% the output noise of the feedforward filter
        hall=(fh-dd);
        msg_rx2(ii,:)=conv(hall(:,ii),msg_tx(ii,:));
        msg_rx_equlized(ii,:)=msg_rx2(ii,gama:N*NoF+gama-1)+nn(ii,:);   
    end  
    %通過一系列的變換得到均衡后的信號
    
    
    
    %fft and s/p
    msg_rx_equalized_fft=fft(msg_rx_equlized);
    msg_rx_noequ_fft=fft(msg_rx1);%沒有均衡的結果
    msg_to_demodmap=reshape(msg_rx_equalized_fft,subCarrNum*NoF*N,1);
%    msg_to_demodmap=reshape(msg_rx_noequ_fft,subCarrNum*NoF*N,1);
    msg_to_demodmap2=[real(msg_to_demodmap) imag(msg_to_demodmap)];
    
    % Demodmap
    msg_gr_demod=demodmap(msg_to_demodmap2,Fd,Fs,modulation,M);
    
    % Gray decode message
    if (strcmp(modulation,'psk'))
        msg_demod = graydecod(msg_gr_demod+1)';
    else
        msg_demod=msg_gr_demod;
    end  
            
    % calculate bit error count, BER, symbol error count and SER, for this iteration.
    [errBit(k) ratBit(k)] = biterr(msg_orig, msg_demod, kk);
    [errSym(k) ratSym(k)] = symerr(msg_orig, msg_demod);
    err(k)=ratSym(k);
end  
aa=ser
bb=err
errors=errBit;
ratio=ratBit;
save FMT_AWGN_BSE_DFE_SER;
% plot the scatterplot
%figure(1);
%if k==length(SNR)
 %   subplot(121)
  %  plot(msg_rx_noequ_fft,'.');
  %  title('received');
  %  axis('square');
  %  xlabel('Re');
  %  ylabel('Im');
  %  subplot(122)
  %  plot(msg_to_demodmap,'.');
  %  title('equalized');
  %  axis('square');
  %  xlabel('Re');
  %  ylabel('Im');

  %end
% plot the ser performance
figure(1);
subplot(2,1,1);
semilogy(SNRpBit,err,'k-*');
xlabel('Eb/N0 (dB)');
ylabel('SER');
title(['simulative SER']);
%figure(3);
grid;
hold on;
subplot(2,1,2);
semilogy(SNRpBit,ser,'b-o');
xlabel('Eb/N0 (dB)');
ylabel('SER');
title(['theoretical SER']);
grid;
figure(3)
reshape(C,64*31,1);
plot(C);
%plot(ifft(Prototype));
%plot(Prototype);