5 bit ramp-up/ramp-down

CommonLibs/BitVector.cpp

@@ -30,6 +30,7 @@
 #include <iostream>
 #include <stdio.h>
 #include <sstream>
+#include <math.h>
 
 using namespace std;
 
@@ -533,7 +534,8 @@ float SoftVector::getEnergy(float *plow) const
 	float avg = 0; float low = 1;
 	for (int i = 0; i < len; i++) {
 		float bit = vec[i];
-		float energy = 2*((bit < 0.5) ? (0.5-bit) : (bit-0.5));
+		float energy = 2*bit-1.0;
+		energy *= energy;
 		if (energy < low) low = energy;
 		avg += energy/len;
 	}

Transceiver52M/Transceiver.cpp

@@ -621,14 +621,28 @@ int Transceiver::detectTSC(TransceiverState *state, signalVector &burst,
   return 1;
 }
 
+void writeToFile(signalVector *burst, const GSM::Time &time, size_t chan, const std::string postfix="")
+{
+  std::ostringstream fname;
+  fname << chan << "_" << time.FN() << "_" << time.TN() << postfix << ".fc";
+  std::ofstream outfile(fname.str().c_str(), std::ofstream::binary);
+  outfile.write((char*)burst->begin(), burst->size() * 2 * sizeof(float));
+  outfile.close();
+}
+
 /*
  * Demodulate GMSK burst using equalization if requested. Otherwise
  * demodulate by direct rotation and soft slicing.
  */
 SoftVector *Transceiver::demodulate(TransceiverState *state,
                                     signalVector &burst, complex amp,
-                                    float toa, size_t tn, bool equalize)
+                                    float toa, size_t tn, bool equalize,
+                                    GSM::Time &wTime, size_t chan,
+                                    Transceiver::BurstQuality *qual)
 {
+  signalVector *aligned, *bit_aligned=NULL;
+  SoftVector *bits;
+
   if (equalize) {
     scaleVector(burst, complex(1.0, 0.0) / amp);
     return equalizeBurst(burst,
@@ -638,17 +652,28 @@ SoftVector *Transceiver::demodulate(TransceiverState *state,
                          *state->DFEFeedback[tn]);
   }
 
-  return demodulateBurst(burst, mSPSRx, amp, toa);
-}
+  aligned = alignBurst(burst, amp, toa);
 
-void writeToFile(radioVector *radio_burst, size_t chan)
-{
-  GSM::Time time = radio_burst->getTime();
-  std::ostringstream fname;
-  fname << chan << "_" << time.FN() << "_" << time.TN() << ".fc";
-  std::ofstream outfile (fname.str().c_str(), std::ofstream::binary);
-  outfile.write((char*)radio_burst->getVector()->begin(), radio_burst->getVector()->size() * 2 * sizeof(float));
-  outfile.close();
+  if (qual) {
+    /* "aligned" burst has samples exactly between bits.
+     * Delay it by 1/2 bit more to get samples aligned to bit positions. */
+    bit_aligned = delayVector(aligned, NULL, 0.5);
+
+    /* Debug: dump bursts to disk */
+    if (needWriteBurstToDisk(wTime, chan))
+      writeToFile(bit_aligned, wTime, chan, "_aligned");
+  }
+
+  bits = demodulateBurst(*aligned, mSPSRx);
+
+  if (qual) {
+    /* Estimate signal quality */
+    estimateBurstQuality(bits->segment(0, gSlotLen).sliced(), bit_aligned, wTime, chan, *qual);
+    delete bit_aligned;
+  }
+
+  delete aligned;
+  return bits;
 }
 
 /*
@@ -657,7 +682,8 @@ void writeToFile(radioVector *radio_burst, size_t chan)
  */
 SoftVector *Transceiver::pullRadioVector(GSM::Time &wTime, double &RSSI, bool &isRssiValid,
                                          double &timingOffset, double &noise,
-                                         size_t chan)
+                                         size_t chan,
+                                         Transceiver::BurstQuality *qual)
 {
   int success;
   bool equalize = false;
@@ -679,10 +705,8 @@ SoftVector *Transceiver::pullRadioVector(GSM::Time &wTime, double &RSSI, bool &i
   CorrType type = expectedCorrType(time, chan);
 
   /* Debug: dump bursts to disk */
-  /* bits 0-7  - chan 0 timeslots
-   * bits 8-15 - chan 1 timeslots */
-  if (mWriteBurstToDiskMask & ((1<<time.TN()) << (8*chan)))
-    writeToFile(radio_burst, chan);
+  if (needWriteBurstToDisk(time, chan))
+    writeToFile(radio_burst->getVector(), time, chan);
 
   /* No processing if the timeslot is off.
    * Not even power level or noise calculation. */
@@ -754,7 +778,7 @@ SoftVector *Transceiver::pullRadioVector(GSM::Time &wTime, double &RSSI, bool &i
   if (equalize && (type != TSC))
     equalize = false;
 
-  bits = demodulate(state, *burst, amp, toa, time.TN(), equalize);
+  bits = demodulate(state, *burst, amp, toa, time.TN(), equalize, time, chan, qual);
 
   delete radio_burst;
   return bits;
@@ -978,6 +1002,125 @@ void Transceiver::driveReceiveRadio()
   }
 }
 
+inline float wrapAngle2Pi(float angle)
+{
+  const float twoPi = 2.0 * M_PI;
+  return angle - twoPi * floor( angle / twoPi );
+}
+
+inline float wrapAnglePi(float angle)
+{
+  const float twoPi = 2.0 * M_PI;
+  return angle - twoPi * floor( (angle+M_PI) / twoPi);
+}
+
+inline float rad2deg(float rad)
+{
+  return rad*180/M_PI;
+}
+
+inline int vectorMaxAbs(const Vector<float> &vec)
+{
+  int max_idx = 0;
+  float max = 0.0;
+  for (size_t i=1; i<vec.size(); i++) {
+    if (fabs(vec[i]) > max) {
+      max_idx = i;
+      max = fabs(vec[i]);
+    }
+  }
+  return max_idx;
+}
+
+inline float vectorRMS(const Vector<float> &vec)
+{
+  float rms = 0;
+  for (size_t i=1; i<vec.size(); i++) {
+    rms += vec[i]*vec[i];
+  }
+  return sqrt(rms/vec.size());
+}
+
+bool vectorLinearFit(const Vector<float> &y, float &slope, float &interceptor)
+{
+  int len_y = y.size();
+  float numerator = 0.0;
+  float denominator = 0.0;
+  float avg_x = len_y/2.0;
+  float avg_y = 0.0;
+
+  if (len_y==0)
+    return false;
+
+  for (int i=0; i<len_y; i++)
+    avg_y += y[i];
+  avg_y /= len_y;
+
+  for (int i=0; i<len_y; i++) {
+    numerator += (i - avg_x) * (y[i] - avg_y);
+    denominator += (i - avg_x) * (i - avg_x);
+  }
+
+  if (denominator == 0.0)
+    return false;
+
+  slope = numerator/denominator;
+  interceptor = avg_y - avg_x*slope;
+
+  return true;
+}
+
+void Transceiver::estimateBurstQuality(const BitVector &wBits, signalVector *received,
+                                       const GSM::Time &wTime, size_t chan,
+                                       Transceiver::BurstQuality &qual)
+{
+  signalVector *burst;
+  float slope, interceptor;
+
+  // this code supports only 4 SPS modulation
+  // we also assume that received vector is 1 SPS
+  assert(mSPSTx==4);
+
+  burst = modulateBurst(wBits, 8 + (wTime.TN() % 4 == 0), mSPSTx);
+
+  /* Debug: dump bursts to disk */
+  if (needWriteBurstToDisk(wTime, chan))
+    writeToFile(burst, wTime, chan, "_demod");
+
+  // flip values to align modulated format with the received format
+  for (size_t i=0; i<burst->size(); i++) {
+      (*burst)[i] = complex((*burst)[i].imag(), -(*burst)[i].real());
+  }
+
+  // calculate phase error for each bit
+  for (size_t i=0; i<qual.phase_err.size(); i++) {
+    float rx_phase = (*received)[i].arg();
+    // modulated data is 4SPS and is 1/4 bit shifted
+    float mod_phase = (*burst)[1+(2+i)*4].arg();
+    qual.phase_err[i] = wrapAnglePi(rx_phase - mod_phase);
+    qual.phase_err_deg[i] = rad2deg(qual.phase_err[i]);
+  }
+
+  // compensate for frequency error
+  if (vectorLinearFit(qual.phase_err, slope, interceptor)) {
+    for (size_t i=0; i<qual.phase_err.size(); i++) {
+      qual.phase_err[i] -= i*slope + interceptor;
+    }
+  }
+
+  // convert to degrees
+  for (size_t i=0; i<qual.phase_err.size(); i++) {
+    qual.phase_err_deg[i] = rad2deg(qual.phase_err[i]);
+  }
+
+  // calculate Peak and RMS values
+  qual.phase_err_max_idx = vectorMaxAbs(qual.phase_err);
+  qual.phase_err_max = qual.phase_err[qual.phase_err_max_idx];
+  qual.phase_err_rms = vectorRMS(qual.phase_err);
+
+  delete burst;
+}
+
 void Transceiver::driveReceiveFIFO(size_t chan)
 {
   SoftVector *rxBurst = NULL;
@@ -988,18 +1131,24 @@ void Transceiver::driveReceiveFIFO(size_t chan)
   double noise; // noise level in dBFS
   GSM::Time burstTime;
   bool isRssiValid; // are RSSI, noise and burstTime valid
+  Transceiver::BurstQuality qual;
 
-  rxBurst = pullRadioVector(burstTime, RSSI, isRssiValid, TOA, noise, chan);
+  rxBurst = pullRadioVector(burstTime, RSSI, isRssiValid, TOA, noise, chan, &qual);
 
   if (rxBurst) { 
     dBm = RSSI+rssiOffset;
     TOAint = (int) (TOA * 256.0 + 0.5); // round to closest integer
 
-	LOG(DEBUG) << std::fixed << std::right
+    LOG(INFO) << std::fixed << std::right
       << " time: "   << burstTime
-	  << " RSSI: "   << std::setw(5) << std::setprecision(1) << RSSI << "dBFS/" << std::setw(6) << -dBm << "dBm"
-	  << " noise: "  << std::setw(5) << std::setprecision(1) << noise << "dBFS/" << std::setw(6) << -(noise+rssiOffset) << "dBm"
-	  << " TOA: "    << std::setw(5) << std::setprecision(2) << TOA
+      << " RSSI: "   << std::setw(5) << std::setprecision(1) << RSSI << "dBFS/" << std::setw(6) << -dBm << "dBm"
+      << " noise: "  << std::setw(5) << std::setprecision(1) << noise << "dBFS/" << std::setw(6) << -(noise+rssiOffset) << "dBm"
+      << " TOA: "    << std::setw(5) << std::setprecision(2) << TOA
+      << " peak: "   << std::setw(6) << std::setprecision(1) << qual.phase_err_deg[qual.phase_err_max_idx]
+//      << " @bit "    << std::setw(3) << qual.phase_err_max_idx
+      << " RMS: "    << std::setw(6) << std::setprecision(1) << rad2deg(qual.phase_err_rms)
+//      << " bits: " << std::setw(5) << std::setprecision(1) << qual.phase_err_deg
+//      << " energy: " << std::setw(5) << std::setprecision(2) << rxBurst->getEnergy()
       << " bits: "   << *rxBurst;
 
     char burstString[gSlotLen+10];

Transceiver52M/Transceiver.h

@@ -157,6 +157,20 @@ public:
   };
 
 private:
+
+  struct BurstQuality {
+    BurstQuality()
+      // 148 bits - burst length including guard bits
+      : phase_err(148), phase_err_deg(148)
+    {}
+
+    Vector<float> phase_err;
+    Vector<float> phase_err_deg;
+    int phase_err_max_idx;
+    float phase_err_max;
+    float phase_err_rms;
+  };
+
   int mBasePort;
   std::string mAddr;
 
@@ -197,7 +211,7 @@ private:
   /** Pull and demodulate a burst from the receive FIFO */
   SoftVector *pullRadioVector(GSM::Time &wTime, double &RSSI, bool &isRssiValid,
                               double &timingOffset, double &noise,
-                              size_t chan = 0);
+                              size_t chan = 0, BurstQuality *qual = NULL);
 
   /** Set modulus for specific timeslot */
   void setModulus(size_t timeslot, size_t chan);
@@ -221,7 +235,9 @@ private:
   /** Demodulat burst and output soft bits */
   SoftVector *demodulate(TransceiverState *state,
                          signalVector &burst, complex amp,
-                         float toa, size_t tn, bool equalize);
+                         float toa, size_t tn, bool equalize,
+                         GSM::Time &wTime, size_t chan,
+                         BurstQuality *qual=NULL);
 
   int mSPSTx;                          ///< number of samples per Tx symbol
   int mSPSRx;                          ///< number of samples per Rx symbol
