An Integer Version of the Demodulator

As the code will be ported to a lowly Arduino using C code, I first looked at if I could use 16 bit signed integers. The answer is strictly yes. Strictly in that the 32 bit result of a 16 bit multiplication is not used for the following division. For assembler you can take advantage of the 32 bit result of a multiplication and division, but "C" does not. This strict requirement limits the variable resolution to +127 and -128 as the high byte is used as intermediate overflow before a resoring division. In simple terms I am using fixed point arithmetic.

Code Optimisation

Note that the code only uses the high 8 bits of the ADC value even though the ADC is 10 bits.

The main drawback with using integer variables is the small signal range. The basic code can detect signals with peak to peak (pp) voltages between 380 mv pp and 5v pp. The correlation function amplitude is a square of signal amplitude.

The demodulator is insensitive to phase and gaussian noise. The demodulator is however very sensitive to harmonic distortion such as clipping and to a lesser extent to attempts to compress the input signal amplitude range (I did try).

The solution is to take the square root of the correlation. Then the correlation function amplitude is linear to the signal amplitude. This extends the detection signal range down to 78 mv pp (i.e. 2 bits).

Different low pass filters were tested with a four pole 900 Hz (-3 dB) determined to work the best.

For maximum dynamic signal range and considering some noise, a signal of 1 v pp should be targeted.

Here is the Excel code:

Option Explicit
Sub AFSKDem()
Application.ScreenUpdating = False
Application.Calculation = xlCalculationManual
Application.EnableEvents = False

Dim Fspace As Integer
Dim Fmark As Integer
Dim Fsample As Integer
Dim Baud As Integer
Dim Ticks As Integer

Dim I As Integer
Dim Xroot As Integer

Dim Time As Integer
Dim Code As Integer
Dim DataIn As Integer
Dim Phase As Integer
Dim Signal As Integer
Dim Correl As Integer
Dim DataOut As Integer

Dim X2 As Integer
Dim X1 As Integer
Dim X0 As Integer
Dim Y2 As Integer
Dim Y1 As Integer
Dim Y0 As Integer
Dim Z2 As Integer
Dim Z1 As Integer
Dim Z0 As Integer
Dim S0 As Integer
Dim S1 As Integer
Dim S2 As Integer
Dim S3 As Integer
Dim S4 As Integer
Dim S5 As Integer
Dim S6 As Integer
Dim iSin(65) As Integer

' AFSK parameters
Fspace = 1200
Fmark = 2200
Fsample = 13200
Baud = 1200
Ticks = Fsample / Fmark * Fsample / Fspace

' Title
ActiveSheet.Cells(1, 1).Value = "iBell 202"
ActiveSheet.Cells(2, 1).Value = "Baud"
ActiveSheet.Cells(3, 1).Value = "Fspace"
ActiveSheet.Cells(4, 1).Value = "Fmark"
ActiveSheet.Cells(5, 1).Value = "Fsample"
ActiveSheet.Cells(6, 1).Value = "Ticks"
ActiveSheet.Cells(7, 1).Value = "Delay"
ActiveSheet.Cells(1, 2).Value = "i900Hz"
ActiveSheet.Cells(2, 2).Value = Baud
ActiveSheet.Cells(3, 2).Value = Fspace
ActiveSheet.Cells(4, 2).Value = Fmark
ActiveSheet.Cells(5, 2).Value = Fsample
ActiveSheet.Cells(6, 2).Value = Ticks
ActiveSheet.Cells(7, 2).Value = 6

