CASDuino recording functionality

Hi, I am trying to get a recording function working on the 328p.
It is all work in progress, but I am working on two tracks. A version with an opamp (so a bit of additional hardware) and a version that uses the ADC of the 328p (2 resistors and a capacitor).
I am working on the ADC version first, if that fails, I will pick up the opamp version.
There is a signal that can be detected with the ADC. Maybe we can work together?

Opamp test-circuit:

I started thinking about the recording function when I saw the Sofacas program. This uses the Goertzel algorithm to decode the data. Not knowing what that was, I started reading and found that it is a very efficient method to detect a few frequencies in a signal.I wanted to create a 'record' function for the Casduino:

• should work on an Arduino (Atmel 328p processor)
• 1200/2400bps standard MSX support
• save to sd-card
• work with or integrate with casduino
• minimum additional hardware

The next chapter in the book I found was about the sliding DFT (Digital Signal Processing by John Proakis). This is quite similar to the Goertzel, but it does not work on a block, but you update it per sample. It is an efficient method to calculate specific DFT bins. This is the method I have been experimenting with.
For the MSX, I wanted to decode both 1200 and 2400 bps signals.

For 1200 bps, the frequency bins of 1200Hz are 2400Hz are needed, for 2400 bps, 2400Hz and 4800Hz are needed. I picked a sampling frequency of 19200Hz for the Arduino. This means that one 1200Hz cycle '0' and the two cycles 2400Hz '1' fit in 16 samples. For 2400bps, this becomes 8 samples. The tape signal from the MSX was measured around 40mV pp. This is not a lot, but this is detectable with the ADC in the Arduino 328p with the reference voltage set to 1.1V.

The algorithm for the sliding DFT is the following:
Xm(n)=exp(i*2*pi*m/N)*[Xm(n-1)+x(n)-x(n-N)]

• Xm(n) is the DFT bin value of bin m at timestamp n.
• x(n) is the sample value at timestamp n.
• N is the number of samples in the window (16 for 1200bps and 8 for 2400bps at Fs=19200Hz sample frequency)
• m is the bin-number. With Fs=19200Hz and N=16, bin (m=1)=1200Hz and bin (m=2)=2400Hz

It is mentioned in the book that this is 'marginally' stable, more on that later... Since we do not have complex numbers in plain C, we have to split this into something like here below. An update for one frequency bin costs 4 multiplications and 4 additions.

#define N 16
#define PI 3.141593
int16_t xin,xout,samplebuf[N];
uint8_t sampleindex;
float tmp_re,tmp_im,bin_re,bin_im,fact_re,fact_im;

// factor (can be pre-calculated)
fact_re=cos(2*PI*m/(float)N);
fact_im=sin(2*PI*m/(float)N);

// get ADC signal
xin=ADCL | (ADCH << 8);

// swap samples in and out of buffer for sdft
xout=samplebuf[sampleindex];
samplebuf[sampleindex]=xin;
sampleindex++;
sampleindex&=(N-1);

// real additions
bin_re+=xin;
bin_re-=xout;

// complex multiply fact*f_bin, temp vars to avoid interaction
tmp_re=fact_re*bin_re-fact_im*bin_im;
tmp_im=fact_re*bin_im+fact_im*bin_re);
bin_re=tmp_re;
bin_im=tmp_im;

