fl2k_fm: add stereo and RDS support

Signed-off-by: Steve Markgraf <steve@steve-m.de>

5 files changed, 629 insertions, 11 deletions

diff --git a/include/rds_mod.h b/include/rds_mod.h
new file mode 100644
index 0000000..682e77c
--- /dev/null
+++ b/include/rds_mod.h
@@ -0,0 +1,38 @@
+/*
+ * RDS Modulator from:
+ * PiFmRds - FM/RDS transmitter for the Raspberry Pi
+ * https://github.com/ChristopheJacquet/PiFmRds
+ * 
+ * Copyright (C) 2014 by Christophe Jacquet, F8FTK
+ *
+ * adapted for use with fl2k_fm:
+ * Copyright (C) 2018 by Steve Markgraf <steve@steve-m.de>
+ *
+ * SPDX-License-Identifier: GPL-3.0+
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef RDS_H
+#define RDS_H
+
+#define RDS_MODULATOR_RATE	(57000 * 4)
+
+void get_rds_samples(double *buffer, int count);
+void set_rds_pi(uint16_t pi_code);
+void set_rds_rt(char *rt);
+void set_rds_ps(char *ps);
+void set_rds_ta(int ta);
+
+#endif /* RDS_H */
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 413348b..2255598 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -76,7 +76,7 @@ endif()
 add_executable(fl2k_file fl2k_file.c)
 add_executable(fl2k_tcp fl2k_tcp.c)
 add_executable(fl2k_test fl2k_test.c)
-add_executable(fl2k_fm fl2k_fm.c)
+add_executable(fl2k_fm fl2k_fm.c rds_waveforms.c rds_mod.c)
 set(INSTALL_TARGETS libosmo-fl2k_shared libosmo-fl2k_static fl2k_file fl2k_tcp fl2k_test fl2k_fm)
 
 target_link_libraries(fl2k_file libosmo-fl2k_shared 
diff --git a/src/fl2k_fm.c b/src/fl2k_fm.c
index 6d619dc..74953ab 100644
--- a/src/fl2k_fm.c
+++ b/src/fl2k_fm.c
@@ -30,6 +30,7 @@
 #ifndef _WIN32
 #include <unistd.h>
 #include <fcntl.h>
+#include <getopt.h>
 #else
 #include <windows.h>
 #include <io.h>
@@ -41,6 +42,7 @@
 #include <pthread.h>
 
 #include "osmo-fl2k.h"
+#include "rds_mod.h"
 
 #define BUFFER_SAMPLES_SHIFT	16
 #define BUFFER_SAMPLES		(1 << BUFFER_SAMPLES_SHIFT)
@@ -66,10 +68,15 @@ int8_t *buf2 = NULL;
 uint32_t samp_rate = 100000000;
 
 /* default signal parameters */
+#define PILOT_FREQ	19000	/* In Hz */
+#define STEREO_CARRIER	38000	/* In Hz */
+
 int delta_freq = 75000;
 int carrier_freq = 97000000;
 int carrier_per_signal;
 int input_freq = 44100;
+int stereo_flag = 0;
+int rds_flag = 0;
 
 double *freqbuf; 
 double *slopebuf; 
@@ -83,7 +90,7 @@ void usage(void)
 		"\t[-d device index (default: 0)]\n"
 		"\t[-c carrier frequency (default: 9.7 MHz)]\n"
 		"\t[-f FM deviation (default: 75000 Hz, WBFM)]\n"
-		"\t[-i input audio sample rate (default: 44100 Hz)]\n"
+		"\t[-i input audio sample rate (default: 44100 Hz for mono FM)]\n"
 		"\t[-s samplerate in Hz (default: 100 MS/s)]\n"
 		"\tfilename (use '-' to read from stdin)\n\n"
 	);
@@ -295,15 +302,16 @@ inline double modulate_sample(int lastwritepos, double lastfreq, double sample)
 	return freq;
 }
 