' Integer Sin() table
iSin(0) = 127: iSin(11) = 237: iSin(22) = 237: iSin(33) = 127: iSin(44) = 17: iSin(55) = 17
iSin(1) = 139: iSin(12) = 243: iSin(23) = 231: iSin(34) = 115: iSin(45) = 11: iSin(56) = 23
iSin(2) = 151: iSin(13) = 247: iSin(24) = 223: iSin(35) = 103: iSin(46) = 7: iSin(57) = 31
iSin(3) = 163: iSin(14) = 251: iSin(25) = 215: iSin(36) = 91: iSin(47) = 3: iSin(58) = 39
iSin(4) = 174: iSin(15) = 253: iSin(26) = 206: iSin(37) = 80: iSin(48) = 1: iSin(59) = 48
iSin(5) = 185: iSin(16) = 254: iSin(27) = 196: iSin(38) = 69: iSin(49) = 0: iSin(60) = 58
iSin(6) = 196: iSin(17) = 254: iSin(28) = 185: iSin(39) = 58: iSin(50) = 0: iSin(61) = 69
iSin(7) = 206: iSin(18) = 253: iSin(29) = 174: iSin(40) = 48: iSin(51) = 1: iSin(62) = 80
iSin(8) = 215: iSin(19) = 251: iSin(30) = 163: iSin(41) = 39: iSin(52) = 3: iSin(63) = 91
iSin(9) = 223: iSin(20) = 247: iSin(31) = 151: iSin(42) = 31: iSin(53) = 7: iSin(64) = 103
iSin(10) = 231: iSin(21) = 243: iSin(32) = 139: iSin(43) = 23: iSin(54) = 11: iSin(65) = 115

' Generate some data to encode
Time = 0
Phase = 0
Signal = 0
Correl = 0
DataOut = 0
' Shift registers
S0 = 0
S1 = 0
S2 = 0
S3 = 0
S4 = 0
S5 = 0
S6 = 0
X0 = 0
X1 = 0
X2 = 0
Y0 = 0
Y1 = 0
Y2 = 0
Z0 = 0
Z1 = 0
Z2 = 0
For Code = 0 To 9
  ' DataIn = Int(Rnd + 0.5)
  ' DataIn = Code Mod 2
  If (Code = 0) Then DataIn = 0
  If (Code = 1) Then DataIn = 1
  If (Code = 2) Then DataIn = 0
  If (Code = 3) Then DataIn = 1
  If (Code = 4) Then DataIn = 1
  If (Code = 5) Then DataIn = 0
  If (Code = 6) Then DataIn = 0
  If (Code = 7) Then DataIn = 1
  If (Code = 8) Then DataIn = 1
  If (Code = 9) Then DataIn = 0
  
  ' Encode
  For I = 1 To Fsample / Baud

    If (DataIn = 0) Then
      ' Fspace
      Signal = iSin(Phase)-128
      Phase = Phase + Fsample / Fmark
      If (Phase >= Ticks) Then Phase = Phase - Ticks
    Else
      ' Fmark
      Signal = iSin(Phase)-128
      Phase = Phase + Fsample / Fspace
      If (Phase >= Ticks) Then Phase = Phase - Ticks
    End If

    ' Decode the data
    Z2 = Z1
    Z1 = Z0
    Y2 = Y1
    Y1 = Y0
    X2 = X1
    X1 = X0
    S6 = S5
    S5 = S4
    S4 = S3
    S3 = S2
    S2 = S1
    S1 = S0
    S0 = Signal
    X0 = S0 * S6
    ' Integer square root
    If (X0 > 1) Then
      Xroot = 17 + X0 / 130
      Xroot = (Xroot + X0 / Xroot) / 2
      Xroot = (Xroot + X0 / Xroot) / 2
      Xroot = (Xroot + X0 / Xroot) / 2
      X0 = Xroot
    ElseIf (X0 < -1) Then
      Xroot = 17 - X0 / 130
      Xroot = (Xroot - X0 / Xroot) / 2
      Xroot = (Xroot - X0 / Xroot) / 2
      Xroot = (Xroot - X0 / Xroot) / 2
      X0 = -Xroot
    End If
    
    ' 600 Hz (-3dB) four pole critically damped low pass filter using approximate fractions
    ' Y0 = X0 * 40 / 573 + X1 * 80 / 573 + X2 * 40 / 589 + Y1 * 83 / 88 - Y2 * 2 / 9
    ' Z0 = Y0 * 40 / 573 + Y1 * 80 / 573 + Y2 * 40 / 589 + Z1 * 83 / 88 - Z2 * 2 / 9
    
    ' 800 Hz (-3dB) four pole critically damped low pass filter using approximate fractions
    ' Y0 = X0 * 91 / 830 + X1 * 91 / 415 + X2 * 91 / 830 + Y1 * 127 / 188 - Y2 * 81 / 710
    ' Z0 = Y0 * 91 / 830 + Y1 * 91 / 415 + Y2 * 91 / 830 + Z1 * 127 / 188 - Z2 * 81 / 710
    
    ' 900 Hz (-3dB) four pole critically damped low pass filter using approximate fractions
    Y0 = X0 * 51 / 388 + X1 * 51 / 194 + X2 * 51 / 388 + Y1 * 11 / 20 - Y2 * 43 / 569
    Z0 = Y0 * 51 / 388 + Y1 * 51 / 194 + Y2 * 51 / 388 + Z1 * 11 / 20 - Z2 * 43 / 569
    
    ' 1000 Hz (-3dB) four pole critically damped low pass filter using approximate fractions
    ' Y0 = X0 * 59 / 382 + X1 * 59 / 191 + X2 * 59 / 382 + Y1 * 107 / 250 - Y2 * 49 / 1070
    ' Z0 = Y0 * 59 / 382 + Y1 * 59 / 191 + Y2 * 59 / 382 + Z1 * 107 / 250 - Z2 * 49 / 1070
    
    Correl = Z0
    If (Correl > 0) Then
      DataOut = 1
    Else
      DataOut = 0
    End If
    
    ActiveSheet.Cells(10, 1).Value = "Time"
    ActiveSheet.Cells(10, 2).Value = "DataIn"
    ActiveSheet.Cells(10, 3).Value = "Signal"
    ActiveSheet.Cells(10, 4).Value = "Correlation"
    ActiveSheet.Cells(10, 5).Value = "DataOut"
    ActiveSheet.Cells(11 + Time, 1).Value = Time
    ActiveSheet.Cells(11 + Time, 2).Value = DataIn
    ActiveSheet.Cells(11 + Time, 3).Value = Signal / 128#
    ActiveSheet.Cells(11 + Time, 4).Value = Correl / 128#
    ActiveSheet.Cells(11 + Time, 5).Value = DataOut - 3
    
    Time = Time + 1
  Next I
