diff --git a/Transceiver52M/laurent.m b/Transceiver52M/laurent.m
new file mode 100644
index 00000000..ef154283
--- /dev/null
+++ b/Transceiver52M/laurent.m
@@ -0,0 +1,83 @@
+%
+% Laurent decomposition of GMSK signals
+% Generates C0, C1, and C2 pulse shapes
+%
+% Pierre Laurent, "Exact and Approximate Construction of Digital Phase
+%   Modulations by Superposition of Amplitude Modulated Pulses", IEEE
+%   Transactions of Communications, Vol. 34, No. 2, Feb 1986.
+%
+% Author: Thomas Tsou <tom@tsou.cc>
+%
+
+% Modulation parameters
+oversamp = 16;
+L = 3;
+f = 270.83333e3;
+T = 1/f;
+h = 0.5;
+BT = 0.30;
+B = BT / T;
+
+% Generate sampling points for L symbol periods
+t = -(L*T/2):T/oversamp:(L*T/2);
+t = t(1:end-1) + (T/oversamp/2);
+
+% Generate Gaussian pulse
+g = qfunc(2*pi*B*(t - T/2)/(log(2)^.5)) - qfunc(2*pi*B*(t + T/2)/(log(2)^.5));
+g = g / sum(g) * pi/2;
+g = [0 g];
+
+% Integrate phase 
+q = 0;
+for i = 1:size(g,2);
+    q(i) = sum(g(1:i));
+end
+
+% Compute two sided "generalized phase pulse" function
+s = 0;
+for i = 1:size(g,2);
+    s(i) = sin(q(i)) / sin(pi*h);
+end
+for i = (size(g,2) + 1):(2 * size(g,2) - 1);
+    s(i) = sin(pi*h - q(i - (size(g,2) - 1))) / sin(pi*h);
+end
+
+% Compute C0 pulse: valid for all L values
+c0 = s(1:end-(oversamp*(L-1)));
+for i = 1:L-1;
+    c0 = c0 .* s((1 + i*oversamp):end-(oversamp*(L - 1 - i)));
+end
+
+% Compute C1 pulse: valid for L = 3 only!
+%   C1 = S0 * S4 * S2
+c1 = s(1:end-(oversamp*(4)));
+c1 = c1 .* s((1 + 4*oversamp):end-(oversamp*(4 - 1 - 3)));
+c1 = c1 .* s((1 + 2*oversamp):end-(oversamp*(4 - 1 - 1)));
+
+% Compute C2 pulse: valid for L = 3 only!
+%   C2 = S0 * S1 * S5
+c2 = s(1:end-(oversamp*(5)));
+c2 = c2 .* s((1 + 1*oversamp):end-(oversamp*(5 - 1 - 0)));
+c2 = c2 .* s((1 + 5*oversamp):end-(oversamp*(5 - 1 - 4)));
+
+% Plot C0, C1, C2 Laurent pulse series
+figure(1);
+hold off;
+plot((0:size(c0,2)-1)/oversamp - 2,c0, 'b');
+hold on;
+plot((0:size(c1,2)-1)/oversamp - 2,c1, 'r');
+plot((0:size(c2,2)-1)/oversamp - 2,c2, 'g');
+
+% Generate OpenBTS pulse
+numSamples = size(c0,2); 
+centerPoint = (numSamples - 1)/2;
+i = ((0:numSamples) - centerPoint) / oversamp;
+xP = .96*exp(-1.1380*i.^2 - 0.527*i.^4);
+xP = xP / max(xP) * max(c0);
+
+% Plot C0 pulse compared to OpenBTS pulse
+figure(2);
+hold off;
+plot((0:size(c0,2)-1)/oversamp, c0, 'b');
+hold on;
+plot((0:size(xP,2)-1)/oversamp, xP, 'r');
diff --git a/Transceiver52M/sigProcLib.cpp b/Transceiver52M/sigProcLib.cpp
index 3f34fa2a..4c7d1e62 100644
--- a/Transceiver52M/sigProcLib.cpp
+++ b/Transceiver52M/sigProcLib.cpp
@@ -420,11 +420,22 @@ void generateGSMPulse(int sps, int symbolLength)
   GSMPulse->empty->isRealOnly(true);
   *(GSMPulse->empty->begin()) = 1.0f;
 
-  len = sps * symbolLength;
-  if (len < 4)
-    len = 4;
+  /*
+   * For 4 samples-per-symbol use a precomputed single pulse Laurent
+   * approximation. This should yields below 2 degrees of phase error at
+   * the modulator output. Use the existing pulse approximation for all
+   * other oversampling factors.
+   */
+  switch (sps) {
+  case 4:
+    len = 16;
+    break;
+  default:
+    len = sps * symbolLength;
+    if (len < 4)
+      len = 4;
+  }
 
-  /* GSM pulse approximation */
   GSMPulse->buffer = convolve_h_alloc(len);
   GSMPulse->gaussian = new signalVector((complex *)
                                         GSMPulse->buffer, 0, len);
@@ -433,12 +444,32 @@ void generateGSMPulse(int sps, int symbolLength)
 
   signalVector::iterator xP = GSMPulse->gaussian->begin();
 
-  center = (float) (len - 1.0) / 2.0;
+  if (sps == 4) {
+    *xP++ = 4.46348606e-03;
+    *xP++ = 2.84385729e-02;
+    *xP++ = 1.03184855e-01;
+    *xP++ = 2.56065552e-01;
+    *xP++ = 4.76375085e-01;
+    *xP++ = 7.05961177e-01;
+    *xP++ = 8.71291644e-01;
+    *xP++ = 9.29453645e-01;
+    *xP++ = 8.71291644e-01;
+    *xP++ = 7.05961177e-01;
+    *xP++ = 4.76375085e-01;
+    *xP++ = 2.56065552e-01;
+    *xP++ = 1.03184855e-01;
+    *xP++ = 2.84385729e-02;
+    *xP++ = 4.46348606e-03;
+    *xP++ = 0.0;
+  } else {
+    center = (float) (len - 1.0) / 2.0;
 
-  for (int i = 0; i < len; i++) {
-    arg = ((float) i - center) / (float) sps;
-    *xP++ = 0.96 * exp(-1.1380 * arg * arg -
-                        0.527 * arg * arg * arg * arg);
+    /* GSM pulse approximation */
+    for (int i = 0; i < len; i++) {
+      arg = ((float) i - center) / (float) sps;
+      *xP++ = 0.96 * exp(-1.1380 * arg * arg -
+			 0.527 * arg * arg * arg * arg);
+    }
   }
 
   float avgAbsval = sqrtf(vectorNorm2(*GSMPulse->gaussian)/sps);