-void modulator(void)
+void fm_modulator_mono(int use_rds)
 {
 	unsigned int i;
 	size_t len;
 	double freq;
 	double lastfreq = carrier_freq;
-	double slope;
 	int16_t audio_buf[AUDIO_BUF_SIZE];
 	uint32_t lastwritepos = writepos;
+	double sample;
+	double rds_samples[AUDIO_BUF_SIZE];
 
 	while (!do_exit) {
 		len = writelen(AUDIO_BUF_SIZE);
@@ -313,10 +321,88 @@ void modulator(void)
 			if (len == 0)
 				do_exit = 1;
 
+			if (use_rds)
+				get_rds_samples(rds_samples, len);
+
 			for (i = 0; i < len; i++) {
+				sample = audio_buf[i] / 32767.0;
+
+				if (use_rds) {
+					sample *= 4;
+					sample += rds_samples[i];
+					sample /= 5;
+				}
+
 				/* Modulate and buffer the sample */
-				lastfreq = modulate_sample(lastwritepos, lastfreq,
-							   audio_buf[i]/32767.0);
+				lastfreq = modulate_sample(lastwritepos, lastfreq, sample);
+				lastwritepos = writepos++;
+				writepos %= BUFFER_SAMPLES;
+			}
+		} else {
+			pthread_cond_wait(&fm_cond, &fm_mutex);
+		}
+	}
+}
+
+void fm_modulator_stereo(int use_rds)
+{
+	unsigned int i;
+	size_t len, sample_cnt;
+	double freq;
+	double lastfreq = carrier_freq;
+	int16_t audio_buf[AUDIO_BUF_SIZE];
+	uint32_t lastwritepos = writepos;
+
+	dds_t pilot, stereo;
+	double L, R, LpR, LmR, sample;
+	double rds_samples[AUDIO_BUF_SIZE];
+
+	/* Prepare stereo carriers */
+	pilot = dds_init(input_freq, PILOT_FREQ, 0);
+	stereo = dds_init(input_freq, STEREO_CARRIER, 0);
+
+	while (!do_exit) {
+		len = writelen(AUDIO_BUF_SIZE);
+		if (len > 1 && !(len % 2)) {
+			len = fread(audio_buf, 2, len, file);
+
+			if (len == 0)
+				do_exit = 1;
+
+			/* stereo => two audio samples per baseband sample */
+			sample_cnt = len/2;
+
+			if (use_rds)
+				get_rds_samples(rds_samples, sample_cnt);
+
+			for (i = 0; i < sample_cnt; i++) {
+				/* Get samples for both channels, and calculate the 
+				* mono (L+R) and the difference signal used to recreate
+				* the stereo data (L-R). */
+				L = audio_buf[i*2] / 32767.0;
+				R = audio_buf[i*2+1] / 32767.0;
+				LpR = (L + R) / 2;
+				LmR = (L - R) / 2;
+
+				/* Create a composite sample consisting of the mono
+				* signal at baseband, a 19kHz pilot and a the difference
+				* signal DSB-SC modulated on a 38kHz carrier */
+				sample = 4.05 * LpR;					/* Mono signal */
+				sample += 0.9 * (dds_real(&pilot)/127.0);		/* Pilot */
+				sample += 4.05 * LmR * (dds_real(&stereo)/127.0);	/* DSB-SC stereo */
+
+				if (use_rds) {
+					/* add RDS signal */
+					sample += rds_samples[i];
+
+					/* Normalize so we get the signal within [-1, 1] */
+					sample /= 10;
+				} else {
+					sample /= 9;
+				}
+
+				lastfreq = modulate_sample(lastwritepos, lastfreq, sample);
+
 				lastwritepos = writepos++;
 				writepos %= BUFFER_SAMPLES;
 			}
@@ -346,13 +432,30 @@ int main(int argc, char **argv)
 	int dev_index = 0;
 	pthread_attr_t attr;
 	char *filename = NULL;
+	int option_index = 0;
+	int input_freq_specified = 0;
 
 #ifndef _WIN32
 	struct sigaction sigact, sigign;
 #endif
 
-	while ((opt = getopt(argc, argv, "d:c:f:i:s:")) != -1) {
+	static struct option long_options[] =
+	{
+		{"stereo", no_argument, &stereo_flag, 1},
+		{"rds",    no_argument, &rds_flag,    1},
+		{0, 0, 0, 0}
+	};
+
+	while (1) {
+		opt = getopt_long(argc, argv, "d:c:f:i:s:", long_options, &option_index);
+
+		/* end of options reached */
+		if (opt == -1)
+			break;
+
 		switch (opt) {
+		case 0:
+			break;
 		case 'd':
 			dev_index = (uint32_t)atoi(optarg);
 			break;
@@ -364,6 +467,7 @@ int main(int argc, char **argv)
 			break;
 		case 'i':
 			input_freq = (uint32_t)atof(optarg);
+			input_freq_specified = 1;
 			break;
 		case 's':
 			samp_rate = (uint32_t)atof(optarg);
@@ -384,7 +488,24 @@ int main(int argc, char **argv)
 		exit(1);
 	}
 
-	if(strcmp(filename, "-") == 0) { /* Read samples from stdin */
+	if (rds_flag && input_freq_specified) {
+		if (input_freq != RDS_MODULATOR_RATE) {
+			fprintf(stderr, "RDS modulator only works with "
+					"228 kHz audio sample rate!\n");
+			exit(1);
+		}
+	} else if (rds_flag && !input_freq_specified) {
+		input_freq = RDS_MODULATOR_RATE;
+	}
+
+	if (stereo_flag && input_freq < (RDS_MODULATOR_RATE/2)) {
+		fprintf(stderr, "Audio sample rate needs to be at least "
+				"114 kHz for stereo FM!\n");
+		exit(1);
+	}
+
+
+	if (strcmp(filename, "-") == 0) { /* Read samples from stdin */
 		file = stdin;
 #ifdef _WIN32
 		_setmode(_fileno(stdin), _O_BINARY);
@@ -452,8 +573,10 @@ int main(int argc, char **argv)
 	/* Calculate needed constants */
 	carrier_per_signal = samp_rate / input_freq;
 
-	carrier_freq = samp_rate - carrier_freq;
-
+	/* Set RDS parameters */
+	set_rds_pi(0x0dac);
+	set_rds_ps("fl2k_fm");
+	set_rds_rt("VGA FM transmitter");
 
 #ifndef _WIN32
 	sigact.sa_handler = sighandler;
@@ -468,7 +591,15 @@ int main(int argc, char **argv)
 	SetConsoleCtrlHandler( (PHANDLER_ROUTINE) sighandler, TRUE );
 #endif
 
-	modulator();
+	if (stereo_flag) {
+		fm_modulator_stereo(rds_flag);
+	} else {
+		if (rds_flag)
+			fprintf(stderr, "Warning: RDS with mono (without 19 kHz pilot"
+					" tone) doesn't work with all receivers!\n");
+
+		fm_modulator_mono(rds_flag);
+	}
 
 out:
 	fl2k_close(dev);
diff --git a/src/rds_mod.c b/src/rds_mod.c
new file mode 100644
index 0000000..865ca11
--- /dev/null
+++ b/src/rds_mod.c
@@ -0,0 +1,296 @@
+/*
+ * RDS Modulator from:
+ * PiFmRds - FM/RDS transmitter for the Raspberry Pi
+ * https://github.com/ChristopheJacquet/PiFmRds
+ * 
+ * Copyright (C) 2014 by Christophe Jacquet, F8FTK
+ *
+ * adapted for use with fl2k_fm:
+ * Copyright (C) 2018 by Steve Markgraf <steve@steve-m.de>
+ *
+ * SPDX-License-Identifier: GPL-3.0+
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdint.h>
+#include <string.h>
+#include <time.h>
+#include <stdlib.h>
+
+#define RT_LENGTH	64
+#define PS_LENGTH	8
+#define GROUP_LENGTH	4
+
+extern double waveform_biphase[576];
+
+struct {
+	uint16_t pi;
+	int ta;
+	char ps[PS_LENGTH];
+	char rt[RT_LENGTH];
+} rds_params = { 0 };
+
+/* The RDS error-detection code generator polynomial is
+   x^10 + x^8 + x^7 + x^5 + x^4 + x^3 + x^0
+*/
+#define POLY		0x1B9
+#define POLY_DEG	10
+#define MSB_BIT		(1 << 15)
+#define BLOCK_SIZE	16
+
+#define BITS_PER_GROUP (GROUP_LENGTH * (BLOCK_SIZE+POLY_DEG))
+#define SAMPLES_PER_BIT 192
+#define FILTER_SIZE (sizeof(waveform_biphase)/sizeof(double))
+#define SAMPLE_BUFFER_SIZE (SAMPLES_PER_BIT + FILTER_SIZE)
+
+
+uint16_t offset_words[] = { 0x0FC, 0x198, 0x168, 0x1B4 };
+// We don't handle offset word C' here for the sake of simplicity
+
+/* Classical CRC computation */
+uint16_t crc(uint16_t block)
+{
+	uint16_t crc = 0;
+	int i, bit, msb;
+
+	for (i = 0; i < BLOCK_SIZE; i++) {
+		bit = (block & MSB_BIT) != 0;
+		block <<= 1;
+
+		msb = (crc >> (POLY_DEG-1)) & 1;
+		crc <<= 1;
+
+		if ((msb ^ bit) != 0)
+			crc = crc ^ POLY;
+	}
+
+	return crc;
+}
+
+/* Possibly generates a CT (clock time) group if the minute has just changed
+   Returns 1 if the CT group was generated, 0 otherwise
+*/
+int get_rds_ct_group(uint16_t *blocks)
+{
+	static int latest_minutes = -1;
+	int l, mjd, offset;
+
+	// Check time
+	time_t now;
+	struct tm *utc;
+
+	now = time(NULL);
+	utc = gmtime(&now);
+
+	if(utc->tm_min != latest_minutes) {
+		// Generate CT group
+		latest_minutes = utc->tm_min;
+
+		l = utc->tm_mon <= 1 ? 1 : 0;
+		mjd = 14956 + utc->tm_mday + 
+		      (int)((utc->tm_year - l) * 365.25) +
+		      (int)((utc->tm_mon + 2 + l*12) * 30.6001);
+
+		blocks[1] = 0x4400 | (mjd>>15);
+		blocks[2] = (mjd<<1) | (utc->tm_hour>>4);
+		blocks[3] = (utc->tm_hour & 0xF)<<12 | utc->tm_min<<6;
+
+		utc = localtime(&now);
+
+		offset = utc->tm_gmtoff / (30 * 60);
+		blocks[3] |= abs(offset);
+		if (offset < 0)
+			blocks[3] |= 0x20;
+
+		return 1;
+	} else {
+		return 0;
+	}
+}
+
+/* Creates an RDS group. This generates sequences of the form 0A, 0A, 0A, 0A, 2A, etc.
+   The pattern is of length 5, the variable 'state' keeps track of where we are in the
+   pattern. 'ps_state' and 'rt_state' keep track of where we are in the PS (0A) sequence
+   or RT (2A) sequence, respectively.
+*/
+void get_rds_group(int *buffer)
+{
+	static int state = 0;
+	static int ps_state = 0;
+	static int rt_state = 0;
+	uint16_t blocks[GROUP_LENGTH] = { rds_params.pi, 0, 0, 0 };
+	uint16_t block, check;
+	int i, j;
+
+	// Generate block content
+	if (!get_rds_ct_group(blocks)) { // CT (clock time) has priority on other group types
+		if (state < 4) {
+			blocks[1] = 0x0400 | ps_state;
+
+			if (rds_params.ta)
+				blocks[1] |= 0x0010;
+
+			blocks[2] = 0xCDCD;	 // no AF
+			blocks[3] = rds_params.ps[ps_state*2] << 8 | rds_params.ps[ps_state*2+1];
+			ps_state++;
+
+			if (ps_state >= 4)
+				ps_state = 0;
+		} else { // state == 5
+			blocks[1] = 0x2400 | rt_state;
+			blocks[2] = rds_params.rt[rt_state*4+0] << 8 | rds_params.rt[rt_state*4+1];
+			blocks[3] = rds_params.rt[rt_state*4+2] << 8 | rds_params.rt[rt_state*4+3];
+			rt_state++;
+			if (rt_state >= 16)
+				rt_state = 0;
+		}
+
+		state++;
+		if (state >= 5)
+			state = 0;
+	}
+	
+	// Calculate the checkword for each block and emit the bits
+	for (i = 0; i < GROUP_LENGTH; i++) {
+		block = blocks[i];
+		check = crc(block) ^ offset_words[i];
+		for (j = 0; j < BLOCK_SIZE; j++) {
+			*buffer++ = ((block & (1 << (BLOCK_SIZE-1))) != 0);
+			block <<= 1;
+		}
+		for (j = 0; j < POLY_DEG; j++) {
+			*buffer++ = ((check & (1 << (POLY_DEG-1))) != 0);
+			check <<= 1;
+		}
+	}
+}
+
+/* Get a number of RDS samples. This generates the envelope of the waveform using
+   pre-generated elementary waveform samples, and then it amplitude-modulates the 
+   envelope with a 57 kHz carrier, which is very efficient as 57 kHz is 4 times the
+   sample frequency we are working at (228 kHz).
+ */
+void get_rds_samples(double *buffer, uint32_t count)
+{
+	static int bit_buffer[BITS_PER_GROUP];
+	static int bit_pos = BITS_PER_GROUP;
+	static double sample_buffer[SAMPLE_BUFFER_SIZE] = {0};
+
+	static int prev_output = 0;
+	static int cur_output = 0;
+	static int cur_bit = 0;
+	static int sample_count = SAMPLES_PER_BIT;
+	static int inverting = 0;
+	static int phase = 0;
+
+	static unsigned int in_sample_index = 0;
+	static unsigned int out_sample_index = SAMPLE_BUFFER_SIZE-1;
+
+	unsigned int i, j, idx;
+	double val, sample;
+	double *src;
+
+	for (i = 0; i < count; i++) {
+		if (sample_count >= SAMPLES_PER_BIT) {
+			if (bit_pos >= BITS_PER_GROUP) {
+				get_rds_group(bit_buffer);
+				bit_pos = 0;
+			}
+
+			// do differential encoding
+			cur_bit = bit_buffer[bit_pos];
+			prev_output = cur_output;
+			cur_output = prev_output ^ cur_bit;
+
+			inverting = (cur_output == 1);
+
+			src = waveform_biphase;
+			idx = in_sample_index;
+
+			for (j = 0; j < FILTER_SIZE; j++) {
+				val = *src++;
+				if (inverting)
+					val = -val;
+
+				sample_buffer[idx++] += val;
+
+				if (idx >= SAMPLE_BUFFER_SIZE)
+					idx = 0;
+			}
+
+			in_sample_index += SAMPLES_PER_BIT;
+			if (in_sample_index >= SAMPLE_BUFFER_SIZE)
+				in_sample_index -= SAMPLE_BUFFER_SIZE;
+
+			bit_pos++;
+			sample_count = 0;
+		}
+
+		sample = sample_buffer[out_sample_index];
+		sample_buffer[out_sample_index] = 0;
+		out_sample_index++;
+		if (out_sample_index >= SAMPLE_BUFFER_SIZE)
+			out_sample_index = 0;
+
+		// modulate at 57 kHz
+		// use phase for this
+		switch (phase) {
+			case 0:
+			case 2: sample = 0; break;
+			case 1: break;
+			case 3: sample = -sample; break;
+		}
+		phase++;
+		if (phase >= 4)
+			phase = 0;
+
+		*buffer++ = sample;
+		sample_count++;
+	}
+}
+
+void set_rds_pi(uint16_t pi_code)
+{
+	rds_params.pi = pi_code;
+}
+
+void set_rds_rt(char *rt)
+{
+	int i;
+
+	strncpy(rds_params.rt, rt, 64);
+
+	for (i = 0; i < 64; i++) {
+		if (rds_params.rt[i] == 0)
+			rds_params.rt[i] = 32;
+	}
+}
+
+void set_rds_ps(char *ps)
+{
+	int i;
+
+	strncpy(rds_params.ps, ps, 8);
+
+	for (i = 0; i < 8; i++) {
+		if (rds_params.ps[i] == 0)
+			rds_params.ps[i] = 32;
+	}
+}
+
+void set_rds_ta(int ta)
+{
+	rds_params.ta = ta;
+}
diff --git a/src/rds_waveforms.c b/src/rds_waveforms.c
new file mode 100644
index 0000000..d508449
--- /dev/null
+++ b/src/rds_waveforms.c
@@ -0,0 +1,153 @@
+/* This file was automatically generated by "generate_waveforms.py" from
+ * https://github.com/ChristopheJacquet/PiFmRds/
+ *
+ * (C) 2014 Christophe Jacquet.
+ * Released under the GNU GPL v3 license.
+*/
+
+double waveform_biphase[] = {
+ 0.00253265133022,  0.00255504491037,  0.00256667102126,  0.00256723854970,
+ 0.00255649674667,  0.00253423716573,  0.00250029547253,  0.00245455311551,
+ 0.00239693884806,  0.00232743009314,  0.00224605414143,  0.00215288917468,
+ 0.00204806510656,  0.00193176423352,  0.00180422168917,  0.00166572569587,
+ 0.00151661760823,  0.00135729174364,  0.00118819499588,  0.00100982622839,
+ 0.00082273544470,  0.00062752273428,  0.00042483699288,  0.00021537441720,
+-0.00000012322530, -0.00022087054977, -0.00044604072817, -0.00067476788077,
+-0.00090614968071, -0.00113925016637, -0.00137310275567, -0.00160671345499,
+-0.00183906425517, -0.00206911670572, -0.00229581565752, -0.00251809316382,
+-0.00273487252813, -0.00294507248686, -0.00314761151410, -0.00334141223473,
+-0.00352540593170, -0.00369853713255, -0.00385976825946, -0.00400808432674,
+-0.00414249766903, -0.00426205268297, -0.00436583056466, -0.00445295402495,
+-0.00452259196407, -0.00457396408696, -0.00460634544047, -0.00461907085337,
+-0.00461153926002, -0.00458321788861, -0.00453364629481, -0.00446244022186,
+-0.00436929526830, -0.00425399034471, -0.00411639090130, -0.00395645190841,
+-0.00377422057251, -0.00356983877090, -0.00334354518872, -0.00309567714275,
+-0.00282667207705, -0.00253706871632, -0.00222750786389, -0.00189873283191,
+-0.00155158949247, -0.00118702593940, -0.00080609175162, -0.00040993684994,
+ 0.00000019005938,  0.00042294345989,  0.00085688342327,  0.00130047843362,
+ 0.00175210863919,  0.00221006953169,  0.00267257605183,  0.00313776711824,
+ 0.00360371057560,  0.00406840855594,  0.00452980324611,  0.00498578305229,
+ 0.00543418915128,  0.00587282241677,  0.00629945070701,  0.00671181649912,
+ 0.00710764485348,  0.00748465169059,  0.00784055236082,  0.00817307048659,
+ 0.00847994705483,  0.00875894973632,  0.00900788240743,  0.00922459484812,
+ 0.00940699258958,  0.00955304688301,  0.00966080476069,  0.00972839915915,
+ 0.00975405907344,  0.00973611971083,  0.00967303261139,  0.00956337570235,
+ 0.00940586325271,  0.00919935569384,  0.00894286927184,  0.00863558549687,
+ 0.00827686035476,  0.00786623324593,  0.00740343561706,  0.00688839925073,
+ 0.00632126417893,  0.00570238618641,  0.00503234387065,  0.00431194522560,
+ 0.00354223371723,  0.00272449381977,  0.00186025598240,  0.00095130099727,
+-0.00000033625901, -0.00099236373728, -0.00202222980620, -0.00308712154498,
+-0.00418396376074, -0.00530941875093, -0.00645988682960, -0.00763150763453,
+-0.00882016223000, -0.01002147601830, -0.01123082247100, -0.01244332768770,
+-0.01365387579080, -0.01485711515770, -0.01604746549480, -0.01721912575140,
+-0.01836608287050, -0.01948212137170, -0.02056083375870, -0.02159563173920,
+-0.02257975824610, -0.02350630024450, -0.02436820230520, -0.02515828092590,
+-0.02586923957590, -0.02649368443980, -0.02702414083260, -0.02745307025480,
+-0.02777288805630, -0.02797598167370, -0.02805472940410, -0.02800151967640,
+-0.02780877077830, -0.02746895099680, -0.02697459912600, -0.02631834529620,
+-0.02549293207490, -0.02449123579000, -0.02330628802350, -0.02193129722250,
+-0.02035967037330, -0.01858503468320, -0.01660125921430, -0.01440247641080,
+-0.01198310346410, -0.00933786345535, -0.00646180621834, -0.00335032886310,
+ 0.00000080409844,  0.00359544108347,  0.00743702428450,  0.01152857108890,
+ 0.01587265612930,  0.02047139402060,  0.02532642284170,  0.03043888841330,
+ 0.03580942942790,  0.04143816348300,  0.04732467406730,  0.05346799854990,
+ 0.05986661721770,  0.06651844340760,  0.07342081477420,  0.08057048573450,
+ 0.08796362112670,  0.09559579111960,  0.10346196740200,  0.11155652068800,
+ 0.11987321955300,  0.12840523064300,  0.13714512025400,  0.14608485732300,
+ 0.15521581782400,  0.16452879059300,  0.17401398458500,  0.18366103756000,
+ 0.19345902621300,  0.20339647772800,  0.21346138276100,  0.22364120983400,
+ 0.23392292112600,  0.24429298964900,  0.25473741777900,  0.26524175712200,
+ 0.27579112968500,  0.28637025030900,  0.29696345035600,  0.30755470256700,
+ 0.31812764709700,  0.32866561863200,  0.33915167458400,  0.34956862427200,
+ 0.35989905906300,  0.37012538339000,  0.38022984661400,  0.39019457563000,
+ 0.40000160819100,  0.40963292682000,  0.41907049336100,  0.42829628390300,
+ 0.43729232419900,  0.44604072538100,  0.45452371993000,  0.46272369782900,
+ 0.47062324279700,  0.47820516854900,  0.48545255496600,  0.49234878413200,
+ 0.49887757611800,  0.50502302444800,  0.51076963117100,  0.51610234143400,
+ 0.52100657748800,  0.52546827205100,  0.52947390092900,  0.53301051483100,
+ 0.53606577028900,  0.53862795962200,  0.54068603984100,  0.54222966045900,
+ 0.54324919009800,  0.54373574185300,  0.54368119733200,  0.54307822931100,
+ 0.54192032294800,  0.54020179549500,  0.53791781445300,  0.53506441412200,
+ 0.53163851050000,  0.52763791447900,  0.52306134330900,  0.51790843029000,
+ 0.51217973265200,  0.50587673761100,  0.49900186656100,  0.49155847739800,
+ 0.48355086494200,  0.47498425946700,  0.46586482332000,  0.45619964562500,
+ 0.44599673508500,  0.43526501088000,  0.42401429168100,  0.41225528278300,
+ 0.39999956139700,  0.38725956008100,  0.37404854845100,  0.36038061302100,
+ 0.34627063539800,  0.33173426875400,  0.31678791267100,  0.30144868639300,
+ 0.28573440054000,  0.26966352734900,  0.25325516949400,  0.23652902755300,
+ 0.21950536619500,  0.20220497914700,  0.18464915302200,  0.16685963008900,
+ 0.14885857004900,  0.13066851092000,  0.11231232909100,  0.09381319865270,
+ 0.07519455007850,  0.05648002834940,  0.03769345061440,  0.01885876347700,
+ 0.00000000000000, -0.01885876347700, -0.03769345061440, -0.05648002834940,
+-0.07519455007850, -0.09381319865270, -0.11231232909100, -0.13066851092000,
+-0.14885857004900, -0.16685963008900, -0.18464915302200, -0.20220497914700,
+-0.21950536619500, -0.23652902755300, -0.25325516949400, -0.26966352734900,
+-0.28573440054000, -0.30144868639300, -0.31678791267100, -0.33173426875400,
+-0.34627063539800, -0.36038061302100, -0.37404854845100, -0.38725956008100,
+-0.39999956139700, -0.41225528278300, -0.42401429168100, -0.43526501088000,
+-0.44599673508500, -0.45619964562500, -0.46586482332000, -0.47498425946700,
+-0.48355086494200, -0.49155847739800, -0.49900186656100, -0.50587673761100,
+-0.51217973265200, -0.51790843029000, -0.52306134330900, -0.52763791447900,
+-0.53163851050000, -0.53506441412200, -0.53791781445300, -0.54020179549500,
+-0.54192032294800, -0.54307822931100, -0.54368119733200, -0.54373574185300,
+-0.54324919009800, -0.54222966045900, -0.54068603984100, -0.53862795962200,
+-0.53606577028900, -0.53301051483100, -0.52947390092900, -0.52546827205100,
+-0.52100657748800, -0.51610234143400, -0.51076963117100, -0.50502302444800,
+-0.49887757611800, -0.49234878413200, -0.48545255496600, -0.47820516854900,
+-0.47062324279700, -0.46272369782900, -0.45452371993000, -0.44604072538100,
+-0.43729232419900, -0.42829628390300, -0.41907049336100, -0.40963292682000,
+-0.40000160819100, -0.39019457563000, -0.38022984661400, -0.37012538339000,
+-0.35989905906300, -0.34956862427200, -0.33915167458400, -0.32866561863200,
+-0.31812764709700, -0.30755470256700, -0.29696345035600, -0.28637025030900,
+-0.27579112968500, -0.26524175712200, -0.25473741777900, -0.24429298964900,
+-0.23392292112600, -0.22364120983400, -0.21346138276100, -0.20339647772800,
+-0.19345902621300, -0.18366103756000, -0.17401398458500, -0.16452879059300,
+-0.15521581782400, -0.14608485732300, -0.13714512025400, -0.12840523064300,
+-0.11987321955300, -0.11155652068800, -0.10346196740200, -0.09559579111960,
+-0.08796362112670, -0.08057048573450, -0.07342081477420, -0.06651844340760,
+-0.05986661721770, -0.05346799854990, -0.04732467406730, -0.04143816348300,
+-0.03580942942790, -0.03043888841330, -0.02532642284170, -0.02047139402060,
+-0.01587265612930, -0.01152857108890, -0.00743702428450, -0.00359544108347,
+-0.00000080409844,  0.00335032886310,  0.00646180621834,  0.00933786345535,
+ 0.01198310346410,  0.01440247641080,  0.01660125921430,  0.01858503468320,
+ 0.02035967037330,  0.02193129722250,  0.02330628802350,  0.02449123579000,
+ 0.02549293207490,  0.02631834529620,  0.02697459912600,  0.02746895099680,
+ 0.02780877077830,  0.02800151967640,  0.02805472940410,  0.02797598167370,
+ 0.02777288805630,  0.02745307025480,  0.02702414083260,  0.02649368443980,
+ 0.02586923957590,  0.02515828092590,  0.02436820230520,  0.02350630024450,
+ 0.02257975824610,  0.02159563173920,  0.02056083375870,  0.01948212137170,
+ 0.01836608287050,  0.01721912575140,  0.01604746549480,  0.01485711515770,
+ 0.01365387579080,  0.01244332768770,  0.01123082247100,  0.01002147601830,
+ 0.00882016223000,  0.00763150763453,  0.00645988682960,  0.00530941875093,
+ 0.00418396376074,  0.00308712154498,  0.00202222980620,  0.00099236373728,
+ 0.00000033625901, -0.00095130099727, -0.00186025598240, -0.00272449381977,
+-0.00354223371723, -0.00431194522560, -0.00503234387065, -0.00570238618641,
+-0.00632126417893, -0.00688839925073, -0.00740343561706, -0.00786623324593,
+-0.00827686035476, -0.00863558549687, -0.00894286927184, -0.00919935569384,
+-0.00940586325271, -0.00956337570235, -0.00967303261139, -0.00973611971083,
+-0.00975405907344, -0.00972839915915, -0.00966080476069, -0.00955304688301,
+-0.00940699258958, -0.00922459484812, -0.00900788240743, -0.00875894973632,
+-0.00847994705483, -0.00817307048659, -0.00784055236082, -0.00748465169059,
+-0.00710764485348, -0.00671181649912, -0.00629945070701, -0.00587282241677,
+-0.00543418915128, -0.00498578305229, -0.00452980324611, -0.00406840855594,
+-0.00360371057560, -0.00313776711824, -0.00267257605183, -0.00221006953169,
+-0.00175210863919, -0.00130047843362, -0.00085688342327, -0.00042294345989,
+-0.00000019005938,  0.00040993684994,  0.00080609175162,  0.00118702593940,
+ 0.00155158949247,  0.00189873283191,  0.00222750786389,  0.00253706871632,
+ 0.00282667207705,  0.00309567714275,  0.00334354518872,  0.00356983877090,
+ 0.00377422057251,  0.00395645190841,  0.00411639090130,  0.00425399034471,
+ 0.00436929526830,  0.00446244022186,  0.00453364629481,  0.00458321788861,
+ 0.00461153926002,  0.00461907085337,  0.00460634544047,  0.00457396408696,
+ 0.00452259196407,  0.00445295402495,  0.00436583056466,  0.00426205268297,
+ 0.00414249766903,  0.00400808432674,  0.00385976825946,  0.00369853713255,
+ 0.00352540593170,  0.00334141223473,  0.00314761151410,  0.00294507248686,
+ 0.00273487252813,  0.00251809316382,  0.00229581565752,  0.00206911670572,
+ 0.00183906425517,  0.00160671345499,  0.00137310275567,  0.00113925016637,
+ 0.00090614968071,  0.00067476788077,  0.00044604072817,  0.00022087054977,
+ 0.00000012322530, -0.00021537441720, -0.00042483699288, -0.00062752273428,
+-0.00082273544470, -0.00100982622839, -0.00118819499588, -0.00135729174364,
+-0.00151661760823, -0.00166572569587, -0.00180422168917, -0.00193176423352,
+-0.00204806510656, -0.00215288917468, -0.00224605414143, -0.00232743009314,
+-0.00239693884806, -0.00245455311551, -0.00250029547253, -0.00253423716573,
+-0.00255649674667, -0.00256723854970, -0.00256667102126, -0.00255504491037,
+};