clc 
close all

fs=10E9 ;
ts=1/fs;
mu=2; % QAM modulation index
m=2^mu; % QAM modulation number
K=128; % subcarrier number
M=5;  % time slot number
CP_rate=1/4;  %
b_length=mu*K*M;  %length of data block
size_over=b_length;
N=K;  %upsample rate (N  must >= K)
ifft_length=N*M;  % fft length
alpha=0.1; % roll-off factor of filter


blockLen = gfdm.M*gfdm.K;
EbN0Vec=(0:1:10);
% for si=1:length(EbN0)
%             for i=1:its
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TRANSMITTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D=zeros(K,Nch*M);
iter=2;
for i=1:iter   
    disp(i);
for kk=1:length(EbN0Vec)
     snrdb=EbN0Vec(kk)+10*log10(mu);
     
 
    numErrsmf=0;
    numErrszf=0;
    numErrsmmse=0;
    numErrsmf1=0;
    numErrszf1=0;
    numErrsmmse1=0;
    
   
%%%%%%%%%%%%%%%%%%GFDM
b=randi([0,1],1,b_length);% generate bit flow

dataIn=reshape(b,length(b)/mu,mu);
dataSymbol=bi2de(dataIn);
dataSymbol=dataSymbol';
d=qammod(dataSymbol,4); % 16QAM moduation to get d(S/P)
%   scatterplot(d);
D=reshape(d,[K,M]); % get data block D


%%% Upsampling
s=zeros(ifft_length,M); % s will be a sampling matrix
for ii=1:ifft_length
    for jj=1:M
        if ii==(jj-1)*N+1
            s(ii,jj)=1;
        end
    end
end


%%% Pulse shaping 

filter_h=rcosdesign(alpha,M,N);
 filter_h=filter_h(1:length(filter_h)-1); % filter_h is the coefficient of the filter
filter_column=filter_h';
G_Tx=circulant(filter_column,1);
gtx=G_Tx*s;

% %%%%%%%%%%%%%%%fft%%%%%%%%%%%%%%%%%%%%

Nn = 1:1:K*M; % row vector for n
Kk = 1:1:K; % ro0w vecor for k
WN1 = exp(-1i*2*pi/(N)); % Twiddle factor (w)
kn = Kk'*Nn; % creates a N by N matrix of nk values
WNnk1 = (WN1 .^ kn);

 W=zeros(K,M*K);
 
for k=1:1:K
    for n=1:1:K*M
        W(k,n)=exp(-1i*2*pi*k*n/N);
    end
end


WNnk1=ctranspose(WNnk1);
W=ctranspose(W);
 
% wh²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²²
% 

A=zeros(M*K);

 for n= 1:1:M*K
%     for j=1:M*K 
    A(n,:)= kron(WNnk1(n,:),gtx(n,:));  
 end




% dd=s2*ctranspose(W);
data=reshape(D,M*K,1);
% d{k}=D{k}(:);
x_trans=A*data;


% 

% CP_adding

x_whole=[x_trans(ifft_length*(1-CP_rate)+1:ifft_length);x_trans];



% create_field('unique',x_whole,[],struct('power','average'));


 Indata=awgn(x_whole,snrdb,'measured');

%%%%%%%%%%%%%%%CP remove

y_cp_awgn=Indata(CP_rate*ifft_length+1:CP_rate*ifft_length+ifft_length);
%  y_cp=y_cp.';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  AWGN
%%%%%% Equalization



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GFDM
%%%Matched Filter

A_MF=ctranspose(A);

d_MF=A_MF*y_cp_awgn;



%%Zero Forcing
a_zf=ctranspose(A)*A;
% a_zf2=inv(a_zf);
A_ZF=a_zf\ctranspose(A);
d_ZF=A_ZF*y_cp_awgn;


%%MMSE

var_data=var(data)^2;
var_noise=(10.^(-snrdb/10))^2;
p1=(var_noise/var_data)*eye(M*K);
p2=ctranspose(A)*A;
p3=(p1+p2)\ctranspose(A);
d_MMSE=p3*y_cp_awgn;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Demodulation  %%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% ZF
%  scatterplot(d_MF);
dem_ZF=qamdemod(d_ZF,4);
b_ZF=de2bi(dem_ZF);
b_ZF=b_ZF(:);
b_ZF=b_ZF.';

%%%%%%%%%%%%%%%%%%%% MF
% scatterplot(d_MF);
dem_MF=qamdemod(d_MF,4);
b_MF=de2bi(dem_MF);
b_MF=b_MF(:);
b_MF=b_MF.';

%%%%%%%%%%%%%%%%%MMSE
dem_MMSE=qamdemod(d_MMSE,4);
b_MMSE=de2bi(dem_MMSE);
b_MMSE=b_MMSE(:);
b_MMSE=b_MMSE.';







% bermfawgn=zeros(1,length(EbN0Vec));
nErrorsmf=biterr(b,b_MF);
numErrsmf=numErrsmf+nErrorsmf;
bermfawgn_gfdm(kk)=numErrsmf/b_length;


% berzfawgn=zeros(1,length(EbN0Vec));
nErrorszf=biterr(b,b_ZF);
numErrszf=numErrszf+nErrorszf;
berzfawgn_gfdm(kk)=numErrszf/b_length;

% bermmseawgn=zeros(1,length(EbN0Vec));
nErrorsmmse=biterr(b,b_MMSE);
numErrsmmse=numErrsmmse+nErrorsmmse;
bermmseawgn_gfdm(kk)=numErrsmmse/b_length;

end

mf_gfdm(1,:)=bermfawgn_gfdm;
mf_gfdm=[mf_gfdm;bermfawgn_gfdm];

zf_gfdm(1,:)=berzfawgn_gfdm;
zf_gfdm=[zf_gfdm;berzfawgn_gfdm];

 mmse_gfdm(1,:)=bermmseawgn_gfdm;
mmse_gfdm=[mmse_gfdm;bermmseawgn_gfdm];

end


mf_gfdm=mf_gfdm(1:iter,:);
zf_gfdm=zf_gfdm(1:iter,:);
mmse_gfdm=mmse_gfdm(1:iter,:);


mean_mf_gfdm=mean(mf_gfdm);
mean_zf_gfdm=mean(zf_gfdm);
mean_mmse_gfdm=mean(mmse_gfdm);



 berth=berawgn(EbN0Vec,'qam',4);
 
 semilogy(EbN0Vec,berth,'b*-')
 hold on

semilogy(EbN0Vec,mean_mf_gfdm,'go-')
 hold on
 semilogy(EbN0Vec,mean_zf_gfdm,'b^-')
  hold on
  semilogy(EbN0Vec,mean_mmse_gfdm,'ro-')
%     hold on
%  semilogy(EbN0Vec,bermffiber,'co-')
%     hold on
%   semilogy(EbN0Vec,berzffiber,'g^-')
  legend('QAM BER Theory','GFDM MF  ','GFDM ZF ','GFDM MMSE ')

%   legend('QAM BER Theory','GFDM MF awgn ','GFDM ZF awgn','GFDM MMSE awgn','GFDM MF fiber 80km ','GFDM ZF fiber 80 km')
title(' QAM  awgn  K=128 M=5 (google patent)');
xlabel('EbN0(db)');
ylabel('BER');

 figure(5)
  plot(mag2db(fftshift(abs(fft((x_trans), 2*length(x_trans)))))/2, 'r');
  