@@ -235,6 +251,15 @@ private:
   unsigned mMaxExpectedDelay;          ///< maximum expected time-of-arrival offset in GSM symbols
   unsigned mWriteBurstToDiskMask;      ///< debug: bitmask to indicate which timeslots to dump to disk
 
+
+  bool needWriteBurstToDisk(const GSM::Time &wTime, size_t chan)
+  {
+    /* Debug: dump bursts to disk */
+    /* bits 0-7  - chan 0 timeslots
+     * bits 8-15 - chan 1 timeslots */
+    return mWriteBurstToDiskMask & ((1<<wTime.TN()) << (8*chan));
+  }
+
   std::vector<TransceiverState> mStates;
 
   /** Start and stop I/O threads through the control socket API */
@@ -245,6 +270,10 @@ private:
   Mutex mLock;
 
 protected:
+
+  /** Estimate received burst quality and print it to debug output */
+  void estimateBurstQuality(const BitVector &wBits, signalVector *received, const GSM::Time &wTime, size_t chan, BurstQuality &qual);
+
   /** drive lower receive I/O and burst generation */
   void driveReceiveRadio();
 

Transceiver52M/sigProcLib.cpp

@@ -45,14 +45,14 @@ using namespace GSM;
 #define CLIP_THRESH		30000.0f
 
 /** Lookup tables for trigonometric approximation */
