#include #include #include #include #include #include #include #include "quisk.h" #include #include "microphone.h" #include "filter.h" #ifdef MS_WINDOWS #include static int mic_cleanup = 0; // must clean up winsock #else #include #include #include #endif #define DEBUG 0 #define DEBUG_LEVEL 0 #if DEBUG_LEVEL || DEBUG_IO || DEBUG static int debug_timer = 1; // count up number of samples #endif // The microphone samples must be 48000 sps or 8000 sps. The output sample // rate is always MIC_OUT_RATE samples per second // FM needs pre-emphasis and de-emphasis. See vk1od.net/FM/FM.htm for details. // For IIR design, see http://www.abvolt.com/research/publications2.htm. // Microhone preemphasis: boost high frequencies 0.00 to 1.00 double quisk_mic_preemphasis; // Microphone clipping; try 3.0 or 4.0 double quisk_mic_clip; // If true, decimate 48000 sps mic samples to 8000 sps for processing #define DECIM_8000 1 struct alc { complex double * buffer; int buf_size; int index; int block_index; int counter; int fault; double max_magn; double gain_now[20]; double gain_max; double gain_min; double gain_change; double next_change; double final_gain; } ; // These are external: int mic_max_display; // display value of maximum microphone signal level 0 to 2**15 - 1 int quiskSpotLevel = -1; // level is -1 for Spot button Off; else the Spot level 0 to 1000. int quiskImdLevel = 500; // level for rxMode IMD, 0 to 1000 int hermes_mox_bit; // mox bit to send to Hermes static SOCKET mic_socket = INVALID_SOCKET; // send microphone samples to a socket static double mic_agc_level = 0.10; // Mic levels below this are noise and are ignored static int mic_level; // maximum microphone signal level for display static int mic_timer; // time to display maximum mic level static int align4; // add two bytes to start of audio samples to align to 4 bytes static double modulation_index = 1.6; // For FM transmit, the modulation index static int is_vox = 0; // Is the VOX level exceeded? static int vox_level = CLIP16; // VOX trigger level as a number 0 to CLIP16 static int timeVOX = 2000; // VOX hang time in milliseconds static int tx_sample_rate = 48000; // Used for SoapySDR static int reverse_tx_sideband; static int PsEnable = 0; static int PsCal = 0; static int doTxCorrect = 0; // Corrections for UDP sample transmit static double TxCorrectLevel; static complex TxCorrectDc; // Used for the Hermes protocol #define HERMES_TX_BUF_SAMPLES 4800 // buffer size in I/Q samples (two shorts) #define HERMES_TX_BUF_SHORTS (HERMES_TX_BUF_SAMPLES * 2) static int hermes_read_index; // index to read from buffer static int hermes_write_index; // index to write to buffer static int hermes_num_samples; // number of samples in the buffer static short hermes_buf[HERMES_TX_BUF_SHORTS]; // buffer to store Tx I/Q samples waiting to be sent at 48 ksps static int hermes_filter_rx; // hermes filter to use for Rx static int hermes_filter_tx; // hermes filter to use for Tx static void serial_key_samples(complex double *, int); static play_state_t last_play_state; #define TX_BLOCK_SHORTS 600 // transmit UDP packet with this many shorts (two bytes) (perhaps + 1) #define MIC_MAX_HOLD_TIME 400 // Time to hold the maximum mic level on the Status screen in milliseconds // If USE_GET_SIN is not zero, replace mic samples with a sin wave at a // frequency determined by the sidetone slider and an amplitude determined // by the Spot button level. LEVEL FAILS // If USE_GET_SIN is 1, pass these samples through the transmit filters. // If USE_GET_SIN is 2, transmit these samples directly. #define USE_GET_SIN 0 // If USE_2TONE is not zero, replace samples with a 2-tone test signal. #define USE_2TONE 0 #if USE_GET_SIN static void get_sin(complex double * cSamples, int count) { // replace mic samples with a sin wave int i; double freq; complex double phase1; // Phase increment static complex double vector1 = CLIP32 / 2; // Use the sidetone slider 0 to 1000 to set frequency //freq = (quisk_sidetoneCtrl - 500) / 1000.0 * MIC_OUT_RATE; freq = quisk_sidetoneCtrl * 5; freq = ((int)freq / 50) * 50; #if USE_GET_SIN == 2 phase1 = cexp(I * 2.0 * M_PI * freq / MIC_OUT_RATE); count *= MIC_OUT_RATE / quisk_sound_state.mic_sample_rate; #else phase1 = cexp(I * 2.0 * M_PI * freq / quisk_sound_state.mic_sample_rate); #endif for (i = 0; i < count; i++) { vector1 *= phase1; cSamples[i] = vector1; } #if DEBUG_IO || DEBUG if (debug_timer == 0) printf ("get_sin freq %.0lf\n", freq); #endif } #endif #if USE_2TONE static void get_2tone(complex double * cSamples, int count) { // replace mic samples int i; static complex double phase1=0, phase2; // Phase increment static complex double vector1; static complex double vector2; if (phase1 == 0) { // initialize phase1 = cexp((I * 2.0 * M_PI * IMD_TONE_1) / quisk_sound_state.mic_sample_rate); phase2 = cexp((I * 2.0 * M_PI * IMD_TONE_2) / quisk_sound_state.mic_sample_rate); vector1 = CLIP32 / 2.0; vector2 = CLIP32 / 2.0; } for (i = 0; i < count; i++) { vector1 *= phase1; vector2 *= phase2; cSamples[i] = (vector1 + vector2); } } #endif static double CcmPeak(double * dsamples, complex double * csamples, int count) { int i, j; complex double csample; double dtmp, dsample, newlevel, oldlevel; static double out_short, out_long; static struct Ccmpr { int buf_size; int index_read; double themax; double level; double * d_samp; complex double * c_samp; double * levl; } dat = {0}; if ( ! dat.buf_size) { // initialize; the sample rate is 8000 dat.buf_size = 8000 * 30 / 1000; // total delay in samples dat.index_read = 0; // index to output; and then write a new sample here dat.themax = 1.0; // maximum level in the buffer dat.level = 1.0; // current output level dat.d_samp = (double *) malloc(dat.buf_size * sizeof(double)); // buffer for double samples dat.c_samp = (complex double *) malloc(dat.buf_size * sizeof(complex double)); // buffer for complex samples dat.levl = (double *) malloc(dat.buf_size * sizeof(double)); // magnitude of the samples for (i = 0; i < dat.buf_size; i++) { dat.d_samp[i] = 0; dat.c_samp[i] = 0; dat.levl[i] = 1.0; } dtmp = 1.0 / 8000; // sample time out_short = 1.0 - exp(- dtmp / 0.010); // short time constant out_long = 1.0 - exp(- dtmp / 3.000); // long time constant return 1.0; } for (i = 0; i < count; i++) { if (dsamples) { dsample = dsamples[i]; dsamples[i] = dat.d_samp[dat.index_read] / dat.level; // FIFO output dat.d_samp[dat.index_read] = dsample; // write new sample at read index newlevel = fabs(dsample); } else { csample = csamples[i]; csamples[i] = dat.c_samp[dat.index_read] / dat.level; // FIFO output dat.c_samp[dat.index_read] = csample; // write new sample at read index newlevel = cabs(csample); } oldlevel = dat.levl[dat.index_read]; dat.levl[dat.index_read] = newlevel; if (newlevel < dat.themax && oldlevel < dat.themax) { // some other sample is the maximum // no change to dat.themax } else if (newlevel > dat.themax && newlevel > oldlevel) { // newlevel is the maximum dat.themax = newlevel; } else { // search for the maximum level dat.themax = 0; // Find the maximim level in the buffer for (j = 0; j < dat.buf_size; j++) { if (dat.levl[j] > dat.themax) dat.themax = dat.levl[j]; } } // Increase dat.level if the maximum level is greater than 1.0; // decrease it slowly back to 1.0 if it is lower. Output is modulated by dat.level. if (dat.themax > 1.0) // increase rapidly to the peak level dat.level = dat.level * (1.0 - out_short) + dat.themax * out_short; else // decrease slowly back to 1.0 dat.level = dat.level * (1.0 - out_long) + 1.0 * out_long; if (++dat.index_read >= dat.buf_size) dat.index_read = 0; } return dat.level; } static void init_alc(struct alc * pt, int size) { // Call first to set the buffer size. Then call to initialize the structure on each key down. int i; if (pt->buffer == NULL) { pt->buf_size = size; // do not change the size pt->buffer = (complex double *)malloc(size * sizeof(complex double)); for (i = 0; i < 20; i++) // initial gain by rx_mode switch(i) { case DGT_U: case DGT_L: case DGT_IQ: case FDV_U: case FDV_L: pt->gain_now[i] = 1.4; break; default: pt->gain_now[i] = 1.0; break; } } pt->index = 0; pt->block_index = 0; pt->counter = 0; pt->fault = 0; pt->max_magn = 0; pt->gain_max = 3.0; pt->gain_min = 0.1; pt->gain_change = 0; pt->next_change = 0; pt->final_gain = 0; for (i = 0; i < pt->buf_size; i++) pt->buffer[i] = 0; } static void process_alc(complex double * cSamples, int count, struct alc * pt, rx_mode_type rx_mode) { // Automatic Level Control (ALC) int i; double d, magn; complex double csamp; for (i = 0; i < count; i++) { csamp = cSamples[i]; // new sample to add to buffer cSamples[i] = pt->buffer[pt->index] * pt->gain_now[rx_mode]; // remove sample from buffer and apply gain #if DEBUG_LEVEL || DEBUG_IO || DEBUG magn = cabs(cSamples[i]); if (magn >= CLIP16) printf("ALC clip gain %9.6f level %9.6f\n", pt->gain_now[rx_mode], magn / CLIP16); #endif pt->buffer[pt->index] = csamp; // add new sample to buffer magn = cabs(csamp); // measure new sample if (magn * (pt->gain_now[rx_mode] + pt->gain_change * pt->buf_size) > (CLIP16 - 10)) { pt->gain_change = ((CLIP16 - 10) / magn - pt->gain_now[rx_mode]) / pt->buf_size; pt->final_gain = pt->gain_now[rx_mode] + pt->gain_change * pt->buf_size; if (pt->final_gain > pt->gain_max) { pt->final_gain = pt->gain_max; pt->gain_change = (pt->final_gain - pt->gain_now[rx_mode]) / pt->buf_size; } else if (pt->final_gain < pt->gain_min) { pt->final_gain = pt->gain_min; pt->gain_change = (pt->final_gain - pt->gain_now[rx_mode]) / pt->buf_size; } pt->block_index = pt->index; pt->counter = 0; pt->fault = 0; pt->next_change = 1E10; //printf("ALC DEC: gain %9.6f change %9.6f final gain %9.6f\n", pt->gain_now[rx_mode], pt->gain_change, pt->final_gain); } else if (pt->index == pt->block_index) { //d = cabs(cSamples[i]) / (CLIP16 - 10); //printf("ALC Fin: gain %9.6f out level %9.6f\n", pt->gain_now[rx_mode], d); #if 0 alc_start_gain = alc_gain; k = alc_block_index; for (j = 1; j < pt->buf_size; j++) { if (++k >= pt->buf_size) k = 0; magn = cabs(alc_buffer[k]); if (magn < 100) { alc_fault++; continue; } else { d = ((CLIP16 - 10) / magn - alc_start_gain) / j; if ( alc_next_change > d) alc_next_change = d; } } #endif d = 5.0; // number of seconds to double gain d = 1.0 / (48000.0 * d); if (pt->next_change > d) pt->next_change = d; if (pt->next_change != 1E10 && pt->fault < pt->buf_size - 10) { pt->gain_change = pt->next_change; } pt->final_gain = pt->gain_now[rx_mode] + pt->gain_change * pt->buf_size; if (pt->final_gain > pt->gain_max) { pt->final_gain = pt->gain_max; pt->gain_change = (pt->final_gain - pt->gain_now[rx_mode]) / pt->buf_size; } else if (pt->final_gain < pt->gain_min) { pt->final_gain = pt->gain_min; pt->gain_change = (pt->final_gain - pt->gain_now[rx_mode]) / pt->buf_size; } //printf("ALC New: gain %9.6f change %9.6f final gain %9.6f\n", pt->gain_now[rx_mode], pt->gain_change, pt->final_gain); pt->fault = 0; pt->counter = 0; pt->next_change = 1E10; } else { if (magn < 100) { pt->fault++; } else { d = ((CLIP16 - 10) / magn - pt->final_gain) / ++pt->counter; if ( pt->next_change > d) pt->next_change = d; } } pt->gain_now[rx_mode] += pt->gain_change; if (++pt->index >= pt->buf_size) pt->index = 0; } #if DEBUG_LEVEL || DEBUG_IO || DEBUG for (i = 0; i < count; i++) { magn = cabs(cSamples[i]); if (pt->max_magn < magn) pt->max_magn = magn; } if (debug_timer == 0) { printf("ALC Out: gain %9.6f max lvl%9.6f final gain %9.6f\n", pt->gain_now[rx_mode], pt->max_magn / CLIP16, pt->final_gain); pt->max_magn = 0; } #endif } static int tx_filter(complex double * filtered, int count) { // Input samples are creal(filtered), output is filtered. The input rate must be 8000 or 48000 sps. int i, is_ssb, mic_interp; int sample_rate = 8000; double dsample, dtmp, magn; complex double csample; static double inMax=0.3; static double x_1=0; static double aaa, bbb, ccc, Xmin, Xmax, Ymax; static int samples_size = 0; static double * dsamples = NULL; static complex double * csamples = NULL; static double time_long, time_short; static struct quisk_dFilter filtDecim, dfiltInterp; static struct quisk_dFilter filtAudio1, filtAudio2, dfiltAudio3; static struct quisk_cFilter cfiltAudio3, cfiltInterp; static struct quisk_dFilter filter1={NULL}, filter2; #if DEBUG_IO || DEBUG char * clip; static double dbOut = 0, Level0 = 0, Level1 = 0, Level2 = 0, Level3 = 0, Level4 = 0; #endif is_ssb = (rxMode == LSB || rxMode == USB); if (!filtered) { // initialization if (! filter1.dCoefs) { quisk_filt_dInit(&filter1, quiskMicFilt8Coefs, sizeof(quiskMicFilt8Coefs)/sizeof(double)); quisk_filt_dInit(&filter2, quiskMicFilt8Coefs, sizeof(quiskMicFilt8Coefs)/sizeof(double)); quisk_filt_dInit(&filtDecim, quiskLpFilt48Coefs, sizeof(quiskLpFilt48Coefs)/sizeof(double)); quisk_filt_dInit(&dfiltInterp, quiskLpFilt48Coefs, sizeof(quiskLpFilt48Coefs)/sizeof(double)); quisk_filt_cInit(&cfiltInterp, quiskLpFilt48Coefs, sizeof(quiskLpFilt48Coefs)/sizeof(double)); quisk_filt_dInit(&filtAudio1, quiskFiltTx8kAudioB, sizeof(quiskFiltTx8kAudioB)/sizeof(double)); quisk_filt_dInit(&filtAudio2, quiskFiltTx8kAudioB, sizeof(quiskFiltTx8kAudioB)/sizeof(double)); quisk_filt_dInit(&dfiltAudio3, quiskFiltTx8kAudioB, sizeof(quiskFiltTx8kAudioB)/sizeof(double)); quisk_filt_cInit(&cfiltAudio3, quiskFiltTx8kAudioB, sizeof(quiskFiltTx8kAudioB)/sizeof(double)); dtmp = 1.0 / sample_rate; // sample time time_long = 1.0 - exp(- dtmp / 3.000); time_short = 1.0 - exp(- dtmp / 0.005); Ymax = pow(10.0, - 1 / 20.0); // maximum y Xmax = pow(10.0, 3 / 20.0); // x where slope is zero; for x > Xmax, y == Ymax Xmin = Ymax - fabs(Ymax - Xmax); // x where slope is 1 and y = x; start of compression aaa = 1.0 / (2.0 * (Xmin - Xmax)); // quadratic bbb = -2.0 * aaa * Xmax; ccc = Ymax - aaa * Xmax * Xmax - bbb * Xmax; #if DEBUG_IO || DEBUG printf("Compress to %.2lf dB from %.2lf to %.2lf dB\n", 20 * log10(Ymax), 20 * log10(Xmin), 20 * log10(Xmax)); #endif } if (is_ssb) { quisk_filt_tune(&filter1, 1650.0 / sample_rate, rxMode != LSB); quisk_filt_tune(&filter2, 1650.0 / sample_rate, rxMode != LSB); } return 0; } mic_interp = MIC_OUT_RATE / sample_rate; // check size of dsamples[] and csamples[] buffer if (count > samples_size) { samples_size = count * 2 * mic_interp; if (dsamples) free(dsamples); if (csamples) free(csamples); dsamples = (double *)malloc(samples_size * sizeof(double)); csamples = (complex double *)malloc(samples_size * sizeof(complex double)); } // copy to dsamples[], normalize to +/- 1.0 for (i = 0; i < count; i++) dsamples[i] = creal(filtered[i]) / CLIP16; // Decimate to 8000 Hz if (quisk_sound_state.mic_sample_rate != sample_rate) count = quisk_dDecimate(dsamples, count, &filtDecim, quisk_sound_state.mic_sample_rate / sample_rate); // restrict bandwidth 300 to 2700 Hz count = quisk_dFilter(dsamples, count, &filtAudio1); #if DEBUG_IO || DEBUG // Measure peak input audio level for (i = 0; i < count; i++) { magn = fabs(dsamples[i]); if (magn > Level0) Level0 = magn; } #endif // high pass filter for preemphasis: See Radcom, January 2010, page 76. // quisk_mic_preemphasis == 1 was measured as 6 dB / octave. // gain at 800 Hz was measured as 0.104672. for (i = 0; i < count; i++) { dtmp = dsamples[i]; dsamples[i] = dtmp - quisk_mic_preemphasis * x_1; x_1 = dtmp; // delayed sample dsamples[i] *= 2; // compensate for loss #if DEBUG_IO || DEBUG magn = fabs(dsamples[i]); if (magn > Level1) Level1 = magn; #endif } if (is_ssb) { // SSB // FIR bandpass filter; separate into I and Q for (i = 0; i < count; i++) { csample = quisk_dC_out(dsamples[i], &filter1) * 2.0; // filter loss 0.5 // Measure average peak input audio level and normalize magn = cabs(csample); if (magn > inMax) inMax = inMax * (1 - time_short) + time_short * magn; else if(magn > mic_agc_level) inMax = inMax * (1 - time_long) + time_long * magn; else inMax = inMax * (1 - time_long) + time_long * mic_agc_level; csample /= inMax; magn /= inMax; #if DEBUG_IO || DEBUG if (magn > Level2) Level2 = magn; #endif // Audio compression. csample *= quisk_mic_clip; magn *= quisk_mic_clip; if (magn > 1.0) csample = csample / magn; dsamples[i] = creal(csample); } } else { // AM and FM // Measure average peak input audio level and normalize for (i = 0; i < count; i++) { dsample = dsamples[i]; magn = fabs(dsample); if (magn > inMax) inMax = inMax * (1 - time_short) + time_short * magn; else if(magn > mic_agc_level) inMax = inMax * (1 - time_long) + time_long * magn; else inMax = inMax * (1 - time_long) + time_long * mic_agc_level; dsample /= inMax; magn /= inMax; #if DEBUG_IO || DEBUG if (magn > Level2) Level2 = magn; #endif // Audio compression. dsample *= quisk_mic_clip; magn *= quisk_mic_clip; if (magn < Xmin) dsamples[i] = dsample; else if (magn > Xmax) dsamples[i] = copysign(Ymax, dsample); else dsamples[i] = copysign(aaa * magn * magn + bbb * magn + ccc, dsample); } } // remove clipping distortion; restrict bandwidth 300 to 2700 Hz count = quisk_dFilter(dsamples, count, &filtAudio2); if (is_ssb) { // SSB // FIR bandpass filter; separate into I and Q for (i = 0; i < count; i++) { csamples[i] = quisk_dC_out(dsamples[i], &filter2) * 2.0; // filter loss 0.5 #if DEBUG_IO || DEBUG magn = cabs(csamples[i]); if (magn > Level3) Level3 = magn; #endif } // round off peaks CcmPeak(NULL, csamples, count); #if DEBUG_IO || DEBUG for (i = 0; i < count; i++) { magn = cabs(csamples[i]); if (magn > Level4) Level4 = magn; } #endif // remove clipping distortion count = quisk_cDecimate(csamples, count, &cfiltAudio3, 1); // Interpolate up to 48000 if (mic_interp != 1) count = quisk_cInterpolate(csamples, count, &cfiltInterp, mic_interp); // convert back to 16 bits for (i = 0; i < count; i++) { filtered[i] = csamples[i] * CLIP16; #if DEBUG_IO || DEBUG magn = cabs(csamples[i]); if (magn > dbOut) dbOut = magn; #endif } } else { // AM and FM #if DEBUG_IO || DEBUG for (i = 0; i < count; i++) { magn = fabs(dsamples[i]); if (magn > Level3) Level3 = magn; } #endif // round off peaks CcmPeak(dsamples, NULL, count); #if DEBUG_IO || DEBUG for (i = 0; i < count; i++) { magn = fabs(dsamples[i]); if (magn > Level4) Level4 = magn; } #endif // remove clipping distortion count = quisk_dFilter(dsamples, count, &dfiltAudio3); // Interpolate up to 48000 if (mic_interp != 1) count = quisk_dInterpolate(dsamples, count, &dfiltInterp, mic_interp); // convert back to 16 bits for (i = 0; i < count; i++) { filtered[i] = dsamples[i] * CLIP16; #if DEBUG_IO || DEBUG magn = fabs(dsamples[i]); if (magn > dbOut) dbOut = magn; #endif } } #if DEBUG_IO || DEBUG if (debug_timer == 0) { if (dbOut > 1.0) clip = "Clip"; else clip = ""; dbOut = 20 * log10(dbOut); printf ("pre %3.1lf dB clip %2.0lf InMax %6.2lf Level0 %6.2lf Level1 %6.2lf Level2 %6.2lf Level3 %6.2lf Level4 %6.2lf dbOut %6.2lf %s\n", quisk_mic_preemphasis, 20 * log10(quisk_mic_clip), 20 * log10(inMax), 20 * log10(Level0), 20 * log10(Level1), 20 * log10(Level2), 20 * log10(Level3), 20 * log10(Level4), dbOut, clip); Level0 = Level1 = Level2 = Level3 = Level4 = dbOut = 0; } //QuiskPrintTime(" tx_filter", 2); #endif return count; } static int tx_filter_digital(complex double * filtered, int count) { // Input samples are creal(filtered), output is filtered. // This filter has minimal processing and is used for digital modes. int i; static int do_init = 1; static struct quisk_dFilter filter1; if (do_init) { // initialization do_init = 0; quisk_filt_dInit(&filter1, quiskDgtFilt48Coefs, sizeof(quiskDgtFilt48Coefs)/sizeof(double)); // pass 1350, stop 1650 } if ( ! filtered) { // Change to rxMode quisk_filt_tune(&filter1, 1650.0 / 48000, rxMode != DGT_L && rxMode != LSB && rxMode != FDV_L); return 0; } for (i = 0; i < count; i++) // FIR bandpass filter; separate into I and Q filtered[i] = quisk_dC_out(creal(filtered[i]), &filter1) * 2.00; // tuned filter loss 0.5 //quisk_calc_audio_graph(CLIP16, filtered, NULL, count, 0); return count; } PyObject * quisk_get_tx_filter(PyObject * self, PyObject * args) { // return the TX filter response to display on the graph // This is for debugging. Change quisk.py to call QS.get_tx_filter() instead // of QS.get_filter(). int i, j, k; int freq, time; PyObject * tuple2; complex double cx; double scale; double * average, * fft_window, * bufI, * bufQ; fftw_complex * samples, * pt; // complex data for fft fftw_plan plan; // fft plan double phase, delta; int nTaps = 325; if (!PyArg_ParseTuple (args, "")) return NULL; // Create space for the fft of size data_width pt = samples = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * data_width); plan = fftw_plan_dft_1d(data_width, pt, pt, FFTW_FORWARD, FFTW_MEASURE); average = (double *) malloc(sizeof(double) * (data_width + nTaps)); fft_window = (double *) malloc(sizeof(double) * data_width); bufI = (double *) malloc(sizeof(double) * nTaps); bufQ = (double *) malloc(sizeof(double) * nTaps); for (i = 0, j = -data_width / 2; i < data_width; i++, j++) // Hanning fft_window[i] = 0.5 + 0.5 * cos(2. * M_PI * j / data_width); for (i = 0; i < data_width + nTaps; i++) average[i] = 0.5; // Value for freq == 0 for (freq = 1; freq < data_width / 2.0 - 10.0; freq++) { //freq = data_width * 0.2 / 48.0; delta = 2 * M_PI / data_width * freq; phase = 0; // generate some initial samples to fill the filter pipeline for (time = 0; time < data_width + nTaps; time++) { average[time] += cos(phase); // current sample phase += delta; if (phase > 2 * M_PI) phase -= 2 * M_PI; } } // now filter the signal using the transmit filter tx_filter(NULL, 0); // initialize scale = 1.0; for (i = 0; i < data_width; i++) if (fabs(average[i + nTaps]) > scale) scale = fabs(average[i + nTaps]); scale = CLIP16 / scale; // limit to CLIP16 for (i = 0; i < nTaps; i++) samples[i] = average[i] * scale; tx_filter(samples, nTaps); // process initial samples for (i = 0; i < data_width; i++) samples[i] = average[i + nTaps] * scale; tx_filter(samples, data_width); // process the samples for (i = 0; i < data_width; i++) // multiply by window samples[i] *= fft_window[i]; fftw_execute(plan); // Calculate FFT // Normalize and convert to log10 scale = 0.3 / data_width / scale; for (k = 0; k < data_width; k++) { cx = samples[k]; average[k] = cabs(cx) * scale; if (average[k] <= 1e-7) // limit to -140 dB average[k] = -7; else average[k] = log10(average[k]); } // Return the graph data tuple2 = PyTuple_New(data_width); i = 0; // Negative frequencies: for (k = data_width / 2; k < data_width; k++, i++) PyTuple_SetItem(tuple2, i, PyFloat_FromDouble(20.0 * average[k])); // Positive frequencies: for (k = 0; k < data_width / 2; k++, i++) PyTuple_SetItem(tuple2, i, PyFloat_FromDouble(20.0 * average[k])); free(bufQ); free(bufI); free(average); free(fft_window); fftw_destroy_plan(plan); fftw_free(samples); return tuple2; } // Send samples using the Metis-Hermes protocol. A frame is 8 bytes: L/R audio and I/Q mic samples. // All samples are 2 bytes. The 1032 byte UDP packet contains 63*2 radio sound samples, and 63*2 mic I/Q samples. // Samples are sent synchronously with the input samples. static void quisk_hermes_tx_reset(void) { // Reset the buffer to half full of zero samples int i; hermes_num_samples = HERMES_TX_BUF_SAMPLES / 2; hermes_read_index = 0; hermes_write_index = hermes_num_samples * 2; for (i = 0; i < HERMES_TX_BUF_SHORTS; i++) // Put zero mic samples into the buffer hermes_buf[i] = 0; //printf("quisk_hermes_tx_reset: hermes_num_samples %d\n", hermes_num_samples); } static void quisk_hermes_tx_add(complex double * cSamples, int tx_count, int key_down) { // Add samples to the Tx buffer. int i; static int hermes_buf_has_samples=1; if (key_down) { // add non-zero samples to buffer hermes_buf_has_samples = 1; } else if (hermes_buf_has_samples) { // key is not down; reset buffer to zero quisk_hermes_tx_reset(); hermes_buf_has_samples = 0; return; } else { // key is not down but buffer is zeroed; just reset pointers hermes_num_samples = HERMES_TX_BUF_SAMPLES / 2; hermes_read_index = 0; hermes_write_index = hermes_num_samples * 2; return; } //printf("hermes_tx_add start: hermes_num_samples %d, tx_count %d\n", hermes_num_samples, tx_count); if (hermes_num_samples + tx_count >= HERMES_TX_BUF_SAMPLES) { // no more space; throw away half the samples quisk_udp_mic_error("Tx hermes buffer overflow"); //printf("Tx hermes buffer overflow: hermes_num_samples %d tx_count %d\n", hermes_num_samples, tx_count); i = hermes_num_samples - HERMES_TX_BUF_SAMPLES / 2; // number of samples to remove hermes_num_samples -= i; hermes_read_index += i * 2; if (hermes_read_index >= HERMES_TX_BUF_SHORTS) hermes_read_index -= HERMES_TX_BUF_SHORTS; } // Send the transmit samples to the predistortion logic. Modify the samples in place. PreDistort(NULL, NULL, 0, cSamples, tx_count); hermes_num_samples += tx_count; for (i = 0; i < tx_count; i++) { // Put transmit mic samples into the buffer hermes_buf[hermes_write_index++] = (short)cimag(cSamples[i]); hermes_buf[hermes_write_index++] = (short)creal(cSamples[i]); if (hermes_write_index >= HERMES_TX_BUF_SHORTS) hermes_write_index = 0; } //printf ("Buffer usage %.1f %%, hermes_num_samples %d, tx_count %d\n", 100.0 * hermes_num_samples / HERMES_TX_BUF_SAMPLES, hermes_num_samples, tx_count); //printf("hermes_tx_add end: hermes_num_samples %d\n", hermes_num_samples); } void quisk_hermes_tx_send(int tx_socket, int * tx_records) { // Send one UDP block of mic samples using the Metis-Hermes protocol. Timing is from blocks received, rate is 48k. // If the key is up we send samples anyway, but the samples are zero. int i, j, offset, sent, ratio; int PsInUse; short s; unsigned char sendbuf[1032]; unsigned char * pt_buf; static unsigned int seq = 0; static unsigned char C0_index = 0; complex double cw_samples[63 * 2]; static double writequeue_time0 = 0; static unsigned char rx1_freq[4]; // frequency of the first receiver //printf("hermes_tx_send start 1: hermes_num_samples %d\n", hermes_num_samples); if (tx_records == NULL) { seq = 0; C0_index = 0; quisk_hermes_tx_reset(); return; } if (quisk_play_state != last_play_state) { last_play_state = quisk_play_state; serial_key_samples(NULL, 0); } PsInUse = PsEnable || PsCal; //printf("hermes_tx_send start 2: hermes_num_samples %d\n", hermes_num_samples); ratio = quisk_sound_state.sample_rate / 48000; // send rate is 48 ksps //printf ("quisk_hermes_tx_send ratio %d count %d\n", ratio, *tx_records); if (*tx_records / ratio < 63 * 2) // tx_records is the number of samples received for each receiver return; // Send 63*2 Tx samples with control bytes if (quisk_play_state == SOFTWARE_CWKEY) { // Send CW samples serial_key_samples(cw_samples, 63 * 2); j = 0; for (i = 0; i < 63 * 2; i++) { hermes_buf[j++] = (short)cimag(cw_samples[i]); hermes_buf[j++] = (short)creal(cw_samples[i]); } hermes_read_index = 0; hermes_num_samples = 63 * 2; } //printf ("Buffer usage %.1f %%\n", 100.0 * hermes_num_samples / HERMES_TX_BUF_SAMPLES); //printf ("Tx quisk_hermes_tx_send ratio %d, count %d, samples %d\n", ratio, *tx_records, hermes_num_samples); *tx_records -= 63 * 2 * ratio; if (hermes_num_samples < 63 * 2) { // Not enough samples to send //printf("Tx hermes buffer underflow: hermes_num_samples %d\n", hermes_num_samples); quisk_udp_mic_error("Tx hermes buffer underflow"); quisk_hermes_tx_reset(); } hermes_num_samples -= 63 * 2; sendbuf[0] = 0xEF; sendbuf[1] = 0xFE; sendbuf[2] = 0x01; sendbuf[3] = 0x02; sendbuf[4] = seq >> 24 & 0xFF; sendbuf[5] = seq >> 16 & 0xFF; sendbuf[6] = seq >> 8 & 0xFF; sendbuf[7] = seq & 0xFF; seq++; sendbuf[8] = 0x7F; sendbuf[9] = 0x7F; sendbuf[10] = 0x7F; offset = C0_index * 4; // offset into quisk_pc_to_hermes is C0[7:1] * 4 sendbuf[11] = C0_index << 1 | hermes_mox_bit; // C0 sendbuf[12] = quisk_pc_to_hermes[offset++]; // C1 sendbuf[13] = quisk_pc_to_hermes[offset++]; // C2 sendbuf[14] = quisk_pc_to_hermes[offset++]; // C3 sendbuf[15] = quisk_pc_to_hermes[offset++]; // C4 switch (C0_index) { case 0: // Do not change receiver count without stopping Hermes and restarting sendbuf[15] = quisk_multirx_count << 3 | 0x04; // Send the old count, not the changed count if (hermes_mox_bit) // send filter selection on J16 sendbuf[13] = hermes_filter_tx << 1; else sendbuf[13] = hermes_filter_rx << 1; break; case 2: // Rx1 frequency rx1_freq[0] = sendbuf[12]; rx1_freq[1] = sendbuf[13]; rx1_freq[2] = sendbuf[14]; rx1_freq[3] = sendbuf[15]; break; case 3: // Rx2 frequency if (PsInUse && hermes_mox_bit) { sendbuf[12] = rx1_freq[0]; sendbuf[13] = rx1_freq[1]; sendbuf[14] = rx1_freq[2]; sendbuf[15] = rx1_freq[3]; } break; case 10: if (PsInUse && hermes_mox_bit) sendbuf[13] |= 0x40; // assert pure signal break; } if (++C0_index > 16) C0_index = 0; pt_buf = sendbuf + 16; for (i = 0; i < 63; i++) { // add 63 samples *pt_buf++ = 0x00; // Left/Right audio sample *pt_buf++ = 0x00; *pt_buf++ = 0x00; *pt_buf++ = 0x00; s = hermes_buf[hermes_read_index++]; *pt_buf++ = (s >> 8) & 0xFF; // Two bytes of I *pt_buf++ = s & 0xFF; s = hermes_buf[hermes_read_index++]; *pt_buf++ = (s >> 8) & 0xFF; // Two bytes of Q *pt_buf++ = s & 0xFF; if (hermes_read_index >= HERMES_TX_BUF_SHORTS) hermes_read_index = 0; } sendbuf[520] = 0x7F; sendbuf[521] = 0x7F; sendbuf[522] = 0x7F; // Changes for HermesLite v2 thanks to Steve, KF7O // Add a delay between ACK requests so two successive requests are spaced. if (writequeue_time0 == 0) writequeue_time0 = QuiskTimeSec() + 0.050; // initial delay if (quisk_hermeslite_writepointer == 1 && QuiskTimeSec() - writequeue_time0 > 0.020) { writequeue_time0 = QuiskTimeSec(); // Only send periodic hermeslite writes in second part of frame sendbuf[523] = quisk_hermeslite_writequeue[0] << 1 | hermes_mox_bit; sendbuf[524] = quisk_hermeslite_writequeue[1]; sendbuf[525] = quisk_hermeslite_writequeue[2]; sendbuf[526] = quisk_hermeslite_writequeue[3]; sendbuf[527] = quisk_hermeslite_writequeue[4]; if ((sendbuf[523] & 0x80) == 0) { // No acknowledge requested so fire and forget quisk_hermeslite_writepointer = 0; } else { quisk_hermeslite_writepointer = 2; } } else { offset = C0_index * 4; // offset into quisk_pc_to_hermes is C0[7:1] * 4 sendbuf[523] = C0_index << 1 | hermes_mox_bit; // C0 sendbuf[524] = quisk_pc_to_hermes[offset++]; // C1 sendbuf[525] = quisk_pc_to_hermes[offset++]; // C2 sendbuf[526] = quisk_pc_to_hermes[offset++]; // C3 sendbuf[527] = quisk_pc_to_hermes[offset++]; // C4 switch (C0_index) { case 0: // Do not change receiver count without stopping Hermes and restarting sendbuf[527] = quisk_multirx_count << 3 | 0x04; // Send the old count, not the changed count if (hermes_mox_bit) // send filter selection on J16 sendbuf[525] = hermes_filter_tx << 1; else sendbuf[525] = hermes_filter_rx << 1; break; case 2: // Rx1 frequency rx1_freq[0] = sendbuf[524]; rx1_freq[1] = sendbuf[525]; rx1_freq[2] = sendbuf[526]; rx1_freq[3] = sendbuf[527]; break; case 3: // Rx2 frequency if (PsInUse && hermes_mox_bit) { sendbuf[524] = rx1_freq[0]; sendbuf[525] = rx1_freq[1]; sendbuf[526] = rx1_freq[2]; sendbuf[527] = rx1_freq[3]; } break; case 10: if (PsInUse && hermes_mox_bit) sendbuf[525] |= 0x40; // assert pure signal break; } if (++C0_index > 16) C0_index = 0; } pt_buf = sendbuf + 528; for (i = 0; i < 63; i++) { // add 63 samples *pt_buf++ = 0x00; // Left/Right audio sample *pt_buf++ = 0x00; *pt_buf++ = 0x00; *pt_buf++ = 0x00; s = hermes_buf[hermes_read_index++]; *pt_buf++ = (s >> 8) & 0xFF; // Two bytes of I *pt_buf++ = s & 0xFF; s = hermes_buf[hermes_read_index++]; *pt_buf++ = (s >> 8) & 0xFF; // Two bytes of Q *pt_buf++ = s & 0xFF; if (hermes_read_index >= HERMES_TX_BUF_SHORTS) hermes_read_index = 0; } sent = send(tx_socket, (char *)sendbuf, 1032, 0); if (sent != 1032) quisk_udp_mic_error("Tx UDP socket error in Hermes"); if (quisk_play_state == SOFTWARE_CWKEY) quisk_hermes_tx_reset(); //printf("hermes_tx_send end: hermes_num_samples %d\n", hermes_num_samples); } // udp_iq has an initial zero followed by the I/Q samples. // The initial zero is sent iff align4 == 1. static void transmit_udp(complex double * cSamples, int count) { // Send count samples using the HiQSDR protocol. Each sample is sent as two shorts (4 bytes) of I/Q data. // Transmission is delayed until a whole block of data is available. int i, sent; static short udp_iq[TX_BLOCK_SHORTS + 1] = {0}; static int udp_size = 1; if (mic_socket == INVALID_SOCKET) return; if ( ! cSamples) { // initialization udp_size = 1; udp_iq[0] = 0; // should not be necessary return; } if (doTxCorrect) { for (i = 0; i < count; i++) cSamples[i] = cSamples[i] * TxCorrectLevel + TxCorrectDc; } for (i = 0; i < count; i++) { // transmit samples udp_iq[udp_size++] = (short)creal(cSamples[i]); udp_iq[udp_size++] = (short)cimag(cSamples[i]); if (udp_size >= TX_BLOCK_SHORTS) { // check count if (align4) sent = send(mic_socket, (char *)udp_iq, udp_size * 2, 0); else sent = send(mic_socket, (char *)udp_iq + 1, --udp_size * 2, 0); if (sent != udp_size * 2) QuiskPrintf("Send socket returned %d\n", sent); udp_size = 1; } } } static void transmit_mic_carrier(complex double * cSamples, int count, double level) { // send a CW carrier instead of mic samples #if 0 // transmit a saw tooth static int up = 1; static int ramp = 0; int i; int max = 5000; static int pause = 5000 / 10; int half_max = (int)(0.1 / 0.2339 * max); for (i = 0; i < count; i++) { cSamples[i] = (double)ramp / max * CLIP16; if (ramp == half_max && pause) { pause--; } else if (up) { if (ramp == max) { up = 0; pause = max / 10; } else { ramp += 1; } } else { if (ramp == 0) { up = 1; pause = max / 10; } else { ramp -= 1; } } //printf("Ramp %d\n", ramp); } #endif #if 1 // transmit a carrier equal to the number of samples int i; for (i = 0; i < count; i++) cSamples[i] = level * CLIP16; #endif #if 0 // replace mic samples with a sin wave int i; double freq; complex double phase1; // Phase increment static complex double vector1 = CLIP16 / 2; // Use the sidetone slider 0 to 1000 to set frequency freq = quisk_sidetoneCtrl * 5; freq = ((int)freq / 50) * 50; phase1 = cexp(I * 2.0 * M_PI * freq / 48000.0); //phase1 = conj(phase1); for (i = 0; i < count; i++) { vector1 *= phase1; cSamples[i] = vector1 * level; } #endif } static void transmit_mic_imd(complex double * cSamples, int count, double level) { // send a 2-tone test signal instead of mic samples int i; complex double v; static complex double phase1=0, phase2; // Phase increment static complex double vector1; static complex double vector2; if (phase1 == 0) { // initialize phase1 = cexp((I * 2.0 * M_PI * IMD_TONE_1) / MIC_OUT_RATE); phase2 = cexp((I * 2.0 * M_PI * IMD_TONE_2) / MIC_OUT_RATE); vector1 = CLIP16 / 2.0; vector2 = CLIP16 / 2.0; } for (i = 0; i < count; i++) { // transmit a carrier equal to the number of samples vector1 *= phase1; vector2 *= phase2; v = level * (vector1 + vector2); cSamples[i] = v; } } int quisk_process_microphone(int mic_sample_rate, complex double * cSamples, int count) { int i, sample, maximum, interp, mic_interp, key_down; double d, ctcss_delta, audio_scale, ctcss_scale; static int key_was_down=0, key_down_counter=0; static struct quisk_cFilter filtInterp={NULL}; static struct quisk_cFilter filtInterp2={NULL}; static struct quisk_cFilter filt240D4 = {NULL}; static struct quisk_cFilter filt300D6 = {NULL}; static struct quisk_cHB45Filter HalfBand = {NULL, 0, 0}; static double ctcss_angle=0; static struct alc tx_alc = {NULL}; #if 0 // Measure soundcard actual sample rate static time_t seconds = 0; static int total = 0; struct timeval tb; static double dtime; gettimeofday(&tb); total += count; if (seconds == 0) { seconds = tb.tv_sec; dtime = tb.tv_sec + 0.000001 * tb.tv_usec; } else if (tb.tv_sec - seconds > 4) { printf("Mic soundcard rate %.3f\n", total / (tb.tv_sec + .000001 * tb.tv_usec - dtime)); seconds = tb.tv_sec; printf("backlog %d, count %d\n", backlog, count); } #endif #if DEBUG_IO || DEBUG //QuiskPrintTime("", -1); #endif #if DEBUG_IO || DEBUG double magn; static double out_max=0; #endif #if DEBUG_LEVEL || DEBUG_IO || DEBUG debug_timer += count; if (debug_timer >= mic_sample_rate) // one second debug_timer = 0; #endif // Microphone sample are input at mic_sample_rate. But after processing, // the output rate is MIC_OUT_RATE. interp = MIC_OUT_RATE / mic_sample_rate; #if USE_GET_SIN get_sin(cSamples, count); // Replace mic samples with a sin wave #endif #if USE_2TONE get_2tone(cSamples, count); // Replace mic samples with a 2-tone test signal #endif // measure maximum microphone level maximum = 1; for (i = 0; i < count; i++) { cSamples[i] *= (double)CLIP16 / CLIP32; // convert 32-bit samples to 16 bits d = creal(cSamples[i]); sample = (int)fabs(d); if (sample > maximum) maximum = sample; } // VOX processing if (maximum > vox_level) { is_vox = mic_sample_rate / 1000 * timeVOX; // reset timer to maximum } else if(is_vox) { is_vox -= count; // decrement timer if (is_vox < 0) is_vox = 0; } // mic display level if (maximum > mic_level) mic_level = maximum; mic_timer -= count; // time out the max microphone level to display if (mic_timer <= 0) { mic_timer = mic_sample_rate / 1000 * MIC_MAX_HOLD_TIME; mic_max_display = mic_level; mic_level = 1; } if ( ! tx_alc.buffer) init_alc(&tx_alc, 960); // at 48000 sps, 960 is 20 msec // quiskTxHoldState is a state machine to implement a pause for a repeater frequency shift for FM key_down = quisk_is_key_down(); if (rxMode == FM || rxMode == DGT_FM) { switch (quiskTxHoldState) { case 0: // Never implement any hold break; case 1: // Start hold when key goes down if (key_down) quiskTxHoldState = 2; break; case 2: // Key down hold is in progress; wait until state changes to 3 break; case 3: // Hold is released; when key goes up, hold starts again if ( ! key_down) quiskTxHoldState = 4; break; case 4: // Key up hold is in progress; wait until state changes to 1 break; } } if (quiskTxHoldState == 2 || quiskTxHoldState == 4) { // don't transmit until the hold is cleared key_down = 0; for (i = 0; i < count; i++) cSamples[i] = 0; } if (key_down != key_was_down) { key_down_counter = quisk_start_ssb_delay * 48; // Key was pressed. Zero the first few samples to clear the buffers. init_alc(&tx_alc, 0); // init ALC for Tx and TMP_RECORD_MIC key_was_down = key_down; } if (quisk_play_state != last_play_state) { last_play_state = quisk_play_state; serial_key_samples(NULL, 0); } if (key_down || DEBUG_MIC) { // create transmit I/Q samples #if USE_GET_SIN == 2 transmit_udp(cSamples, count * interp); #else if (quiskSpotLevel >= 0) { // Spot is in use count *= interp; transmit_mic_carrier(cSamples, count, quiskSpotLevel / 1000.0); } else if (tci_tx_audio_client) { count = tx_filter_digital(cSamples, count); // filter samples, minimal processing process_alc(cSamples, count, &tx_alc, DGT_U); } else switch (rxMode) { case LSB: // LSB case USB: // USB if (quisk_record_state == TMP_PLAY_SPKR_MIC) count = tx_filter_digital(cSamples, count); // filter samples, minimal processing else count = tx_filter(cSamples, count); // filter samples process_alc(cSamples, count, &tx_alc, rxMode); break; case AM: // AM if (quisk_record_state != TMP_PLAY_SPKR_MIC) // no audio processing for recorded sound count = tx_filter(cSamples, count); for (i = 0; i < count; i++) // transmit (0.5 + ampl/2, 0) cSamples[i] = (creal(cSamples[i]) + CLIP16) * 0.5; process_alc(cSamples, count, &tx_alc, rxMode); break; case FM: // FM if (quisk_record_state != TMP_PLAY_SPKR_MIC) // no audio processing for recorded sound count = tx_filter(cSamples, count); if (quisk_ctcss_freq > 9) { ctcss_delta = 2.0 * M_PI / MIC_OUT_RATE * quisk_ctcss_freq; ctcss_scale = 450.0 * modulation_index / quisk_ctcss_freq; // for 15% of total deviation audio_scale = 0.85 * modulation_index / CLIP16; for (i = 0; i < count; i++) { ctcss_angle += ctcss_delta; if (ctcss_angle >= 2.0 * M_PI) ctcss_angle -= 2.0 * M_PI; cSamples[i] = CLIP16 * cexp(I * (audio_scale * creal(cSamples[i]) + ctcss_scale * sin (ctcss_angle))); } } else { audio_scale = modulation_index / CLIP16; for (i = 0; i < count; i++) // this is phase modulation == FM and 6 dB /octave preemphasis cSamples[i] = CLIP16 * cexp(I * audio_scale * creal(cSamples[i])); } process_alc(cSamples, count, &tx_alc, rxMode); break; case DGT_U: // external digital modes case DGT_L: case DGT_IQ: case DGT_FM: case FDV_U: // FreeDV case FDV_L: count = tx_filter_digital(cSamples, count); // filter samples, minimal processing process_alc(cSamples, count, &tx_alc, rxMode); break; case IMD: // transmit IMD 2-tone test count *= interp; transmit_mic_imd(cSamples, count, quiskImdLevel / 1000.0); break; default: break; } #if DEBUG_IO || DEBUG for (i = 0; i < count; i++) { magn = cabs(cSamples[i]); if (out_max < magn) out_max = magn; } if (debug_timer == 0) { printf("Max cSamples[] Output Level %9.6f\n", out_max / CLIP16); out_max = 0; } #endif if (key_down_counter > 0) { // zero the first few samples to clear the buffers. for (i = 0; i < count && key_down_counter > 0; i++, key_down_counter--) { if (key_down_counter < 480) cSamples[i] *= (1.0 - key_down_counter / 480.0); // slow increase at end else cSamples[i] = 0; // initial samples are zero } } #endif } if (reverse_tx_sideband) for (i = 0; i < count; i++) cSamples[i] = conj(cSamples[i]); if (quisk_use_rx_udp == 10) ; // hermes handles this itself else if (quisk_play_state == SOFTWARE_CWKEY) serial_key_samples(cSamples, count); if(quisk_pt_sample_write) { // Used for SoapySDR // Interpolate the mic samples to the Tx sample rate //printf("Tx sample rate %i\n", tx_sample_rate); switch (tx_sample_rate) { case 100000: count = quisk_cInterp2HB45(cSamples, count, &HalfBand); // Fall through case 50000: if (! filt240D4.dCoefs) quisk_filt_cInit(&filt240D4, quiskFilt240D4Coefs, sizeof(quiskFilt240D4Coefs)/sizeof(double)); if (! filt300D6.dCoefs) quisk_filt_cInit(&filt300D6, quiskFilt300D6Coefs, sizeof(quiskFilt300D6Coefs)/sizeof(double)); count = quisk_cInterpDecim(cSamples, count, &filt240D4, 5, 4); // 60 kSps count = quisk_cInterpDecim(cSamples, count, &filt300D6, 5, 6); // 50 kSps break; case 96000: if (! filtInterp2.dCoefs) quisk_filt_cInit(&filtInterp2, quiskFilt48dec24Coefs, sizeof(quiskFilt48dec24Coefs)/sizeof(double)); count = quisk_cInterpolate(cSamples, count, &filtInterp2, 2); break; case 192000: if (! filtInterp2.dCoefs) quisk_filt_cInit(&filtInterp2, quiskFilt48dec24Coefs, sizeof(quiskFilt48dec24Coefs)/sizeof(double)); count = quisk_cInterp2HB45(cSamples, count, &HalfBand); count = quisk_cInterpolate(cSamples, count, &filtInterp2, 2); break; default: // 48000 break; } (*quisk_pt_sample_write)(cSamples, count); } else if (quisk_use_rx_udp == 10) { // Send Hermes mic samples when key is up or down if ( ! quisk_rx_udp_started) ; else { if ((rxMode == CWL || rxMode == CWU) && quiskSpotLevel < 0) // CW and no Spot for (i = 0; i < count; i++) cSamples[i] = 0; quisk_hermes_tx_add(cSamples, count, key_down); } } else if (quisk_use_rx_udp && key_down) { // Send mic samples to UDP when key is down transmit_udp(cSamples, count); } if (quisk_record_state == TMP_RECORD_MIC) { switch (rxMode) { case LSB: // LSB case USB: // USB count = tx_filter(cSamples, count); // filter samples process_alc(cSamples, count, &tx_alc, rxMode); break; case AM: // AM count = tx_filter(cSamples, count); process_alc(cSamples, count, &tx_alc, rxMode); break; case FM: // FM count = tx_filter(cSamples, count); process_alc(cSamples, count, &tx_alc, rxMode); break; default: for (i = 0; i < count; i++) cSamples[i] = 0; break; } // Perhaps interpolate the mic samples back to the sound play rate mic_interp = quisk_sound_state.playback_rate / MIC_OUT_RATE; if (mic_interp > 1) { if (! filtInterp.dCoefs) quisk_filt_cInit(&filtInterp, quiskFilt12_19Coefs, sizeof(quiskFilt12_19Coefs)/sizeof(double)); count = quisk_cInterpolate(cSamples, count, &filtInterp, mic_interp); } quisk_tmp_record(cSamples, count, (double)CLIP32 / CLIP16); // convert 16 to 32 bits } #if DEBUG_IO || DEBUG //QuiskPrintTime(" process_mic", 1); #endif return count; } PyObject * quisk_get_tx_audio(PyObject * self, PyObject * args) { const char * name; if (!PyArg_ParseTuple (args, "s", &name)) return NULL; if (strcmp(name, "PsEnable") == 0) return PyInt_FromLong(PsEnable); if (strcmp(name, "PsCal") == 0) return PyInt_FromLong(PsCal); Py_INCREF (Py_None); return Py_None; } PyObject * quisk_set_tx_audio(PyObject * self, PyObject * args, PyObject * keywds) { /* Call with keyword arguments ONLY; change Tx audio parameters */ static char * kwlist[] = {"vox_level", "vox_time", "mic_clip", "mic_preemphasis", "tx_sample_rate", "reverse_tx_sideband", "PsEnable", "PsCal", NULL} ; int vlevel = -9999, clevel = -9999; if (!PyArg_ParseTupleAndKeywords (args, keywds, "|iiidiiii", kwlist, &vlevel, &timeVOX, &clevel, &quisk_mic_preemphasis, &tx_sample_rate, &reverse_tx_sideband, &PsEnable, &PsCal)) return NULL; if (vlevel != -9999) vox_level = (int)(pow(10.0, vlevel / 20.0) * CLIP16); // Convert dB to 16-bit sample if (clevel != -9999) quisk_mic_clip = pow(10.0, clevel / 20.0); // Convert dB to factor Py_INCREF (Py_None); return Py_None; } PyObject * quisk_set_udp_tx_correct(PyObject * self, PyObject * args) // Called from GUI thread { double DcI, DcQ, level; if (!PyArg_ParseTuple (args, "ddd", &DcI, &DcQ, &level)) return NULL; if (DcI == 0 && DcQ == 0 && level == 1.0){ doTxCorrect = 0; } else { doTxCorrect = 1; TxCorrectDc = (DcI + I * DcQ) * CLIP16; DcI = fabs(DcI); DcQ = fabs(DcQ); if (DcI > DcQ) TxCorrectLevel = 1.0 - DcI; else TxCorrectLevel = 1.0 - DcQ; TxCorrectLevel *= level; } Py_INCREF (Py_None); return Py_None; } PyObject * quisk_is_vox(PyObject * self, PyObject * args) { /* return the VOX state */ if (!PyArg_ParseTuple (args, "")) return NULL; return PyInt_FromLong(is_vox); } PyObject * quisk_set_hermes_filter(PyObject * self, PyObject * args) { if (!PyArg_ParseTuple (args, "ii", &hermes_filter_rx, &hermes_filter_tx)) return NULL; Py_INCREF (Py_None); return Py_None; } void quisk_close_mic(void) { if (mic_socket != INVALID_SOCKET) { close(mic_socket); mic_socket = INVALID_SOCKET; } #ifdef MS_WINDOWS if (mic_cleanup) WSACleanup(); #endif } void quisk_open_mic(void) { struct sockaddr_in Addr; int sndsize = 48000; #if DEBUG_IO || DEBUG int intbuf; #ifdef MS_WINDOWS int bufsize = sizeof(int); #else socklen_t bufsize = sizeof(int); #endif #endif #ifdef MS_WINDOWS WORD wVersionRequested; WSADATA wsaData; #endif modulation_index = QuiskGetConfigDouble("modulation_index", 1.6); mic_agc_level = QuiskGetConfigDouble("mic_agc_level", 0.1); if (quisk_sound_state.tx_audio_port == 0x553B) align4 = 0; // Using old port: data starts at byte 42. else align4 = 1; // Start data at byte 44; align to dword if (quisk_sound_state.mic_ip[0]) { #ifdef MS_WINDOWS wVersionRequested = MAKEWORD(2, 2); if (WSAStartup(wVersionRequested, &wsaData) != 0) return; // failure to start winsock mic_cleanup = 1; #endif mic_socket = socket(PF_INET, SOCK_DGRAM, 0); if (mic_socket != INVALID_SOCKET) { setsockopt(mic_socket, SOL_SOCKET, SO_SNDBUF, (char *)&sndsize, sizeof(sndsize)); Addr.sin_family = AF_INET; // This is the UDP port for TX microphone samples, and must agree with the microcontroller. Addr.sin_port = htons(quisk_sound_state.tx_audio_port); #ifdef MS_WINDOWS Addr.sin_addr.S_un.S_addr = inet_addr(quisk_sound_state.mic_ip); #else inet_aton(quisk_sound_state.mic_ip, &Addr.sin_addr); #endif if (connect(mic_socket, (const struct sockaddr *)&Addr, sizeof(Addr)) != 0) { close(mic_socket); mic_socket = INVALID_SOCKET; } else { #if DEBUG_IO || DEBUG if (getsockopt(mic_socket, SOL_SOCKET, SO_SNDBUF, (char *)&intbuf, &bufsize) == 0) printf("UDP mic socket send buffer size %d\n", intbuf); else printf ("Failure SO_SNDBUF\n"); #endif } } } } void quisk_set_tx_mode(void) // called when the mode rxMode is changed { tx_filter(NULL, 0); tx_filter_digital(NULL, 0); transmit_udp(NULL, 0); serial_key_samples(NULL, 0); } #define CW_RISE_MSEC 5 #define CW_DELAY_SAMPLES (START_CW_DELAY_MAX * 48) static void serial_key_samples(complex double * cSamples, int count) // called from the sound thread { // Internal CW from serial port key - fill cSamples with CW samples int i, key_down, delay; static double ampl = 0; static char delay_line[CW_DELAY_SAMPLES]; // play CW a fixed delay after the key static int delay_index; double delta = 1.0 / (CW_RISE_MSEC * 48.0); double themax = 1.0; if ( ! cSamples) { // initialize for (i = 0; i < CW_DELAY_SAMPLES; i++) delay_line[i] = 0; delay_index = 0; return; } delay = quisk_start_cw_delay * 48; for (i = 0; i < count; i++) { key_down = delay_line[delay_index]; delay_line[delay_index] = QUISK_CWKEY_DOWN; if (++delay_index >= delay) { delay_index = 0; } if (key_down) { // key is down if (ampl < themax) { ampl += delta; if (ampl > themax) ampl = themax; } } else { // key is up if (ampl > 0.0) { ampl -= delta; if (ampl < 0.0) ampl = 0.0; } } cSamples[i] = ampl * CLIP16; } } static void spline(double x[], double y[], int n, double yp1, double ypn, double y2[]) { int i,k; double p,qn,sig,un,*u; u = (double *)malloc(n * sizeof(double)); if (yp1 > 0.99e30) y2[1]=u[1]=0.0; else { y2[1] = -0.5; u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1); } for (i=2;i<=n-1;i++) { sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]); p=sig*y2[i-1]+2.0; y2[i]=(sig-1.0)/p; u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]); u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p; } if (ypn > 0.99e30) qn=un=0.0; else { qn=0.5; un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1])); } y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0); for (k=n-1;k>=1;k--) y2[k]=y2[k]*y2[k+1]+u[k]; free(u); } static void polint(double xa[], double ya[], double x, double *y) { int i,m,ns=1; double den,dif,dift,ho,hp,w, dy; double c[4], d[4]; dif=fabs(x-xa[1]); for (i=1;i<=3;i++) { if ( (dift=fabs(x-xa[i])) < dif) { ns=i; dif=dift; } c[i]=ya[i]; d[i]=ya[i]; } *y=ya[ns--]; for (m=1;m<3;m++) { for (i=1;i<=3-m;i++) { ho=xa[i]-x; hp=xa[i+m]-x; w=c[i+1]-d[i]; den = ho - hp; if ( (den=ho-hp) == 0.0) printf("Error in routine polint\n"); den=w/den; d[i]=hp*den; c[i]=ho*den; } *y += (dy=(2*ns < (3-m) ? c[ns+1] : d[ns--])); } } static void splint(double xa[], double ya[], double y2a[], int n, double x, double *y) { static int klo = 1; int khi,k; double h,b,a; khi = klo + 1; if (x < xa[klo] || x > xa[khi]) { klo=1; khi=n; while (khi-klo > 1) { k=(khi+klo) >> 1; if (xa[k] > x) khi=k; else klo=k; } } h=xa[khi]-xa[klo]; if (h == 0.0) printf("Bad xa input to routine splint\n"); a=(xa[khi]-x)/h; b=(x-xa[klo])/h; *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0; } #define NUM_BINS 21 #define BIN_SPACING (1.0 / (NUM_BINS - 1)) #define BIN_TOL (BIN_SPACING * 0.20) #define BLOCK_MSECS 2000 // Time of one data block. #define MILLISECS(x) (48 * (x)) // Convert milliseconds to samples at 48000 samples per second. #define DEBUG_PS 0 #define HOLD 0 // Make calibration run permanent #define RX1_GAIN 0.238 void PreDistort(complex double * amp_in_samples, complex double * amp_out_samples, int nSamples, complex double * tx_samples, int num_tx) { int i, bin, num_bin_data; double G0, phase, mag_in, mag_out, new_mag, aaa, ppp; double error, xp, xpg, xp0, xpg0; double quadX[4], quadY[4], quadZ[4]; complex double in_sample, out_sample; static int state = 0; static int Fault; static int isCalibrated = 0; static int ramp_counter, ramp_count0; static int input_samples = 0; static double tx_to_RX1_gain = 0; static double max_out_mag, starting_phase; #if DEBUG_PS static double clip_tx, clip_mag_in, clip_aaa; #endif #if DEBUG_PS > 2 static int old_state; static double time0; #endif static int spline_P_size; static double spline_P_X[NUM_BINS + 1]; static double spline_P_mag[NUM_BINS + 1]; static double spline_P_mag2[NUM_BINS + 1]; static double spline_P_phase[NUM_BINS + 1]; static double spline_P_phase2[NUM_BINS + 1]; struct BPD { int XP_count; double table_XP; double XP_mag; double mag_X; double mag_P; double phase_P; complex double XP_gain; double XP_gain_mag; double XP_gain_phase; }; static struct BPD Bins[NUM_BINS]; static enum { NothingHappens, SendCW, StopCW } tx_state = NothingHappens; if (quisk_multirx_count == 0) { state = 0; return; } if (PsCal) { if (state == 0) { isCalibrated = 0; state = 10; #if DEBUG_PS > 2 time0 = QuiskTimeSec(); #endif } } else if ( ! PsEnable) { state = 0; return; } else if ( ! hermes_mox_bit) { state = 0; return; } if (num_tx) { input_samples += num_tx; if (isCalibrated) { // perform predistortion of Tx samples; multiply input X by P for (i = 0; i < num_tx; i++) { mag_in = cabs(tx_samples[i]) / CLIP16; splint(spline_P_X, spline_P_mag, spline_P_mag2, spline_P_size, mag_in, &aaa); splint(spline_P_X, spline_P_phase, spline_P_phase2, spline_P_size, mag_in, &ppp); tx_samples[i] *= aaa * cexp(I * ppp); new_mag = mag_in * aaa; if (new_mag > 1.0) { tx_samples[i] /= new_mag; #if DEBUG_PS if (new_mag > 1.0001 && new_mag > clip_tx) { clip_tx = new_mag; clip_mag_in = mag_in; clip_aaa = aaa; } #endif } } } switch (tx_state) { case NothingHappens: default: break; case SendCW: // Start sending a CW carrier. for (i = 0; i < num_tx; i++) { if (ramp_counter > 0) // slow increase tx_samples[i] = CLIP16 * (1.0 - (double)ramp_counter-- / ramp_count0); else tx_samples[i] = CLIP16; } break; case StopCW: // Stop sending a CW carrier. for (i = 0; i < num_tx; i++) { if (ramp_counter > 0) // slow decrease tx_samples[i] = (double)CLIP16 * ramp_counter-- / ramp_count0; else tx_samples[i] = 0; } break; } } if (nSamples <= 0) return; #if DEBUG_PS > 2 if (state != old_state) { aaa = cabs(amp_in_samples[0] / CLIP32); printf("State %2d time %6.3lf sample %12.6f\n", old_state, QuiskTimeSec() - time0, aaa); old_state = state; time0 = QuiskTimeSec(); } #endif switch (state) { case 0: // Start a timer tx_state = NothingHappens; input_samples = 0; starting_phase = 0; if ( ! HOLD) state++; break; case 1: if (input_samples >= MILLISECS(200)) { #if DEBUG_PS clip_tx = 0.0; #endif for (bin = 0; bin < NUM_BINS; bin++) { // initialize Bins[bin].XP_count = 0; Bins[bin].XP_mag = 0; Bins[bin].XP_gain = 0; } input_samples = 0; state++; } break; case 2: case 15: if (input_samples >= MILLISECS(BLOCK_MSECS)) { num_bin_data = 0; Fault = 0; for (bin = 0; bin < NUM_BINS; bin++) { if (Bins[bin].XP_count) { num_bin_data++; Bins[bin].XP_mag /= Bins[bin].XP_count; Bins[bin].XP_gain /= Bins[bin].XP_count; if (starting_phase == 0 && bin >= NUM_BINS / 2) { starting_phase = carg(Bins[bin].XP_gain); #if DEBUG_PS printf("starting_phase %.3f\n", starting_phase); #endif } } } if (PsCal && num_bin_data < NUM_BINS - 4) { printf("PS calibration failed to fill all bins "); for (bin = 0; bin < NUM_BINS; bin++) if (Bins[bin].XP_count) printf("1"); else printf("0"); printf("\n"); Fault = 1; } state++; } else { for (i = 0; i < nSamples; i++) { in_sample = amp_in_samples[i] / CLIP32 / RX1_GAIN; mag_in = cabs(in_sample); if (tx_to_RX1_gain < mag_in) { tx_to_RX1_gain = mag_in; #if DEBUG_PS > 2 printf("tx RX1 gain %.6f, max_out_mag %.6f\n", tx_to_RX1_gain, max_out_mag); #endif } bin = (int)(mag_in / BIN_SPACING + 0.5); if ( ! bin) continue; error = fabs(mag_in - Bins[bin].table_XP); if (error < BIN_TOL) { out_sample = amp_out_samples[i] / CLIP32; Bins[bin].XP_gain += out_sample / in_sample; Bins[bin].XP_mag += mag_in; Bins[bin].XP_count++; } } if (state == 15) { if (tx_state != StopCW) { ramp_counter = ramp_count0 = MILLISECS(BLOCK_MSECS - 500); tx_state = StopCW; } } } break; case 3: case 16: #if DEBUG_PS if (clip_tx != 0.0) { printf ("clip output sample %.6f, mag_in %.6f, predistort %.6f\n", clip_tx, clip_mag_in, clip_aaa); clip_tx = 0.0; } #endif for (bin = 1; bin < NUM_BINS - 1; bin++) { if (Bins[bin].XP_count > 0) { // This bin has data. Bins[bin].XP_gain_mag = cabs(Bins[bin].XP_gain); phase = carg(Bins[bin].XP_gain) - starting_phase; if (phase < - M_PI) phase += 2 * M_PI; else if (phase > M_PI) phase -= 2 * M_PI; Bins[bin].XP_gain_phase = phase; } else { Bins[bin].XP_gain_mag = 1.0; Bins[bin].XP_gain_phase = 0.0; } } Bins[0].XP_count = 1; Bins[0].XP_mag = 0.0; i = 1; for (bin = 1; bin < 6; bin++) { if (Bins[bin].XP_count) { quadX[i] = Bins[bin].XP_mag; quadY[i] = Bins[bin].XP_gain_mag; quadZ[i] = Bins[bin].XP_gain_phase; if (i == 3) break; i++; } } if (i == 3) { polint(quadX, quadY, 0.0, &Bins[0].XP_gain_mag); polint(quadX, quadZ, 0.0, &Bins[0].XP_gain_phase); } else { Bins[0].XP_gain_mag = 1.0; Bins[0].XP_gain_phase = 0.0; Fault = 1; if (DEBUG_PS || PsCal) printf("Failure to find G0 at zero\n"); } Bins[NUM_BINS - 1].XP_count = 1; Bins[NUM_BINS - 1].XP_mag = 1.0; i = 3; for (bin = NUM_BINS - 2; bin >= NUM_BINS - 6; bin--) { if (Bins[bin].XP_count) { quadX[i] = Bins[bin].XP_mag; quadY[i] = Bins[bin].XP_gain_mag; quadZ[i] = Bins[bin].XP_gain_phase; if (--i == 0) break; } } if (i == 0) { polint(quadX, quadY, 1.0, &Bins[NUM_BINS - 1].XP_gain_mag); polint(quadX, quadZ, 1.0, &Bins[NUM_BINS - 1].XP_gain_phase); } else { Bins[NUM_BINS - 1].XP_gain_mag = 1.0; Bins[NUM_BINS - 1].XP_gain_phase = 0.0; Fault = 1; if (DEBUG_PS || PsCal) printf("Failure to find G0 at one\n"); } state++; break; case 4: case 17: G0 = Bins[NUM_BINS - 1].XP_gain_mag; num_bin_data = 0; xp0 = -1.0; xpg0 = -1.0; for (bin = 0; bin < NUM_BINS; bin++) { if (bin == NUM_BINS - 1) { num_bin_data++; Bins[NUM_BINS - 1].mag_P = 1.0; Bins[NUM_BINS - 1].phase_P = - Bins[NUM_BINS - 1].XP_gain_phase; Bins[NUM_BINS - 1].mag_X = 1.0; } else if (Bins[bin].XP_count) { num_bin_data++; Bins[bin].mag_P = G0 / Bins[bin].XP_gain_mag; Bins[bin].phase_P = - Bins[bin].XP_gain_phase; Bins[bin].mag_X = Bins[bin].XP_mag / Bins[bin].mag_P; xp = Bins[bin].mag_X * Bins[bin].mag_P; if (xp < xp0) { if (DEBUG_PS || PsCal) printf("xp is not monotonic at bin %d\n", bin); Fault = 1; } xpg = xp * Bins[bin].XP_gain_mag; if (xpg < xpg0) { if (DEBUG_PS || PsCal) printf("xpg is not monotonic at bin %d\n", bin); Fault = 1; } xp0 = xp; xpg0 = xpg; } else { Bins[bin].mag_P = 0.0; Bins[bin].phase_P = 0.0; Bins[bin].mag_X = 0.0; } } if (num_bin_data < NUM_BINS * 7 / 10) { if (DEBUG_PS || state == 17) printf("Not enough data for corrections %d\n", num_bin_data); Fault = 1; } if ( ! Fault) { #if DEBUG_PS printf("Create a new predistortion\n"); #endif spline_P_size = 0; for (bin = 0; bin < NUM_BINS; bin++) { if (Bins[bin].XP_count) { spline_P_size++; spline_P_X[spline_P_size] = Bins[bin].mag_X; spline_P_mag[spline_P_size] = Bins[bin].mag_P; spline_P_phase[spline_P_size] = Bins[bin].phase_P; } } spline(spline_P_X, spline_P_mag, spline_P_size, 2E30, 2E30, spline_P_mag2); spline(spline_P_X, spline_P_phase, spline_P_size, 2E30, 2E30, spline_P_phase2); } state++; break; case 5: case 18: #if DEBUG_PS if (Fault) printf("Data for predistortion was rejected\n"); printf("\n"); #endif input_samples = 0; state++; break; case 6: if ( ! HOLD) state = 1; input_samples = 0; break; case 10: isCalibrated = 0; input_samples = 0; quisk_set_key_down(1); max_out_mag = 0; tx_to_RX1_gain = 0; starting_phase = 0; ramp_counter = ramp_count0 = MILLISECS(10); #if DEBUG_PS clip_tx = 0.0; #endif for (bin = 0; bin < NUM_BINS; bin++) { Bins[bin].table_XP = BIN_SPACING * bin; Bins[bin].XP_mag = 0; Bins[bin].XP_count = 0; Bins[bin].XP_gain = 0; } tx_state = SendCW; state++; break; case 11: if (input_samples >= MILLISECS(400)) { input_samples = 0; state++; } if (input_samples >= MILLISECS(50) && cabs(amp_in_samples[0]) / CLIP32 / RX1_GAIN > 0.1) { input_samples = 0; state++; } break; case 12: if (input_samples >= MILLISECS(200)) { state++; } else { for (i = 0; i < nSamples; i++) { in_sample = amp_in_samples[i] / CLIP32; mag_in = cabs(in_sample); if (tx_to_RX1_gain < mag_in) tx_to_RX1_gain = mag_in; out_sample = amp_out_samples[i] / CLIP32; mag_out = cabs(out_sample); if (max_out_mag < mag_out) max_out_mag = mag_out; } } break; case 13: #if DEBUG_PS printf("tx RX1 gain %.5f, max_out_mag %.6f\n", tx_to_RX1_gain, max_out_mag); #endif if (fabs(tx_to_RX1_gain - RX1_GAIN) / RX1_GAIN > 1E-2) printf("Measured %.6f and assumed %.6f RX1 gain differ\n", tx_to_RX1_gain, RX1_GAIN); input_samples = 0; state++; break; case 14: input_samples = 0; state++; break; case 19: quisk_set_key_down(0); tx_state = NothingHappens; PsCal = 0; if ( ! Fault) isCalibrated = 1; input_samples = 0; state = 0; break; } }