Generating GMSK signal in MATLAB

[deboleena18]

So i have written the code to generate GMSK signal, but i m not sure if it is correct
it generates a low pass envelope over the PSK signal
Also i wanted to ask if the way i m integrating to obtain phase is correct.
I am doing summation

Please Help Me

f_simulation = 20*10; %MHz
fc = 2.4*1000; %MHz 
bitrate = 1; %Mbps
bitlength = 10;
Eb = 1; %bitenergy
ep = 0; %receiver side
Ts = 1/f_sampling; %receiver side
BT = 0.5; %for gaussian filter
h_index = 0.5; %MSK modulation index




t = 0:1/f_simulation:bitlength-1/f_simulation; %time us
r = zeros(1,bitlength*f_simulation);
n = zeros(1,bitlength*f_simulation);

b = zeros(1,bitlength);
%b = randi([0,1], [1,bitlength]);
b = [1 1 0 0 1 0 1 0 1 1];
data_signal= zeros(1,bitlength*f_simulation);

k = 1;
for k = 1:bitlength
 start_index = (k-1)*f_simulation + 1;
 end_index = start_index + f_simulation - 1;
 data_signal(start_index:end_index) = 2*b(k) - 1;
end

h = gaussfir(BT, 3,f_simulation); %gaussian pule shaping filter

g = conv(data_signal, h); %g(t) the pulse signal

figure;
plot(g);

t = 1;
int = 0;
phi = zeros(1,bitlength*f_simulation);
s = zeros(1,bitlength*f_simulation);
for t = 1:bitlength*f_simulation 
    int = int + g(t)*1/f_simulation; %basically integrating upto latest time
    phi(t) = 2*pi*fc*t/f_simulation + pi*h_index*int; %phase generation
    s(t) = cos(phi(t));
    
end
figure;
plot(s);

- - - Updated - - -

okay so instead of normalized time, i decided to use real time equivalent

it is giving me an envelope free gmsk signal. I have a sampling rate of 2000 samples per bit duration ( i need the signal to look analog when i plot it )

anyways , the issue i am facing now is that the signal is clipped at the bottom

that is wrapped phase there are missing values in the range 140 to 180 degrees
this range reduces as i increase sampling rate to 20,000 per symbol

but if it was due to less samples, why did it not affect the positive peak of the signal

any help would be greatly appreciated

f_simulation = 2e9; %Hz
fc = 2.4e9; %Hz 
bitrate = 1e6; %Mbps
bitlength = 5;
bitdur=1/bitrate;
totdur=bitlength*bitdur;%%sec
totsamp=totdur*f_simulation;%%
Eb = 1e-6; %bitenergy
BT = 0.5; %for gaussian filter
h_index = 0.5; %MSK modulation index

r = zeros(1,totsamp);
n = zeros(1,totsamp);

b = rand(1,bitlength)>0.5; %% little trick 
data_signal= zeros(1,totsamp);
%b = [1 1 0 0 1 0 1 0 1 1];%% if u need fix seq

for k = 1:bitlength
 start_index = (k-1)*f_simulation*bitdur + 1;
 end_index = start_index + f_simulation*bitdur - 1;
 data_signal(start_index:end_index) = 2*b(k) - 1;
end

h = gaussfir(BT, 3,f_simulation*bitdur); %gaussian pulse shaping filter

g = conv(data_signal, h); % final signal

figure(1)
plot(g);

sum = 0;
phi = zeros(1,totsamp);
s = zeros(1,totsamp);
for t = 1:totsamp
    sum = sum + g(t+(length(h)-1)/2 )*totdur/totsamp; %basically integrating upto latest time
    phi(t) = 2*pi*fc*t*totdur/totsamp + pi*h_index*sum; %phase generation
    s(t) = cos(phi(t));
    
end
figure(2)
plot(s);