function RayleighChannel = my_rayleigh(f_max,SymbolNum,snr, signalpower)
% f_max = 50;         % 最大多普勒頻偏
% SymbolNum =1;
% snr =20;

fband = 20*10^6;    % 基帶頻率20MHz
T_s = 1/fband;      % 基帶采樣時鐘周期間隔
% T_symbol = 64*T_s;  % 符號周期
% T_sim = SymbolNum*T_symbol  %仿真時間
K = 80*SymbolNum;   % 采樣點數


t = 0 : T_s :(K-1)*T_s;  


% f1, c1, th1: discrete Doppler frequencies, Doppler coefficients,and Doppler phases of mu_1(t)
% f2, c2, th2: discrete Doppler frequencies, Doppler coefficients,and Doppler phases of mu_2(t)
N_i = [64 64]; % N_i(1)、N_i(2)分別為實部虛部高斯過程的諧波個數  
%sigma_0 = sqrt(signalpower*10^(-snr/10));
%--------------------------參數產生-----------------------------------------
                   % Method of exact Doppler spread (MEDS)
% n=(1:N_i(1))';
% f1 =f_max*sin(pi/(2*N_i(1))*(n-1/2));
% c1 =sigma_0*sqrt(2/(N_i(1)))*ones(size(f1)); 
% Z=rand(size(n));
% [dummy,I]=sort(Z);
% % th1 = 2*pi*n(I)/(N_i(1)+1);
% th1 = rand(N_i(1),1)*2*pi;
                   % Method of equal distances (MED)
% n = 1:N_i(1);
% f1 = f_max/(2*N_i(1))*(2*n-1);
% c1 = 2*sigma_0/sqrt(pi)*(asin(n/N_i(1))-asin((n-1)/N_i(1))).^0.5;
% th1 = rand(N_i(1),1)*2*pi;
                  % Jakes method
n = 1:N_i(1)-1;
f1 = f_max*[[cos(pi*n/(2*(N_i(1)-1/2))),1]',...
                  [cos(pi*n/(2*(N_i(1)-1/2))),1]'];
%c1 = 2*sigma_0/sqrt(N_i(1)-1/2)*[[sin(pi*n/(N_i(1)-1)),1/2]',...
%                                 [cos(pi*n/(N_i(1)-1)),1/2]'];

c1 = 2/sqrt(N_i(1)-1/2)*[[sin(pi*n/(N_i(1)-1)),1/2]',...
                                 [cos(pi*n/(N_i(1)-1)),1/2]'];

th1 = zeros(size(f1));
%----------------------
                 % Method of exact Doppler spread (MEDS)
% n=(1:N_i(2))';
% f2 =f_max*sin(pi/(2*N_i(2))*(n-1/2));
% c2 =sigma_0*sqrt(2/(N_i(2)))*ones(size(f2));
% Z=rand(size(n));
% [dummy,I]=sort(Z);
% % th2 = 2*pi*n(I)/(N_i(2)+1);
% th2 = rand(N_i(2),1)*2*pi;
                 % Method of equal distances (MED)
% n = 1:N_i(2);
% f2 = f_max/(2*N_i(2))*(2*n-1);
% c2 = 2*sigma_0/sqrt(pi)*(asin(n/N_i(2))-asin((n-1)/N_i(2))).^0.5;
% th2 = rand(N_i(2),1)*2*pi;
%                  Jakes method 
n = 1:N_i(2)-1;
f2 = f_max*[[cos(pi*n/(2*(N_i(2)-1/2))),1]',...
                  [cos(pi*n/(2*(N_i(2)-1/2))),1]'];
%c2 = 2*sigma_0/sqrt(N_i(2)-1/2)*[[sin(pi*n/(N_i(2)-1)),1/2]',...
%                                 [cos(pi*n/(N_i(2)-1)),1/2]'];
c2 = 2/sqrt(N_i(2)-1/2)*[[sin(pi*n/(N_i(2)-1)),1/2]',...
                                 [cos(pi*n/(N_i(2)-1)),1/2]'];


th2 = zeros(size(f2));
%---------------------------生成兩個高斯過程--------------------------------
mu_1_t=0;
for k=1:length(f1),
    mu_1_t=mu_1_t+c1(k)*cos(2*pi*f1(k)*t+th1(k));
end
%----------------------
mu_2_t=0;
for k=1:length(f2),
    mu_2_t=mu_2_t+c2(k)*cos(2*pi*f2(k)*t+th2(k));
end
%----------------------------合成瑞利分布-----------------------------------
 RayleighChannel = (mu_1_t+j*mu_2_t)/sqrt(2);
 
% 畫功率譜密度 
% y=zeros(1,K);
% for n=1:K
%     y(n)= RayleighChannel*[zeros(1,n),RayleighChannel(1:K-n)].';
% end
% psd = fft(y,K);
% stem(1:K,abs(psd));


%  figure()
%  subplot(3,1,1),plot(t,abs(mu_1_t));
%  subplot(3,1,2),plot(t,abs(mu_2_t));
%  subplot(3,1,3),plot(t,abs(RayleighChannel));



%%=======================參考程序============================================
% M=16;
% N=4*M+2;
% 
% fband=20*10^6;  %%基帶采樣頻率5MHz
% Wm=2*pi*f_max;    %%最大多普勒頻移
% A(1:M)=2*pi/N.*(1:M);  %%入射波功率均勻分布在2*pi內
% Wn(1:M)=Wm*cos(A(1:M));  %%入射波多普勒頻移
% B(1:M)=pi/M*(1:M);
% B0=pi/4;%%B(M+1)
% C=2*sqrt(1/N);
% 
% Tsample=1/fband; %%基帶采樣時鐘周期間隔
% Tsymbol=80*Tsample; %%符號周期
% 
% t=0:Tsample:SymbolNum*Tsymbol;  
% K=80*SymbolNum;     %%采樣點數
% MatrixB1(1:M)=cos(B(1:M));
% MatrixB2(1:M)=sin(B(1:M));
% for index2=1:K;
%    for index1=1:M;
%         MatrixWn(index1,index2)=cos(Wn(index1)*t(index2));
%     end;
% end;
% 
% in_temp1=MatrixB1*MatrixWn;
% in_temp2(1:K)=sqrt(2)*cos(B0)*cos(Wm*t(1:K));
% qua_temp1=MatrixB2*MatrixWn;
% qua_temp2(1:K)=sqrt(2)*sin(B0)*cos(Wm*t(1:K));
% 
% Xc=2/(sqrt(N))*(in_temp2+2*in_temp1);
% Xs=2/(sqrt(N))*(qua_temp2+2*qua_temp1);
% 
% RayleighChannel=Xc+j*Xs;