Next Code

Application.ScreenUpdating = True
Application.Calculation = xlCalculationAutomatic
Application.EnableEvents = True
End Sub

And here is the results:

Arduino Code

Analog Input Conversion Time

The ADC conversion (i.e. analogRead()) needs to speed up from 112 us (104 us theoretic minimum conversion time). This can be done (at the cost of reduced accuracy) by increasing the clock from 125 kHz to 1 MHz. The conversion time is now 17.2 us and the accuracy is still at least 8 bits.

The code to do this is:

// Define various ADC prescaler
const unsigned char PS_16=(1<<ADPS2);
const unsigned char PS_32=(1<<ADPS2)|(1<<ADPS0);
const unsigned char PS_64=(1<<ADPS2)|(1<<ADPS1);
const unsigned char PS_128=(1<<ADPS2)|(1<<ADPS1)|(1<<ADPS0);

  // Set ADC prescaler 
  ADCSRA&=~PS_128;    // Remove bits set by Arduino library
  // ADCSRA|=PS_64;   // 64 prescaler (250 kHz assuming a 16MHz clock)
  // ADCSRA|=PS_32;   // 32 prescaler (500 kHz assuming a 16MHz clock)
  ADCSRA|=PS_16;      // 16 prescaler (1 MHz assuming a 16MHz clock)

Generating High Quality Audio

To do this we need two timers, one to generate the sample ticks at 13200 Hz and one to generate a fast PWM. Here is some code that I spent most of the day working on:

// Bell 202 modem
#define Fspace     1200
#define Fmark      2200
#define Fsample   13200
#define Baud       1200
#define sampleTicks  66     // sampleTicks = Fsample / Fmark * Fsample / Fspace
#define baudTicks    44     // baudTicks = Fsample / Baud

// Sine table for Fsample equal to 13200 Hz 
#include <avr/pgmspace.h> 
const byte iSin[sampleTicks] PROGMEM = {
  128,140,152,164,175,186,197,207,216,224,232,
  238,244,248,252,254,255,255,254,252,248,244,
  238,232,224,216,207,197,186,175,164,152,140,
  128,116,104, 92, 81, 70, 59, 49, 40, 32, 24,
   18, 12,  8,  4,  2,  1,  1,  2,  4,  8, 12,
   18, 24, 32, 40, 49, 59, 70, 81, 92,104,116,
};
 
