%[c] = fading_channel_jtc(f_d, f_s, n)
%
% Calculates coefficients of Rayleigh fading waveform based on 
% Noise Shaping method.
% Based on JTC fader:
%  http://www.3gpp2.org/public_html/specs/C.R1002-0_v1.0_041221.pdf
%
% Arguments:
%  f_d - doppler frequency [Hz]
%  f_s - baseband sampling frequency [Hz]
%  n   - vector of sample indices for generated coefficients
%
% Returns:
%  c   - channel coefficients

% Copyright 2017 Grzegorz Cisek (grzegorzcisek@gmail.com)

function [c] = fading_channel_jtc(f_d, f_s, n)
  assert(f_d < f_s, 'f_d must be lower than f_s');

  n_spread = max(n) - min(n) + 1;

  FIR_6  = [-0.040357044 0.086013602 0.454477919 0.454477919 0.086013602 -0.040357044];
  FIR_8  = [-0.020403315 -0.027828149 0.137993831 0.410212183 0.410212183 0.137993831 -0.027828149 -0.020403315];
  FIR_11 = [0.003612442 -0.012663859 -0.043131662 0.041846468 0.289519221 0.440818845 0.289519221 0.041846468 -0.043131662 -0.012663859 0.003612442];
  FIR_28 = [0.000134449 -0.001035837 -0.001266781 0.003144467 0.005446165 -0.005705115 -0.015697801 0.005592858 0.035808300 0.004331282 -0.073120692 -0.044582508 0.174645729 0.412153786 0.412153786 0.174645729 -0.044582508 -0.073120692 0.004331282 0.035808300 0.005592858 -0.015697801 -0.005705115 0.005446165 0.003144467 -0.001266781 -0.001035837 0.000134449];

  IIR_Den = [1.00000000000000e+000  -1.25846028151720e+001  8.37812490946412e+001  -3.87987037298430e+002  1.39276627266371e+003  -4.10390303053792e+003  1.02785179975452e+004  -2.23937486340491e+004  4.31338094397904e+004  -7.43192825675541e+004  1.15546040416494e+005  -1.63156800062187e+005  2.10262682146075e+005  -2.48183426008384e+005  2.68980386935004e+005  -2.68097215859525e+005  2.45933660734731e+005  -2.07631089086483e+005  1.61205272092231e+005  -1.14921034341049e+005  7.50416867691390e+004  -4.47318413308728e+004  2.42311152054052e+004  -1.18575082160823e+004  5.20138376926972e+003  -2.02468555919711e+003  6.90055166145190e+002  -2.02201318021456e+002  4.96491885381974e+001  -9.83333040020794e+000  1.47702790399200e+000  -1.50054529262584e-001  7.76285888645037e-003];
  IIR_Num = [6.52480599001352e-002  -5.69082890145800e-001  2.74804511668832e+000  -9.47731351802883e+000  2.57864829961265e+001  -5.82410973113121e+001  1.12471736576870e+002  -1.89048422331328e+002  2.79362373053450e+002  -3.64186311941129e+002  4.17156042029811e+002  -4.13206041327530e+002  3.39016596630252e+002  -2.00592879602055e+002  2.37345458189663e+001  1.53639128020074e+002  -2.94241547288374e+002  3.73595960603745e+002  -3.86429884358901e+002  3.45215057141779e+002  -2.72650557597993e+002  1.92305359245628e+002  -1.21539806306980e+002  6.87739305748592e+001  -3.46961260604939e+001  1.54891344545904e+001  -6.04953831961436e+000  2.03326796798172e+000  -5.74041571016860e-001  1.31218471232963e-001 -2.28674870420246e-002 2.71184861349873e-003 -1.63712912272200e-004];

  prepad = length(IIR_Den) + length(IIR_Num) + 28 + 11 + 8 + 6;

  x  = randn(prepad + ceil(n_spread/16), 1) + j * randn(prepad + ceil(n_spread/16), 1);

  x = filter(IIR_Num, IIR_Den, x);
  x = filter(FIR_28, 1, upsample(x, 2));
  x = filter(FIR_11, 1, upsample(x, 2));
  x = filter(FIR_8,  1, upsample(x, 2));
  x = filter(FIR_6,  1, upsample(x, 2));

  c = resample(16 * x, f_s, 64 * f_d);
  c = c(end-n_spread+1:end);
  c = c(n - min(n) + 1);
end