function y=ChnlGen(x, snr, chnl_mode, multipath_mode,OR,Fs,Fm)
%==============================%
%====   Author: WirelessDream     ====%
%====   Ref:IEEE 802.16.3c-01       ====%
%==============================%
% x: input signal
% snr: signal to noise ratio
% chnl_mode: channel types
%                  0: AWGN
%                 1~6:SUI 1~6
%                 7:GSM Hilly
%                 8:GSM urban
%                 9:ITU-PA
%                 10:ITU-PB
%                 11:ITU-VA
%                 12:ITU-VB
% multipath_mode: multipath generation method
%                     0:static multipath
%                     1:dynamic multipath
% OR: observation rate in Hz
% Fs: sampling rate in Hz
% Fm: maximum doppler frequency shift in Hz

M=1024;%256;     %number of taps of the Doppler filter
Dop_res=0.1;     %Doppler resolution of SUI parameter in Hz (used in resampling-process)
res_accu=20;     %accuracy of resampling process
N = length(x);%No. of independent random realizations

switch chnl_mode
    case 0
        tau = 0;
        pwr = [0.0];
        phase = 0;
    case 1               % SUI 1
        tau = [0 0.4 0.9]; % tap delay in \mu second
        pwr = [0 -15 -20]; % tap power in dB
        K = [4 0 0]; % ricean factor K (90%, omnidirectional antenna)
        Dop = [0.4 0.3 0.5]; % Doppler maximal frequency in Hz
        Fnorm = -0.1771; % gain normalization factor in dB (omnidirectional antenna)
    case 2               % SUI 2
        tau = [0 0.4 1.1]; % tap delay in \mu second
        pwr = [0 -12 -15]; % tap power in dB
        K = [2 0 0]; % ricean factor K (90%, omnidirectional antenna)
        Dop = [0.2 0.15 0.25]; % Doppler maximal frequency in Hz
        Fnorm = -0.3930; % gain normalization factor in dB (omnidirectional antenna)
    case 3               % SUI 3
        tau = [0 0.4 0.9]; % tap delay in \mu second
        pwr = [0 -5 -10]; % tap power in dB
        K = [1 0 0]; % ricean factor K (90%, omnidirectional antenna)
        Dop = [0.4 0.3 0.5]; % Doppler maximal frequency in Hz
        Fnorm = -1.5113; % gain normalization factor in dB (omnidirectional antenna)
    case 4               % SUI 4
        tau = [0 1.5 4]; % tap delay in \mu second
        pwr = [0 -4 -8]; % tap power in dB
        K = [0 0 0]; % ricean factor K (90%, omnidirectional antenna)
        Dop = [0.2 0.15 0.25]; % Doppler maximal frequency in Hz
        Fnorm = -1.9218; % gain normalization factor in dB (omnidirectional antenna)
    case 5               % SUI 5
        tau = [0 4 10]; % tap delay in \mu second
        pwr = [0 -5 -10]; % tap power in dB
        K = [0 0 0]; % ricean factor K (90%, omnidirectional antenna)
        Dop = [2 1.5 2.5]; % Doppler maximal frequency in Hz
        Fnorm = -1.5113; % gain normalization factor in dB (omnidirectional antenna)
    case 6               % SUI 6
        tau = [0 14 20]; % tap delay in \mu second
        pwr = [0 -10 -14]; % tap power in dB
        K = [0 0 0]; % ricean factor K (90%, omnidirectional antenna)
        Dop = [0.4 0.3 0.5]; % Doppler maximal frequency in Hz
        Fnorm = -0.5683; % gain normalization factor in dB (omnidirectional antenna)
    case 7    %GSM Hilly
        tau = [0.0 0.1 0.3 0.5 0.7 1.0 1.3 15.0 15.2 15.7 17.2 20.0];
        pwr = [10.0 8.0 6.0 4.0 0 0 4.0 8.0 9.0 10.0 12.0 14.0];
        phase = zeros(1,length(tau));
    case 8    %GSM Urban
        tau = [0.0 0.1 0.3 0.5 0.8 1.1 1.3 1.7 2.3 3.1 3.2 5.0];
        pwr = [ 4.0 3.0 0.0 2.6 3.0 5.0 7.0 5.0 6.5 8.6 11.0 10.0];
        phase = zeros(1,length(tau));        
    case 9       %ITU-PA
        tau = [0 0.11 0.19 0.41];
        pwr = [0 -9.7 -19.2 -22.8];
        K = zeros(1,length(tau));
        Dop = Fm*ones(1,length(tau));
        Fnorm = 0;                
    case 10      %ITU-PB
        tau = [0 0.2 0.8 1.2 2.3 3.7];
        pwr = [0 -0.9 -4.9 -8.0 -7.8 -23.9];
        K = zeros(1,length(tau));
        Dop = Fm*ones(1,length(tau));
        Fnorm = 0;
    case 11      %ITU-VA
        tau = [0 0.31 0.71 1.09 1.73 2.51];
        pwr = [0 -1 -9 -10 -15 -20];
        K = zeros(1,length(tau));
        Dop = Fm*ones(1,length(tau));
        Fnorm = 0;
    case 12      %ITU-VB
        tau = [0 0.3 8.9 12.9 17.1 20];
        pwr = [-2.5 0 -12.8 -10 -25.2 -16];
        K = zeros(1,length(tau));
        Dop = Fm*ones(1,length(tau));
        Fnorm = 0;