// Define various ADC prescaler
const unsigned char PS_16=(1<<ADPS2);
const unsigned char PS_32=(1<<ADPS2)|(1<<ADPS0);
const unsigned char PS_64=(1<<ADPS2)|(1<<ADPS1);
const unsigned char PS_128=(1<<ADPS2)|(1<<ADPS1)|(1<<ADPS0);

void setup() {
  // Set ADC prescaler (assume a 16 MHz clock)
  ADCSRA&=~PS_128;                                 // Remove bits set by Arduino library
  // ADCSRA|=PS_64;                                // 64 prescaler (250 kHz)
  // ADCSRA|=PS_32;                                // 32 prescaler (500 kHz)
  ADCSRA|=PS_16;                                   // 16 prescaler (1 MHz)
  // About 17.2 us to return an ADC value but only good to 8 bits
  
  cli();                                           // Disable interrupts
  // ATmega48A/PA/88A/PA/168A/PA/328/P
  // Use Timer 0 for sample rate (13.2 kHz)
  // Note delay() and millis() will no longer work (use microseconds() instead) 
  TIMSK0 = 0;                                      // Timer interrupts OFF
  TCCR0A = (2 << WGM00);                           // CTC mode
  TCCR0B = (2 << CS00);                            // prescaler 8, 2 MHz clock
  TIMSK0 = (1 << OCIE0A);                          // COMPA interrupt
  OCR0A = 151;                                     // Sample rate: 2 MHz/152 = 13.2 kHz
  
  // Use Timer 2 for Audio PWM
  // ATmega48A/PA/88A/PA/168A/PA/328/P
  // Note tone library will no longer work
  TIMSK2 = 0;                                      // Timer interrupts OFF
  TCCR2A = (2 << COM2A0)|(3 << WGM20);             // Fast PWM mode, toggle output on OC2A (PB3/D11)
  TCCR2B = (0 << WGM22)|(1 << CS20);               // 16 MHz clock
  OCR2A = 128;                                     // Duty cycle set in OCR2A  (128 = 0v) 
  sei();                                           // Enable interrupts 

  pinMode(11,OUTPUT);                              // Fast PWM
  pinMode(LED_BUILTIN,OUTPUT);                          
}


volatile byte Phase=0;
volatile byte Timer=0;
volatile byte DataPtr=0;
ISR(TIMER0_COMPA_vect) {
  OCR2A = pgm_read_byte(&(iSin[Phase]));
  Phase+=6;
  if (Phase>=sampleTicks) Phase-=sampleTicks;
  Timer++;
  if (Timer>=baudTicks) {
    Timer=0;
    DataPtr++;
  }
}

void loop() {
  
}

Currently the code generates a nice steady sine wave. Next is to add the data stream.

The PWM Filter

I had a play with a 3rd order Sallen Key:

This particular filter is -3 dB down on both the Mark and Space frequencies and -80dB down on the 31.25kHz PWM frequency. I will need to use an OpAmp on the input side anyway so a single supply LM324 should work okay.

I thought about using a transistor for the OpAmp but the gain is insufficient except to act as a buffer for three sequential RC low pass filters (no fancy Chebyshev filters here!).

Demodulator Frequency Response

Before moving on I checked the input frequency response by hacking the code. Although not a proper analysis the results appear realistic:

So the Space frequency response is pretty sharp (bandwidth about 550 Hz). The demodulator is insensitive to the bandwidth below about 1000 Hz which is probably why the Bell 202 protocol can use 150 baud and 387/487 Hz FSK as the return communications channel.

The Mark frequency response is open above 2200 Hz. A 2400 Hz low pass filter (in software?) may be of merit if high frequency noise may be present.

Demondulator Simplification

