%Rayleigh channel function
%ch is the CSI matrix of size 90 by chnums
%chnums is the number of channel realizations.
%fd is the doppler frequency.
%
%by Yuanye Wang
%CTIF, Aalborg University
%yywang@kom.aau.dk
%01-08-2006
function ch=raych2(chnums,fd,t_rms,sampling_rate,N,cp_l,BW)
%BW=5e6;%Ts=1.008e-4;
switch t_rms
    case 2
        % this is SUI-4 model, interpolated
        t = [0,1.5,4,17.8]*1e-6; % maximum: 17.8us 
        p = [0,-4,-8,-30]; % power in dB
        % ?RMS = 1.257 ?s (before interpolation); 2.6 ?s (after interpolation)
    case 1
        % this is done by stretching SUI-3 model and interpolating it
        t = [0,1.1,2.2,6.0]*1e-6; % 
        p = [0,-5,-10,-30]; % power in dB
%         p = 10.^(p/10);
        % ?RMS = 0.997 ?s (after interpolation)       
    case 0.5
        % this is SUI-3 model, interpolated
        t = [0,0.4,0.9,3.0]*1e-6; 
        p = [0,-5,-10,-30]; % power in dB
        % ?RMS = 0.2 ?s (before interpolation); 0.4 ?s (after
        % interpolation)
end
tau = 0:1/BW/sampling_rate:t(end);
PdB = interp1(t,p,tau);     % creates the power delay profile of the interpolated channel
T_Symbol = (N+cp_l)/(BW*sampling_rate);
% F_Sampling = BW*sampling_rate;
% T_Sampling = 1/F_Sampling;
% T_Symbol = (N+cp_l)*T_Sampling;
h = rayleighchan(T_Symbol, fd, tau, PdB);  % Create channel object
%h = rayleighchan(T_sample_time, fd, tau, PdB);  % Create channel object
%reset(h);
%reset(h,round(rand(2,1)*1000));
% rand('state',sum(100*clock));
 reset(h,round(rand*1e5));
h.NormalizePathGains = true;
h.ResetBeforeFiltering = false;
h.StoreHistory = true;
x = ones(chnums, 1);
y = filter(h, x);
ch=h.PathGains';
% g=sum(abs(ch).^2,1);
% figure;plot(g);