end   

y_noisefree = x;

tau_nor = tau*1e-6*OR; % Normalized delay time
tau_nor_min = floor(min(tau_nor));
tau2 = tau_nor-tau_nor_min;
L = length(tau);       % No. of taps

Lx = length(x);
x1 = zeros(L,Lx);
for i=1:L
    if mod(tau2(i),1)==0
        x1(i,:)=cshift(x,tau2(i)); %%circular right shift
    else
        x1(i,:)=ifft(fft(x).*exp(-j*2*pi*tau2(i)/N*[0:N-1]));%!!
    end
end
    
if (multipath_mode == 0)                       %static multipath
        P = 10.^(pwr/10);
        if (0)
            paths = sqrt(P);
        else
            s2 = P./(K+1);
            m2 = P.*K./(K+1);
            m = sqrt(m2);
            
            paths_r = sqrt(1/2)*(randn(1,L)+j*randn(1,L)).*(sqrt(s2));
            %paths_r = sqrt(1/2)*randn_complex.*(sqrt(pwr_r)); %shaoxun
            paths_c = m.*ones(1,L);
            paths = paths_r + paths_c;
        end
        paths = paths*10^(Fnorm/20);        
        y_noisefree=sum(diag(paths)*x1);   
        Chnl.PathsPwr=paths;