Could it really be simplified even more? Well yes, thinking about the radio mixer analogy I realised I could achieve a similar to the square and square root just by using the sign of the delayed signal (basically a switching mixer) . Here is the demodulator code:

    ' Decode the data
    S6 = S5
    S5 = S4
    S4 = S3
    S3 = S2
    S2 = S1
    S1 = S0
    S0 = Signal
    Z2 = Z1
    Z1 = Z0
    Y2 = Y1
    Y1 = Y0
    X2 = X1
    X1 = X0
    X0 = 0
    If (S6 > 0) Then X0 = S0
    If (S6 < 0) Then X0 = -S0

    ' 900 Hz (-3dB) four pole critically damped low pass filter using approximate fractions
    Y0 = X0 * 51 / 388 + X1 * 51 / 194 + X2 * 51 / 388 + Y1 * 11 / 20 - Y2 * 43 / 569
    Z0 = Y0 * 51 / 388 + Y1 * 51 / 194 + Y2 * 51 / 388 + Z1 * 11 / 20 - Z2 * 43 / 569
    
    DataOut = 0
    If (Z0 > 0) Then DataOut = 1

Yes, the four pole filter is the most complicated bit.

Now the benefit of this code is that the full ADC precision (10 bits) can be used without integer overflow.

Here is the results (including some corrections to the low pass filter):

---

Need More Speed

Even with my speed ups the interrupt service request (ISR) time estimate was 75 us. At 13200 Hz the total available time is 75 us, so I need to do better. The low pass filter cost about 25 us per stage so I dropped it back to one stage and lowered the cut off frequency. An few hour later I managed to get the Nano to generate a test pattern, feed it through an external two stage RC filter and then back again for decoding. Success, here is the screen shot:

There is about a 50 us second delay between the input signal (yellow) and the digital output (blue). The input signal mark frequency is attenuated due to the simple two stage RC filter.

I am pretty pleased with this!

Here is the Arduino code:

// Bell 202 modem
#define Fspace     1200
#define Fmark      2200
#define Fsample   13200
#define Baud       1200
#define sampleTicks  66     // sampleTicks = Fsample / Fmark * Fsample / Fspace

// Sine table for Fsample equal to 13200 Hz 
const byte iSin[sampleTicks] = {
  128,140,152,164,175,186,197,207,216,224,232,
  238,244,248,252,254,255,255,254,252,248,244,
  238,232,224,216,207,197,186,175,164,152,140,
  128,116,104, 92, 81, 70, 59, 49, 40, 32, 24,
   18, 12,  8,  4,  2,  1,  1,  2,  4,  8, 12,
   18, 24, 32, 40, 49, 59, 70, 81, 92,104,116,
};
 
// Define various ADC prescaler
const byte PS_16=(1<<ADPS2);
const byte PS_32=(1<<ADPS2)|(1<<ADPS0);
const byte PS_64=(1<<ADPS2)|(1<<ADPS1);
const byte PS_128=(1<<ADPS2)|(1<<ADPS1)|(1<<ADPS0);

void setup() {
  pinMode(11,OUTPUT);                              // Fast PWM
  pinMode(12,OUTPUT);                              // DataOut
  pinMode(A0,INPUT);                               // DataIn
  pinMode(LED_BUILTIN,OUTPUT);                     // Test loop speed

  // Set ADC prescaler (assume a 16 MHz clock)
  ADCSRA&=~PS_128;                                 // Remove bits set by Arduino library
  ADCSRA|=PS_16;                                   // 16 prescaler (1 MHz)
  // About 17.2 us to return an ADC value but only good for 8 bits
  
  cli();                                           // Disable interrupts
  // ATmega48A/PA/88A/PA/168A/PA/328/P
  // Use Timer 0 for sample rate (13.2 kHz)
  // Note delay() and millis() will no longer work (use microseconds() instead) 
  TIMSK0 = 0;                                      // Timer interrupts OFF
  TCCR0A = (2 << WGM00);                           // CTC mode
  TCCR0B = (2 << CS00);                            // prescaler 8, 2 MHz clock
  TIMSK0 = (1 << OCIE0A);                          // COMPA interrupt
  OCR0A = 151;                                     // Sample rate: 2 MHz/152 = 13.2 kHz
  
  // Use Timer 2 for Audio PWM
  // ATmega48A/PA/88A/PA/168A/PA/328/P
  TIMSK2 = 0;                                      // Timer interrupts OFF
  TCCR2A = (2 << COM2A0)|(3 << WGM20);             // Fast PWM mode, toggle output on OC2A (PB3/D11)
  TCCR2B = (0 << WGM22)|(1 << CS20);               // 16 MHz clock
  OCR2A = 128;                                     // Duty cycle set in OCR2A  (128 = 0v) 
  sei();                                           // Enable interrupts 
}