-float cosTable[TABLESIZE+1]; // add 1 element for wrap around
-float sinTable[TABLESIZE+1];
-float sincTable[TABLESIZE+1];
+double cosTable[TABLESIZE+1]; // add 1 element for wrap around
+double sinTable[TABLESIZE+1];
+double sincTable[TABLESIZE+1];
 
 /** Constants */
-static const float M_PI_F = (float)M_PI;
-static const float M_2PI_F = (float)(2.0*M_PI);
-static const float M_1_2PI_F = 1/M_2PI_F;
+static const double M_PI_F = M_PI;
+static const double M_2PI_F = (2.0 * M_PI);
+static const double M_1_2PI_F = (1.0 / M_2PI_F);
 
 /* Precomputed rotation vectors */
 static signalVector *GMSKRotationN = NULL;
@@ -270,8 +270,8 @@ complex expjLookup(float x)
 /** Library setup functions */
 void initTrigTables() {
   for (int i = 0; i < TABLESIZE+1; i++) {
-    cosTable[i] = cos(2.0*M_PI*i/TABLESIZE);
-    sinTable[i] = sin(2.0*M_PI*i/TABLESIZE);
+    cosTable[i] = cos(2.0 * M_PI * (double) i / TABLESIZE);
+    sinTable[i] = sin(2.0 * M_PI * (double) i / TABLESIZE);
   }
 }
 