elseif (multipath_mode == 1)      %  dynamic multipath 
            if (1)%OnlineGen==1
                    P = 10.^(pwr/10); % calculate linear power
                    s2 = P./(K+1); % calculate variance
                    m2 = P.*(K./(K+1)); % calculate constant power
                    m = sqrt(m2);

                    %Additional Info: RMS delay spread    
                    Chnl.rmsdel = sqrt( sum(P.*(tau.^2))/sum(P) - (sum(P.*tau)/sum(P))^2 );
                    if(0)
                        fprintf('rms delay spread %6.3f μs\n', Chnl.rmsdel);
                    end

                    %create the Ricean channel coefficients with the specified powers.
                    L = length(P); % number of taps
                    paths_r = sqrt(1/2)*(randn(L,N) + j*randn(L,N)).*((sqrt(s2))' * ones(1,N));%LxN
                    paths_c = m' * ones(1,N);    
                    save paths_r paths_r paths_c;
                    paths_r = sqrt(1/2)*(randn(L,M) + j*randn(L,M)).*((sqrt(s2))' * ones(1,M));%LxM
                    paths_c = m' * ones(1,M);

                    % white spectrum is shaped according to the Doppler PSD function.
                    %Calculate time-domain filter coefficients first. The filter is then normalized in time-domain.
                    SR = max(Dop)*4; % implicit sample rate!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    for p = 1:L
                        D = Dop(p) / SR; % Doppler freq. normalized to sampling rate (4*highest Doppler)
                        f0 = [0:M*D]/(M*D)-eps; % frequency vector
                        if (0)
                            PSD = 0.785*f0.^4 - 1.72*f0.^2 + 1.0; % PSD approximation    %1x513
                        else
                            PSD = 1./sqrt(1-f0.^2);%JAKE's model
                        end
                        filt = [ PSD(1:end-1) zeros(1,M-2*floor(M*D)) PSD(end:-1:2) ]; %S(f) 1xM
                        filt = sqrt(filt); % from S(f) to |H(f)|
                        if (0)      
                            load paths_r;%%LxN
                            filt = ifftshift(ifft(filt)); % get time domain impulse response
                            filt = real(filt); %get real-valued filter
                            filt = filt / sqrt(sum(filt.^2)); % normalize filter
                            path = fftfilt(filt, [ paths_r(p,:) zeros(1,M) ]);%filters [ paths_r(p,:) zeros(1,M) ] with the FIR filter filt using the overlap/add method.
                            paths_r(p,:) = path(1+M/2:end-M/2);%!! paths_r L*N        
                        else
                            filt = filt / sqrt(sum(filt.^2));%get frequency domain transfer function
                            paths_r(p,:)=ifft(fft(paths_r(p,:)).*(filt))*sqrt(M); %!! paths_r L*M  
                        end
                    end;
                        paths = paths_r + paths_c;%paths L*M 

                    %apply the normalization factor
                    paths = paths * 10^(Fnorm/20); % multiply all coefficients with F
                    %plot(10*log10(abs(paths(1,:))));


                    %average total tap power
                    Chnl.AverTapPwr = 10*log10(mean(abs(paths).^2, 2));
                    Chnl.PSD=psd(paths(1,:), 512, max(Dop));
                    if(0)
                    fprintf('tap mean power level: %0.2f dB\n', Chnl.AverTapPwr);
                    %spectral power distribution
                    figure, Chnl.PSD;
                    end
                    % interpolate the current rate to the specified observation rate. In order to use the Matlab polyphase
                    % implementation resample, we need the resampling factor F specified as a fraction F = P/Q.
                    m = lcm(SR/Dop_res, OR/Dop_res);
                    P = m/SR*Dop_res; % find nominator (OR/Dop_res)
                    Q = m/OR*Dop_res; % find denominator  (SR/Dop_res)
                    %!!!!!!
                    N1=ceil(N*SR/OR);
                    L2=ceil(N1/2)+res_accu;%1
                    PathIdx=M/2-L2:M/2+L2;%43
                    Pathlng=length(PathIdx);
                    %!!!!!!    
                     paths_OR = zeros(L,N);
                    for p=1:L
                        if (1)
                            path_re= resample(paths(p,PathIdx), P, Q, res_accu);%1*(Pathlng*P/Q)
                        else
                            path_re= resample(paths(p,:), P, Q, res_accu);%1*(N*P/Q)  %resamples the sequencepaths(p,:) at P/Q times the original sample rate using a polyphase implementation
                        end
                        L1 = floor(length(path_re)/2-N/2);
                        paths_OR(p,:) = path_re(L1+1:L1+N);%L*N
                    end
                    save paths_OR_10dB_ITUVA_100Hz2 paths_OR; 
    else 
                load paths_OR_10dB_ITUVA_100Hz2
    end
    y_noisefree=sum(paths_OR.*x1);
    Chnl.PathsPwr=paths_OR;                
end %end of "elseif (multipath_mode == 1)"
    
ly = length(y_noisefree);
%%%%y_noisefree=y_noisefree/(norm(y_noisefree)/sqrt(ly));
Ey = norm(y_noisefree)^2 / ly;
Noise_sigma = sqrt(Ey/(10^(snr/10))/2);
%sgma = sqrt(1/(10^(snr/10))/2);
noise = Noise_sigma * (randn(1,ly) + sqrt(-1)*randn(1,ly));
y = y_noisefree + noise;