In order to compare (is it a '1' or a '0'?), we need the magnitude of the bins. This can be calculated: sqrt(bin_re^2+bin_im^2). Testing this (in Python on a PC) showed this worked, but now it needed to work on the Arduino. If we want to update 2 bins real-time, everything needs to fit well within 1/19200Hz =~52us. Based on the information on speed of math operations (link) it became clear that it was never going to work with floats, it might work with (16 bit) ints.(multiply ~1.4us, addition ~0.8us). So I re-wrote the routines with 16-bit 'fixed-point'. I used a 'divisor' of 32768, so this limits the values of a 16 bit signed int to [-1,1> Additions work just as normal, multiplications need to be 'adjusted' by dividing the result with the divisor value.

Next, I learned that I overlooked something significant: integer promotion in C. By definition, int16_t*int16_t returns an int16_t value. In order to get the 32 bit value, you need to cast the inputs to (int32_t). But now the entire calculation is done as a 32 bit multiply and it will be too slow... With no solution available in C, I started looking into avr assembly and found the AVR201 note from Atmel and a site discussing these multiplications (link). The 328p processor has a hardware multiplier and support for fractional multiply (fmuls). This fits perfectly with the 16 bit fixed-point integers. I modified it a bit by removing the 'mul' part of the lower bits. It is now a signed 16*16 bit fractional multiply routine that returns a 16 bit value (the upper 16 bits of the 32 bit result). This does introduce a rounding error, but it is very fast and does not need any additional shifts or division.

#define fmuls16x16_16h(result, multiplicand, multiplier) \
__asm__ __volatile__ ( \
"  clr r2 \n\t" \
" fmuls %B2, %B1 \n\t" /* (signed)ah * (signed)bh*/ \
" movw  %A0, r0 \n\t" \
" fmulsu  %B2, %A1 \n\t" /* (signed)ah * bl*/ \
" sbc %B0, r2 \n\t" \
" add %A0, r1 \n\t" \
" adc %B0, r2 \n\t" \
" fmulsu  %B1, %A2 \n\t" /* (signed)bh * al*/ \
" sbc %B0, r2 \n\t" \
" add %A0, r1 \n\t" \
" adc %B0, r2 \n\t" \
" clr r1 \n\t" \
: "=&a" (result) \
: "a" (multiplicand),  "a" (multiplier) \
: "r2" \
);

For the magnitude calculation, I approximated this by simply adding the absolute values of the real and imaginairy part. This fits in an unsigned 16 bit value 'mag':

// magnitude approximation (fails only if both values are -32768)
mag=abs(bin_re)+abs(bin_im);

Now I was finally able to test this in the Arduino only to be dissapointed... It was not stable and values were all over the place. The book mentions several options to resolve this. The one with the least impact on performance is to add a scaling factor to the largest component of the factor, real or imaginairy (fact). A factor of 0.9998 did work with the current settings. Now it was stable, but there was no signal... Boosting the signal (volume) did the trick:

//multiply for signal resolution (always safe as ADC value is 0-1023)
xin*=16;

Now I was able to detect and (manually) decode part of the signal! Blue is the 1200Hz '0' bin magnitude and in red the 2400Hz '1' magnitude value.

The next steps are a robust decoder and a routine to write the date to the sd-card. SD cards appear to have something called write latency. This means that the card will not accept data for a certain amount of time while it is internally reorganising itself. This time can be up to 250 or even 500ms depending on the card type . The only way to deal with this is to buffer the incoming data. At 19200 samples/s we quickly run out of RAM on the 328p Arduino (2k). Even if we crop the sample into a byte, the maximum buffer time is only ~100ms. So buffering the ADC values does not appear to be an option. The only alternative is to buffer the decoded data and have the complete decoding inside the interrupt routine. This would require a buffer of at least 2400bps/11*500ms=110 bytes.

This is where I am now. I am working on a robust decoder for the incoming data to properly detect the header and decode the data bytes at 1200 bps. Next steps are support for 2400bps, auto detection of 1200/2400 bps and buffered storage to the sd-card. My 'prototype code' is available here: link. Please be aware that this is not a finished product. It is very much work in progress.

Small update:
I was able to decode a 16k file from the MSX at 1200bps without a single bit error. So the decoding seems to work.
The output was via the serial port of the Arduino. Next step: writing the recorded data to an SD card.
I will update the code on github.

Nice progress.

The noise was caused by a ground loop (my mistake...) and csave, bsave and ascii seem to work.
Even 2400 bps works, although this is near the detection limit. 1200 bps is more robust. If anyone wants to give it a try, please do!
Next step, integration in casduino and see if it still fits in the 32k of flash...