volatile byte DataIn=0;
volatile byte DataOut=0;
volatile byte baudTick=11;
ISR(TIMER0_COMPA_vect) {
  static byte Phase=0;
  static int S0=0,S1=0,S2=0,S3=0,S4=0,S5=0,S6=0;
  static int X0=0,X1=0,X2=0;
  static int Y0=0,Y1=0,Y2=0;
  
  // Timer 2 controls PWM output
  OCR2A = iSin[Phase];
  if (DataIn==0) {
    Phase+=6;
  } else {
    Phase+=11;
  }
  if (Phase>=sampleTicks) Phase-=sampleTicks;
  // Done when baudTicks==0
  if (baudTick>0) baudTick--;
  
  // Decode the data
  S6=S5;S5=S4;S4=S3;S3=S2;S2=S1;S1=S0;
  X2=X1;X1=X0;
  Y2=Y1;Y1=Y0;
  S0=analogRead(A0)-512;                  // 18 us
  X0=(S6>0)?S0:-S0;
  // 600 Hz low pass filter
  Y0=(2*(X0+X1+X1+X2)+Y1*74-Y2*23)/59;    // 23 us
  DataOut=(Y0>0)?1:0;
}

void loop() {
  static byte lastBaudTick=0;
  static byte temp=0;

  if (lastBaudTick!=baudTick) {
    digitalWrite(12,DataOut);
    if (baudTick==0) {
      baudTick=11;
      DataIn=temp;
    } else if (baudTick==1) {
      temp=1-temp;
    }    
    lastBaudTick=baudTick;
    // Output loop toggle: should be almost 6600 Hz and steady
    digitalWrite(13,!digitalRead(13)); 
  }
}

It should not worry the low pass filter too much if I replace it with:

• Y0=(X0+X1+X1+X2+Y1*40-Y2*12)>>5; // 8 us

I think the next step to to interface a speaker and a microphone to the Ardunio!

After that I want to consider the 150 baud answer channel.

AlanX

Internet Research

Protocols

After a day or so of researching the Internet it seems as if there are two protocols I should look at:

1. Bell 103A: 300 baud full duplex using 1070/1270 Hz and 2070/2270 Hz
2. Bell 202: 1200 baud half duplex or 1200/150 baud full duplex using 1200/2200 Hz and 387/487 Hz

Regardless of the protocol the possible encoding/decoding technologies are:

• Separate frequency filters and envelop detectors and NRZ using a comparator
• PLL and/or tone decoders (lock time likely too slow for most protocols)
• Zero crossing and period/frequency measurement
• DFT
• Goertzel
• Cross-correlation
• and more.

If your familiar with radio receiver technology much of the above is rather similar.

Cross-correlations seems to be the go!

This article by Denis Sequine simplifies the process down to basics:

https://dx.eng.uiowa.edu/eedesign/fskcorr.pdf

Another requirement is digital filter construction and I found this site to be useful;

http://unicorn.us.com/alex/allpolefilters.html

My First Software Encode and Decode for the Bell 202 Protocol

I often prototype software using Excel Macros. Its quick and easy to debug.

So the first task was to encode a binary sequence.

The first issue was how to join the two frequencies up as 2200 Hz does not match the the 1200 baud rate? The answer is to keep track of the phase.

Here is an example of the encoded output:

Note that the output has been smoothed by excel. The article by Denis Sequine suggests using a 13.2 kHz sample clock, so I went with that. There are only 11 sample points per bit.

Next is to read back the data and calculate the cross-correlation for different delay steps to find the best case. The best case was a delay of 6. Refer to Denis Sequine article as to what this all means.

Finally a low pass filter (green) and a schmitt trigger was used to recover the data (bright red). I found a four pole low pass with a corner frequency of 900 Hz worked well.

