Realisation GMSK modulation in Matlab

Question

I implement a GMSK modulation scheme presented here (p 64 in pdf and copied below)

As gaussian density function I use the standard Maltab function gaussfir and convolve the rect-function with it.

My rect function:

a = randi([0,1],N,1);                       % Random data,  N x 1
ak = 2*a-1;                                 % NRZ data: +/- 1, N 
ak_1 =  kron(ones(length(ak), 1), 0);  
ak_1(1) = ak(end);
for i = 1: length(ak)-1
    ak_1(i+1) = ak(i);
end
for i = 1 : length(ak)
    akp = (-1)^i .* ak.ak_1;    % (-1)^k  d_k * d_{k-1} 
end
ak_rect = kron(akp,ones(M,1));    % M*N x 1

I have a doubt in the realisation the following equation:

My attempt:

nrz_ov_f = conv(H,ak_rect,'full');
nrz_ov_f = nrz_ov_f/max(abs(nrz_ov_f));
%integrate to get phase information
phi = filter(1,[1,-1],nrz_ov_f*Ts);  % Ts = Tb/M; M = 4- oversampling
phi = phi *0.5*pi/2;  % Tb = 2

In the equation there is a product with a_k. Does it mean I have to implement phi as follow:

phi  = phi *ak
% phi = ak * pi/2 * cumsum(phi)

?

Answer 1

The phase function is simply the Gaussian shape accumulated, and using the filter function as the OP has done is a reasonable approach to making an accumulator (digital integrator). The scaling should be such that the waveform transitions from $0$ to $\pi/2$ radians over the duration of one symbol period.

Below shows the expected waveforms of the frequency vs time (normalized Gaussian function) and phase vs time (in radians over one symbol duration) for one GMSK symbol.

The above was oversampled with 100 samples per symbol to show the underlying function. If we were to sample at 8 samples per symbol the results would appear as in the plot below:

The MATLAB code for this was:

Ts= 1; 
B = 1; 
sigma = sqrt(log(2))/(2*pi*B*Ts);
fs = 8;  % samples per symbol
t = [-1:1/fs:1-1/fs];
rect = ones(1,fs);
h = 1/(sqrt(2*pi)*sigma*Ts) .* exp(-t.^2./(2*sigma^2*Ts^2)); 
g = conv(rect/fs, h);
ph = (pi/2*Ts) cumsum(g)/fs;

As to further details into the significance that the phase transition over $\pi/2$ radians over one symbol period (either in the positive or negative direction), please refer to this post.