How can I make a mmWave frequency selective channel in matlab? $$H_d=\sum_{\mathcal{\ell}}^{N_p}\alpha_\ell p(dT_s-\tau_\ell)\mathbb{a}_R(\phi_\ell)\mathbb{a}^*_T(\theta_\ell) $$

I'm trying for simulation of a frequency selective channel in matlab. my problem is the raised cosine filter. without the raised cosine filter the matlab code is as follows. but in presence of pulse shaping filter I confused. can any body help me?

EDITED:

I've done all of Mr. Carlos Danger advises and this is my final code. but I got some errors. can any body help me?
 this is my RRC code:

function [res]=Root_Raised_Cosine(Sample_Rate,Roll_Off_Factor,time)

Sample_Interval = 1/Sample_Rate;
a = Roll_Off_Factor;
index = [-30*Sample_Interval:Sample_Interval:30*Sample_Interval]; 
%This can be increased or decreased according to the requirement

p = zeros(1,length(index));
for i = 1:1:length(index)
    if index(i) == 0
        p(i) = Sample_Rate*((1-a)+4*a/pi);
    else if index(i) == Sample_Interval/(4*a) || index(i)==-Sample_Interval/(4*a)
           p(i) = Sample_Rate*a/sqrt(2)*((1+2/pi)*sin(pi/(4*a))+(1-2/pi)*cos(pi/(4*a)));
          else
            p(i) = Sample_Rate*((sin(pi*index(i)*(1-a))+4*a*index(i)...
                .*cos(pi*index(i)*(1+a)))./(pi*index(i).*(1-(4*a*index(i)).^2)));
         end
    end
end
response = p./sqrt(sum(p.^2)); %Normalization to unit energy
t = find(abs(index-time)== Sample_Interval) ;
res = response(t);
end


and this is my channel.

clc; clear all;

%--------------------------System Parameters------------------------------%

Num_BS_Antennas = 16; % BS antennas
Num_BS_RFchains = 4;  % BS RF chains
BS_Antennas_Index = 0:1:Num_BS_Antennas-1; % Indices of the BS Antennas

Num_MS_Antennas = 16; % MS antennas
Num_MS_RFchains = 4;  %MS RF chains
MS_Antennas_Index = 0:1:Num_MS_Antennas-1; % Indices of the MS Antennas

Num_Qbits = 7; % No. of Quantization bits
Step = 2^Num_Qbits;
Step_Index = 0:1:Step-1;
PS_Quantized_Angle = Step_Index*(2*pi)/Step; % Quantized angles in PSs

%---------------------- Hybrid Precoders and Combiner---------------------%

for i = 1 : 1 : Num_BS_Antennas
for j = 1 : 1 : Num_BS_RFchains
    uniform_index = ceil(length(PS_Quantized_Angle)*rand);
    F(i,j) = 
sqrt(1/Num_BS_Antennas)*exp(1j*PS_Quantized_Angle(uniform_index)); % Hybrid 
Precoder 
end
end
clear i j
for i = 1 : 1 : Num_MS_Antennas
for j = 1 : 1 : Num_MS_RFchains
    uniform_index = ceil(length(PS_Quantized_Angle)*rand);
    W(i,j) = 
sqrt(1/Num_MS_Antennas)*exp(1j*PS_Quantized_Angle(uniform_index)); % Hybrid 
Combiner 
end
end

%--------------------------Channel Parameters ----------------------------%

Delay_Tap_Length = 4;
Sample_Rate = 1760*10^6; % Sampling rate set to be 1760 MHz
Sample_Interval = 1/Sample_Rate;
delay = 0 : Sample_Interval : (Delay_Tap_Length-1)*Sample_Interval;
Num_paths=3; % Number of channel paths
Roll_Off_Factor = .8;
% Channel parameters (angles of arrival and departure and path gains)
AoD=rand(1,Num_paths); %Sin of AoD in the range og [0,1)
AoA=rand(1,Num_paths); %Sin of AoA in the range og [0,1)
alpha=(sqrt(1/2)*sqrt(1/Num_paths)* 
(randn(1,Num_paths)+1j*randn(1,Num_paths)));

% Channel construction
 d = 3;
Channel=zeros(Num_MS_Antennas,Num_BS_Antennas);
for l=1:1:Num_paths
Abh(:,l)=sqrt(1/Num_BS_Antennas)*exp(1j*pi*BS_Antennas_Index*AoD(l));
Amh(:,l)=sqrt(1/Num_MS_Antennas)*exp(1j*pi*MS_Antennas_Index*AoA(l));
res = Root_Raised_Cosine(Sample_Rate,Roll_Off_Factor,d*Sample_Interval- 
delay(l));
Channel=Channel+ 
  (sqrt(Num_BS_Antennas*Num_MS_Antennas)*alpha(l)*Amh(:,l)*Abh(:,l)');

end
Your code appears to be summing up amplitudes over all the paths to simulate Rayleigh (flat) fading.  The assumption of flat fading is that the delay spread of the channel is much less than the symbol period, so all the symbols add up roughly at the same time (1). 

To make the leap to frequency-selective fading, you need to account for the pulse shape in your summation, i.e. you need to have some vector of sample indices n, a function p(n) that returns the pulse amplitudes at times n*Ts, and multipath delay ndelay. Something like this (pseudo-MATLAB code):

n = 0:maxn;
Npath = 4;
ndelay = [2, 5, 3, 7]; % length(ndelay) == Npath
sumval = 0;
for i = 1:Npath
    sumval = sumval + otherstuff*p(n - ndelay(i));
end


You must specify ndelay ahead of time manually or generate it according to some probability distribution.

References: 
(1) Goldsmith, Wireless Communications, chapters 2-3, 2005.