Here is an example of the decoded output (note the decoded output is shifted to the right approximately one bit:

Noise

Of more interest is the effect of phase, amplitude and general noise on the signal recovery:

Now I don't know about you but I am impressed by how well this simple code works!

Here is the code (it still needs work for all the parameters to be effective):

Option Explicit
Sub AFSKDem()
Application.ScreenUpdating = False
Application.Calculation = xlCalculationManual
Application.EnableEvents = False

Const Pi As Double = 3.14159265358979
Dim Fspace As Double
Dim Fmark As Double
Dim Fsample As Double
Dim Baud As Double
Dim Delay As Long
Dim I As Long
Dim Time As Long
Dim Code As Long
Dim Data As Long
Dim Phase As Double
Dim Magn As Double
Dim MagnN As Double
Dim TestX2 As Double
Dim TestX1 As Double
Dim TestX0 As Double
Dim TestY2 As Double
Dim TestY1 As Double
Dim TestY0 As Double
Dim TestZ2 As Double
Dim TestZ1 As Double
Dim TestZ0 As Double

' Get FSK parameters
Fspace = ActiveSheet.Cells(1, 2).Value
Fmark = ActiveSheet.Cells(2, 2).Value
Fsample = ActiveSheet.Cells(3, 2).Value
Baud = ActiveSheet.Cells(4, 2).Value
Delay = ActiveSheet.Cells(5, 2).Value

' Generate some data to encode
Time = 0
Magn = 0#
Phase = 0#
ActiveSheet.Cells(7, 1).Value = "Time"
ActiveSheet.Cells(7, 2).Value = "Code"
ActiveSheet.Cells(7, 3).Value = "Data"
For Code = 0 To 7
  'Data = Int(Rnd + 0.5)
  Data = Code Mod 2
  
  ' Encode
  For I = 1 To Int(Fsample / Baud)
    ActiveSheet.Cells(8 + Time, 1).Value = Time
    ActiveSheet.Cells(8 + Time, 2).Value = Data
    ActiveSheet.Cells(8 + Time, 3).Value = Magn
    If (Data = 0) Then
      ' Fspace
      Magn = Sin(2 * Pi * (Phase + (Rnd - 0.5) * 0.2)) * (0.5 + Rnd) + Rnd * Rnd * (Rnd - 0.5)
      Phase = Phase + Fspace / Fsample
    Else
      ' Fmark
      Magn = Sin(2 * Pi * (Phase + (Rnd - 0.5) * 0.2)) * (0.5 + Rnd) + Rnd * Rnd * (Rnd - 0.5)
      Phase = Phase + Fmark / Fsample
    End If
    Time = Time + 1
  Next I
Next Code

' Decode the data
ActiveSheet.Cells(7, 4).Value = "Test"
Time = 0
Magn = 0#
MagnN = 0#
TestX0 = 0#
TestX1 = 0#
TestX2 = 0#
TestY0 = 0#
TestY1 = 0#
TestY2 = 0#
TestZ0 = 0#
TestZ1 = 0#
TestZ2 = 0#

For Code = 0 To 7
  For I = 1 To Int(Fsample / Baud)
    Magn = ActiveSheet.Cells(8 + Time, 3).Value
    If (Time >= Delay) Then
      MagnN = ActiveSheet.Cells(8 + Time - Delay, 3).Value
    End If
    TestZ2 = TestZ1
    TestZ1 = TestZ0
    TestY2 = TestY1
    TestY1 = TestY0
    TestX2 = TestX1
    TestX1 = TestX0
    TestX0 = Magn * MagnN
    ' 900 Hz Four Pole Low Pass Filter using approximate fractions
    TestY0 = TestX0 * 110 / 1667 + TestX1 * 123 / 932 + TestX2 * 110 / 1667 + TestY1 * 106 / 109 - TestY2 * 61 / 258
    TestZ0 = TestY0 * 110 / 1667 + TestY1 * 123 / 932 + TestY2 * 110 / 1667 + TestZ1 * 106 / 109 - TestZ2 * 61 / 258
    ActiveSheet.Cells(8 + Time, 4).Value = TestZ0
    
    Time = Time + 1
  Next I
Next Code

Application.ScreenUpdating = True
Application.Calculation = xlCalculationAutomatic
Application.EnableEvents = True
End Sub

AlanX