%相關衰落信道和噪聲--Sfunction
function [sys,x0,str,ts]=cornaka_noise_threeTR(t,x,u,flag,L,N,BitRate,m,Rou_gamma_Tx,Rou_gamma_Rx,omegap,sgma)
j=sqrt(-1);
%Rou_gamma_Tx: correlation coefficient matrix of transmitter
%Rou_gamma_Rx: correlation coefficient matrix of receiver
%m: fading parameter
%omegap:  mean fading power omegap=E(x^2)
%N: number of transmit branches
%L: number of receive branches
%sgma： noise 標準差 
%h=[ h(1,1) h(1,2) ... h(1,N)
%    h(2,1) h(2,2) ... h(2,N)
%    ... ... ... ...
%    h(L,1) ... ... h(L,N)]
%g--Gaussian gam--Gamma x--Nakagami(|h|)
tic;
Rou_gamma=kron(Rou_gamma_Tx,Rou_gamma_Rx);%MIMO correlation coefficient matrix ((N*L)*(N*L))
%correlation coefficient matrix of x_vector=[x(1,1) x(2,1) ...x(L,1) x(1,2)...x(L,2)......x(1,N) ... x(L,N)]'
C_gamma=Rou_gamma*(omegap^2/m);
%covariance matrix of Gamma
C_g=sqrt(C_gamma)/(2*sqrt(m));
%covariance matrix of Gaussian

L1=chol(C_g)';%cholesky factorization.
%lower triangular matrix L

if (2*m-floor(2*m))==0 %%2*m is integer
    g=L1*randn(L*N,2*m*BitRate);%2*m*BitRate組獨立的高斯隨機向量(L*N*(2*m*BitRate))（每一組向量的協方差矩陣為C_g）
  for i=1:L*N
       gam(i,:)=sum(reshape(g(i,:).^2,2*m,BitRate));%Gamma RVs 
  end

elseif (floor(2*m)~=1)
  p=floor(2*m);
  alpha=(2*p*m+sqrt(2*p*m*(p+1-2*m)))/(p*(p+1));
  beta=2*m-p*alpha;
  g=L1*randn(L*N,p*BitRate);%p*BitRate組獨立的高斯隨機向量(L*N*(p*BitRate))（每一組向量的協方差矩陣為C_g）
  for i=1:L
       gam1(i,:)=alpha*sum(reshape(g(i,:).^2,p,BitRate));%the first of Gamma RVs 
  end
 gam=gam1+beta*(L1*randn(L*N,BitRate)).^2;%Gamma RVs
else
    p=1;
    alpha=(2*p*m+sqrt(2*p*m*(p+1-2*m)))/(p*(p+1));
    beta=2*m-p*alpha;
    g=L1*randn(L*N,BitRate);%BitRate組獨立的高斯隨機向量((L*N)*BitRate)（每一組向量的協方差矩陣為C_g）
    gam=alpha*g.^2+beta*(L1*randn(L*N,BitRate)).^2;%Gamma RVs((L*N)*BitRate) 
end

delta=2*sum((reshape(gam',L*BitRate,N))');%(1*(L*BitRate))
epsilon=normrnd(0,sqrt(delta)*sgma)+j*normrnd(0,sqrt(delta)*sgma);%(1*(L*BitRate))(noise)

Nakanew=[delta';(real(epsilon))';(imag(epsilon))'];%((3*L*BitRate)*1)符合S函數的模式

switch flag,
    case 0
        [sys,x0,str,ts]=mdlInitializeSizes(Nakanew);%初始化模塊
    case 3
        sys=mdlOutputs(t,x,u,Nakanew);%計算輸出
    case {1,2,4,9}
        sys=[];%無用的標志
    otherwise
        error(['Unhandled flag =',num2str(flag)]);%錯誤處理
end
%模塊結束
%%%%%%%初始化狀態、采樣周期、Size結構%%%%%%%%%
function [sys,x0,str,ts]=mdlInitializeSizes(Nakanew)
%調用simsizes創建sizes結構
sizes=simsizes;
%填寫sizes結構
sizes.NumContStates=0;
sizes.NumDiscStates=0;
sizes.NumOutputs=length(Nakanew);
sizes.NumInputs=1;
sizes.DirFeedthrough=1;
sizes.NumSampleTimes=1;
sys=simsizes(sizes);
x0=[];%沒有連續狀態
str=[];%沒有狀態次序
ts=[0 0];%繼承性的采樣周期

%%%%%模塊輸出%%%%%%%%%%
function sys=mdlOutputs(t,x,u,Nakanew)
sys=Nakanew.*u;