@@ -281,11 +281,13 @@ void initGMSKRotationTables(int sps)
   GMSKReverseRotationN = new signalVector(157 * sps);
   signalVector::iterator rotPtr = GMSKRotationN->begin();
   signalVector::iterator revPtr = GMSKReverseRotationN->begin();
-  float phase = 0.0;
+  double phase = 0.0;
+
   while (rotPtr != GMSKRotationN->end()) {
-    *rotPtr++ = expjLookup(phase);
-    *revPtr++ = expjLookup(-phase);
-    phase += M_PI_F / 2.0F / (float) sps;
+    *rotPtr++ = complex(cos(phase), sin(phase));
+    *revPtr++ = complex(cos(-phase), sin(-phase));
+    phase += M_PI_F / 8.0;
+
   }
 
   GMSKRotation1 = new signalVector(157);
@@ -703,7 +705,7 @@ static signalVector *modulateBurstLaurent(const BitVector &bits,
   int burst_len;
   float phase;
   signalVector *c0_pulse, *c1_pulse, *c0_burst;
-  signalVector *c1_burst, *c0_shaped, *c1_shaped;
+  signalVector *c1_burst, *c0_shaped, *c1_shaped, *c2_shaped;
   signalVector::iterator c0_itr, c1_itr;
 
   /*
@@ -716,7 +718,17 @@ static signalVector *modulateBurstLaurent(const BitVector &bits,
   c0_pulse = GSMPulse->c0;
   c1_pulse = GSMPulse->c1;
 
-  burst_len = sps * (bits.size() + guard_len);
+  int i = 0, head = 4, tail = 4;
+  BitVector _bits = BitVector(148 + head + tail);
+
+  for (; i < head; i++)
+    _bits[i] = 1;
+  for (; i < 148 + head; i++)
+    _bits[i] = bits[i - head];
+  for (; i < 148 + head + tail; i++)
+    _bits[i] = 1;
+
+  burst_len = 625 + (head + tail) * sps;
 
   c0_burst = new signalVector(burst_len, c0_pulse->size());
   c0_burst->isReal(true);
@@ -727,12 +739,12 @@ static signalVector *modulateBurstLaurent(const BitVector &bits,
   c1_itr = c1_burst->begin();
 
   /* Padded differential start bits */
-  *c0_itr = 2.0 * (0x01 & 0x01) - 1.0;
+  *c0_itr = 2.0 * (0x00 & 0x01) - 1.0;
   c0_itr += sps;
 
   /* Main burst bits */
-  for (unsigned i = 0; i < bits.size(); i++) {
-    *c0_itr = 2.0 * (bits[i] & 0x01) - 1.0;
+  for (unsigned i = 0; i < _bits.size(); i++) {
+    *c0_itr = 2.0 * (_bits[i] & 0x01) - 1.0;
     c0_itr += sps;
   }
 
@@ -754,8 +766,8 @@ static signalVector *modulateBurstLaurent(const BitVector &bits,
   c1_itr += sps;
 
   /* Generate C1 phase coefficients */
-  for (unsigned i = 2; i < bits.size(); i++) {
-    phase = 2.0 * ((bits[i - 1] & 0x01) ^ (bits[i - 2] & 0x01)) - 1.0;
+  for (unsigned i = 2; i < _bits.size(); i++) {
+    phase = 2.0 * ((_bits[i - 1] & 0x01) ^ (_bits[i - 2] & 0x01)) - 1.0;
     *c1_itr = *c0_itr * Complex<float>(0, phase);
 
     c0_itr += sps;
@@ -763,13 +775,14 @@ static signalVector *modulateBurstLaurent(const BitVector &bits,
   }
 
   /* End magic */
-  int i = bits.size();
-  phase = 2.0 * ((bits[i-1] & 0x01) ^ (bits[i-2] & 0x01)) - 1.0;
+  i = _bits.size();
+  phase = 2.0 * ((_bits[i-1] & 0x01) ^ (_bits[i-2] & 0x01)) - 1.0;
   *c1_itr = *c0_itr * Complex<float>(0, phase);
 
   /* Primary (C0) and secondary (C1) pulse shaping */
   c0_shaped = convolve(c0_burst, c0_pulse, NULL, START_ONLY);
   c1_shaped = convolve(c1_burst, c1_pulse, NULL, START_ONLY);
+  c2_shaped = new signalVector(625);
 
   /* Sum shaped outputs into C0 */
   c0_itr = c0_shaped->begin();
@@ -777,11 +790,49 @@ static signalVector *modulateBurstLaurent(const BitVector &bits,
   for (unsigned i = 0; i < c0_shaped->size(); i++ )
     *c0_itr++ += *c1_itr++;
 
+  /*
+   * Generate shaping mask with squared-cosine pulse. Only 4 samples-per-symbol
+   * is supported here so use length of 8 samples or 2 symbols.
+   */
+  int len = 20;
+  float mask[len];
+  for (i = 0; i < len; i++)
+    mask[i] = 0.5 * (1.0 - cos(M_PI * (float) i / len));
+
+  /*
+   * Ramp-up mask components:
+   *     C0 filter group delay is 7.5 samples
+   *     Subtract added head bits
+   *     Subtract length of shaping mask
+   *     Delay ramp by 1 sample
+   */
+  i = 0;
+  int start = 8 + head * sps - len + 1 + 4;
+  for (;i < len; i++)
+    c2_shaped->begin()[i] = mask[i] * c0_shaped->begin()[start + i];
+
+  for (; i < 625; i++)
+    c2_shaped->begin()[i] = c0_shaped->begin()[start + i];
+
+  /*
+   * Ramp-down mask components:
+   *     Length of ramp-up mask
+   *     148 useful bits
+   */
+  int j;
+  int end = len + 148 * sps - 8;
+  assert(end + len < 625);
+  for (i = end, j = 0; i < end + len; i++, j++)
+    c2_shaped->begin()[i] *= mask[len - j - 1];
+  for (; i < 625; i++)
+    c2_shaped->begin()[i] = 0;
+
   delete c0_burst;
   delete c1_burst;
+  delete c0_shaped;
   delete c1_shaped;
 
-  return c0_shaped;
+  return c2_shaped;
 }
 
 static signalVector *modulateBurstBasic(const BitVector &bits,
@@ -1523,35 +1574,40 @@ signalVector *decimateVector(signalVector &wVector, size_t factor)
   return dec;
 }
 
-SoftVector *demodulateBurst(signalVector &rxBurst, int sps,
-                            complex channel, float TOA)
+signalVector *alignBurst(signalVector &rxBurst, complex channel, float TOA)
 {
-  signalVector *delay, *dec = NULL;
-  SoftVector *bits;
+  signalVector *delay;
 
   scaleVector(rxBurst, ((complex) 1.0) / channel);
   delay = delayVector(&rxBurst, NULL, -TOA);
 
+  return delay;
+}
+
+SoftVector *demodulateBurst(signalVector &rxBurst, int sps)
+{
+  signalVector *burst, *dec = NULL;
+  SoftVector *bits;
+
   /* Shift up by a quarter of a frequency */
-  GMSKReverseRotate(*delay, sps);
+  GMSKReverseRotate(rxBurst, sps);
 
   /* Decimate and slice */
   if (sps > 1) {
-     dec = decimateVector(*delay, sps);
-     delete delay;
-     delay = NULL;
+     dec = decimateVector(rxBurst, sps);
+     burst = dec;
   } else {
-     dec = delay;
+     burst = &rxBurst;
   }
 
-  vectorSlicer(dec);
+  vectorSlicer(burst);
 
-  bits = new SoftVector(dec->size());
+  bits = new SoftVector(burst->size());
 
   SoftVector::iterator bit_itr = bits->begin();
-  signalVector::iterator burst_itr = dec->begin();
+  signalVector::iterator burst_itr = burst->begin();
 
-  for (; burst_itr < dec->end(); burst_itr++)
+  for (; burst_itr < burst->end(); burst_itr++)
     *bit_itr++ = burst_itr->real();
 
   delete dec;

Transceiver52M/sigProcLib.h

@@ -229,16 +229,21 @@ int analyzeTrafficBurst(signalVector &rxBurst,
 signalVector *decimateVector(signalVector &wVector, size_t factor);
 
 /**
-        Demodulates a received burst using a soft-slicer.
-	@param rxBurst The burst to be demodulated.
-        @param gsmPulse The GSM pulse.
-        @param sps The number of samples per GSM symbol.
+        Applies time of arrival to align burst with bit positions
+        @param rxBurst The burst to be aligned
         @param channel The amplitude estimate of the received burst.
         @param TOA The time-of-arrival of the received burst.
+        @return The aligned burst.
+*/
+signalVector *alignBurst(signalVector &rxBurst, complex channel, float TOA);
+
+/**
+        Demodulates a received burst using a soft-slicer.
+        @param rxBurst The burst to be demodulated.
+        @param sps The number of samples per GSM symbol.
         @return The demodulated bit sequence.
 */
-SoftVector *demodulateBurst(signalVector &rxBurst, int sps,
-                            complex channel, float TOA);
+SoftVector *demodulateBurst(signalVector &rxBurst, int sps);
 
 /**
 	Design the necessary filters for a decision-feedback equalizer.