I was able to eventually figure this out. I happened to have both books by Rappaport and Anderson, Aulin, and Sundberg so I was able to see the derivations for both pulse shapes. The answer for how they are linked came from stumbling across a key section in Marvin Simon's Bandwidth-Efficient Digital Modulation with Application to Deep-Space Communications paper for JPL.
The definition in Anderson, et al. is:
While the definition in Rappaport is:
The key difference as noted by Simon (2.8.2.1) is that g(t) is the generating frequency pulse while h(t) is the filter impulse response resulting from an NRZ data stream "i.e., the frequency pulse stream that would ordinarily be inputted to the FM modulator in MSK".
The frequency pulses with proper scaling are:
Which the cumulative-sum shows both end at 1/2, which is the definition for phase functions in CPM:
So the question of which one to use depends on how the source data is formatted. If the source is a impulse train (i.e. bit stream), then use comm.GMSKModulator. If you have a NRZ stream at the Ts rate then use gaussdesign. Here is an example comparing them all:
bt = 0.5; % 3-dB bandwidth-symbol time
L = 12; % samples per symbol
N = 4; % filter spans 4 symbols
h = gaussdesign(bt,N,L);
h = h/2;
% g was obtained from cpmmodparams.m
H = comm.GMSKModulator;
H.BandwidthTimeProduct = 0.5;
H.PulseLength = 4;
H.SamplesPerSymbol = 12;
H.BitInput = false; % impulses of +1/-1
data = 2*randi([0 1],10*L,1)-1; % generate the NRZ data stream
data2 = zoh(data,L); % zero-order hold at Ts
x = filter(h,1,data2)/L; % filter with the Rappaport pulse
y = upfirdn(data,g,L,1); % upsample the impulses with the Anderson pulse
modSignal = step(H, data); % modulate using comm.GMSKModulator
z = cpfsk_demod(modSignal)/pi; % demod to get the frequency signal
subplot(3,1,1); plot(x(1:1000)); title('Rappaport w/ NRZ data');
subplot(3,1,2); plot(y(1:1000)); title('Anderson w/ impulse stream');
subplot(3,1,3); plot(z(1:1000)); title('comm.GMSKModulator');