DD4WH
Well-known member
Frequency shift with the Weaver method ("third method") using quadrature oscillators
I have always been planning to use the Weaver method of SSB demodulation for a software defined radio. The Weaver method has become famous as the "third method" in the analog ham radio scene.The new Teensy 3.5 with its floating point unit makes this possible.
The method itself is a little complex and others can surely explain it much better than I can.
http://csoundjournal.com/ezine/summer2000/processing/
http://www.pa3ect.eu/start/weaver-derde-methode-ssb-visueel-uitgelegd/
Summers, H. (2015): Weaver article library. – online: http://www.hanssummers.com/weaver/weaverlib [2015-12-06]
Weaver, D. (1956): A Third Method of Generation of Single-Sideband Signals. – Proceedings of the IRE, Dec 1956.
It can be used for several purposes:
- frequency shifting: adding, not multiplying! Normal frequency shifting preserves the harmonics, because frequencies are multiplied. The frequency shifting by the Weaver method adds and therefore does not preserve the harmonic structure. Sounds strange, maybe useful for audio synthesizer freaks ;-)
- SSB (single sideband) demodulation with very good mirror rejection (no FIR Hilberts needed)
- frequency shifting in software defined radios by shifting one of the quadrature oscillators instead of tuning the hardware local oscillator, could be used for passband tuning
- etc. . . .
I have not been able to find code in C on the web for this method, so I tried to implement it by combining code of others on the Teensy, thanks guys! Thanks also to Pete for help in optimizing the floating point queue code!
It uses a mono input (left LINE channel). The 1st quadrature oscillator produces two signals with a frequency of (sampling_rate / 16), I & Q. Q is 90 degrees out of phase. These two signals are then lowpass-filtered by two identical IIR filters. The second quadrature oscillator has a frequency of (sampling rate/16) plus a desired offset (the amount of frequency we want to shift our original signal) and again has two outputs which are 90 degrees apart. The lowpass-filtered I & Q signals are mixed with the two outputs of the second quadrature oscillator. After that we simply add I & Q to get our original mono signal shifted by the desired offset frequency! I had to implement two types of quadrature oscillators, one efficient one, and one processor-intensive one (see the code for the reasons why). The whole audio chain works in 44.1kHz without decimation/interpolation and is implemented in floating point and needs 440µsec to process a block of 128 samples.
Here is the script for the Teensy, tested on the 3.5 with the new CMSIS lib V4.5. Try setting NCO_FREQ = AUDIO_SAMPLE_RATE_EXACT / 16 + 20; that will twist your ears ;-).
It also has a little spectrum analyser for an ILI9341 and lots of filters that could be switched in, but that is not necessary for the Weaver audio chain, but helps in making experiments with filters.
Suggestions for modifications and optimisations welcome!
Frank
Filter_coeffs.h
I have always been planning to use the Weaver method of SSB demodulation for a software defined radio. The Weaver method has become famous as the "third method" in the analog ham radio scene.The new Teensy 3.5 with its floating point unit makes this possible.
The method itself is a little complex and others can surely explain it much better than I can.
http://csoundjournal.com/ezine/summer2000/processing/
http://www.pa3ect.eu/start/weaver-derde-methode-ssb-visueel-uitgelegd/
Summers, H. (2015): Weaver article library. – online: http://www.hanssummers.com/weaver/weaverlib [2015-12-06]
Weaver, D. (1956): A Third Method of Generation of Single-Sideband Signals. – Proceedings of the IRE, Dec 1956.
It can be used for several purposes:
- frequency shifting: adding, not multiplying! Normal frequency shifting preserves the harmonics, because frequencies are multiplied. The frequency shifting by the Weaver method adds and therefore does not preserve the harmonic structure. Sounds strange, maybe useful for audio synthesizer freaks ;-)
- SSB (single sideband) demodulation with very good mirror rejection (no FIR Hilberts needed)
- frequency shifting in software defined radios by shifting one of the quadrature oscillators instead of tuning the hardware local oscillator, could be used for passband tuning
- etc. . . .
I have not been able to find code in C on the web for this method, so I tried to implement it by combining code of others on the Teensy, thanks guys! Thanks also to Pete for help in optimizing the floating point queue code!
It uses a mono input (left LINE channel). The 1st quadrature oscillator produces two signals with a frequency of (sampling_rate / 16), I & Q. Q is 90 degrees out of phase. These two signals are then lowpass-filtered by two identical IIR filters. The second quadrature oscillator has a frequency of (sampling rate/16) plus a desired offset (the amount of frequency we want to shift our original signal) and again has two outputs which are 90 degrees apart. The lowpass-filtered I & Q signals are mixed with the two outputs of the second quadrature oscillator. After that we simply add I & Q to get our original mono signal shifted by the desired offset frequency! I had to implement two types of quadrature oscillators, one efficient one, and one processor-intensive one (see the code for the reasons why). The whole audio chain works in 44.1kHz without decimation/interpolation and is implemented in floating point and needs 440µsec to process a block of 128 samples.
Here is the script for the Teensy, tested on the 3.5 with the new CMSIS lib V4.5. Try setting NCO_FREQ = AUDIO_SAMPLE_RATE_EXACT / 16 + 20; that will twist your ears ;-).
It also has a little spectrum analyser for an ILI9341 and lots of filters that could be switched in, but that is not necessary for the Weaver audio chain, but helps in making experiments with filters.
Suggestions for modifications and optimisations welcome!
Frank
Code:
/***********************************************************************
*
* Weaver-style audio chain
*
* Frank DD4WH 2016_10_30
*
* first experiment
*/
// with optimizations by Pete El Supremo 2016_10_27, thanks Pete!
//
// this is a setup to test whether we can do real time audio processing
// with the Teensy 3.5 in floating point math
// it takes the LINE INPUT audio, converts it from int16_t to float32_t,
// does the audio processing in float32_t with the NEW ARM CMSIS lib (see https://forum.pjrc.com/threads/38325-Excellent-results-with-Floating-Point-FFT-IFFT-Processing-and-Teensy-3-6?p=119797&viewfull=1#post119797),
// and see here: https://forum.pjrc.com/threads/38325-Excellent-results-with-Floating-Point-FFT-IFFT-Processing-and-Teensy-3-6?p=120177&viewfull=1#post120177
// . . . converts the audio back to int16_t and gives it back to the headphone output
//
// MIT licence
#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <SD.h>
#include <Metro.h>
#include "font_Arial.h"
#include <ILI9341_t3.h>
#include <arm_math.h>
#include <arm_const_structs.h>
#include "Filter_coeffs.h"
#define BACKLIGHT_PIN 0
#define TFT_DC 20
#define TFT_CS 21
#define TFT_RST 32 // 255 = unused. connect to 3.3V
#define TFT_MOSI 7
#define TFT_SCLK 14
#define TFT_MISO 12
//void sinewave(void);
ILI9341_t3 tft = ILI9341_t3(TFT_CS, TFT_DC, TFT_RST, TFT_MOSI, TFT_SCLK, TFT_MISO);
Metro five_sec=Metro(2000); // Set up a 0.5 second Metro
// this audio comes from the codec by I2S2
AudioInputI2S i2s_in;
AudioRecordQueue Q_in_L;
AudioRecordQueue Q_in_R;
AudioPlayQueue Q_out_L;
AudioPlayQueue Q_out_R;
AudioAnalyzeFFT256 myFFT;
AudioOutputI2S i2s_out;
AudioConnection patchCord1(i2s_in, 0, Q_in_L, 0);
AudioConnection patchCord2(i2s_in, 1, Q_in_R, 0);
AudioConnection patchCord5(Q_out_R,0,myFFT,0);
AudioConnection patchCord3(Q_out_L, 0, i2s_out, 1);
AudioConnection patchCord4(Q_out_R, 0, i2s_out, 0);
AudioControlSGTL5000 sgtl5000_1; //xy=265.212
int idx_t = 0;
int idx = 0;
int64_t sum;
float32_t mean;
int n_L;
int n_R;
long int n_clear;
int peak[512];
int barm[512];
ulong samp_ptr = 0;
bool FFT_state = false;
const int myInput = AUDIO_INPUT_LINEIN;
// We're only processing one buffer at a time so the
// number of samples is fixed at 128
#define BUFFER_SIZE 128
// buffers for quadrature oscillator 1
float32_t Osc1_Q_buffer [BUFFER_SIZE];
float32_t Osc1_I_buffer [BUFFER_SIZE];
float32_t a1_buffer [BUFFER_SIZE];
float32_t b1_buffer [BUFFER_SIZE];
float32_t c1_buffer [BUFFER_SIZE];
float32_t d1_buffer [BUFFER_SIZE];
float32_t e1_buffer [BUFFER_SIZE];
float32_t f1_buffer [BUFFER_SIZE];
bool sine_flag1 = false;
int16_t IF_FREQ1 = AUDIO_SAMPLE_RATE_EXACT / 16;
bool dir1 = false;
// definitions for quadrature oscillator 2
float32_t NCO_FREQ = AUDIO_SAMPLE_RATE_EXACT / 16; // + 20;
float32_t NCO_INC = 2 * PI * NCO_FREQ / AUDIO_SAMPLE_RATE_EXACT;
float32_t OSC_COS = cos (NCO_INC);
float32_t OSC_SIN = sin (NCO_INC);
float32_t Osc_Vect_Q = 1.0;
float32_t Osc_Vect_I = 0.0;
float32_t Osc_Gain = 0.0;
float32_t Osc_Q = 0.0;
float32_t Osc_I = 0.0;
float32_t i_temp = 0.0;
float32_t q_temp = 0.0;
//float32_t Osc2_Q_buffer [BUFFER_SIZE];
//float32_t Osc2_I_buffer [BUFFER_SIZE];
/*
bool sine_flag2 = false;
//int16_t IF_FREQ2 = 5515;
int16_t IF_FREQ2 = AUDIO_SAMPLE_RATE_EXACT / 16;
bool dir2 = true;
*/
// Decimation and Interpolation factor
#define DECIMATION_FACTOR 4
float32_t float_buffer_L [BUFFER_SIZE];
float32_t float_buffer_R [BUFFER_SIZE];
float32_t float_buffer_L_3 [BUFFER_SIZE];
float32_t float_buffer_R_3 [BUFFER_SIZE];
// This holds the samples output by the decimation
// and therefore there are a factor of 4 fewer samples
float32_t float_buffer_L_2 [BUFFER_SIZE / DECIMATION_FACTOR];
float32_t float_buffer_R_2 [BUFFER_SIZE / DECIMATION_FACTOR];
// determine the number of decimation and interpolation coefficients
// from the size of their arrays
const int N_DEC_COEFFS = sizeof(FIR_dec1_coeffs)/sizeof(FIR_dec1_coeffs[0]);
const int N_INT_COEFFS = sizeof(FIR_int1_coeffs)/sizeof(FIR_int1_coeffs[0]);
// decimation with FIR lowpass
arm_fir_decimate_instance_f32 FIR_dec1_L;
float32_t FIR_decim_state_L [N_DEC_COEFFS + BUFFER_SIZE - 1];
arm_fir_decimate_instance_f32 FIR_dec1_R;
float32_t FIR_decim_state_R [N_DEC_COEFFS + BUFFER_SIZE - 1];
// interpolation with FIR lowpass
arm_fir_interpolate_instance_f32 FIR_int1_L;
float32_t FIR_interp_state_L [(N_INT_COEFFS / DECIMATION_FACTOR) + BUFFER_SIZE - 1];
arm_fir_interpolate_instance_f32 FIR_int1_R;
float32_t FIR_interp_state_R [(N_INT_COEFFS / DECIMATION_FACTOR) + BUFFER_SIZE - 1];
// FIR filter
// determine the number of filter coeffs from the size of their arrays
const int N_FIR1_COEFFS = sizeof(FIR_lowpass1_coeffs)/sizeof(FIR_lowpass1_coeffs[0]);
arm_fir_instance_f32 FIR_lowpass1_L;
float32_t FIR_lowpass1_state_L [N_FIR1_COEFFS + BUFFER_SIZE - 1];
arm_fir_instance_f32 FIR_lowpass1_R;
float32_t FIR_lowpass1_state_R [N_FIR1_COEFFS + BUFFER_SIZE - 1];
/****************************************************************************************
* init IIR filters
****************************************************************************************/
// 2-pole biquad IIR - definitions and Initialisation
const int N_stages_biquad_lowpass1 = sizeof(biquad_lowpass1_coeffs)/sizeof(biquad_lowpass1_coeffs[0]) / 5;
float32_t biquad_lowpass1_state_L [N_stages_biquad_lowpass1 * 4];
float32_t biquad_lowpass1_state_R [N_stages_biquad_lowpass1 * 4];
arm_biquad_casd_df1_inst_f32 biquad_lowpass1_L = {N_stages_biquad_lowpass1, biquad_lowpass1_state_L, biquad_lowpass1_coeffs};
arm_biquad_casd_df1_inst_f32 biquad_lowpass1_R = {N_stages_biquad_lowpass1, biquad_lowpass1_state_R, biquad_lowpass1_coeffs};
// 10-pole biquad IIR - definitions and Initialisation
const int N_stages_biquad_lowpass2 = sizeof(biquad_lowpass2_coeffs)/sizeof(biquad_lowpass2_coeffs[0]) / 5;
float32_t biquad_lowpass2_state_L [N_stages_biquad_lowpass2 * 4];
float32_t biquad_lowpass2_state_R [N_stages_biquad_lowpass2 * 4];
arm_biquad_casd_df1_inst_f32 biquad_lowpass2_L = {N_stages_biquad_lowpass2, biquad_lowpass2_state_L, biquad_lowpass2_coeffs};
arm_biquad_casd_df1_inst_f32 biquad_lowpass2_R = {N_stages_biquad_lowpass2, biquad_lowpass2_state_R, biquad_lowpass2_coeffs};
// complex FFT with the new library CMSIS V4.5
const static arm_cfft_instance_f32 *S;
const int length_FFT = 256;
float32_t FFT_buffer [length_FFT * 2] __attribute__ ((aligned (4)));
void setup() {
Serial.begin(115200);
delay(1000);
// Audio connections require memory. and the record queue
// uses this memory to buffer incoming audio.
AudioMemory(100);
// Enable the audio shield. select input. and enable output
sgtl5000_1.enable();
sgtl5000_1.inputSelect(myInput);
sgtl5000_1.volume(0.5);
sgtl5000_1.adcHighPassFilterDisable(); // does not help too much!
/* // Initialize the SD card
SPI.setMOSI(7);
SPI.setSCK(14);
if (!(SD.begin(10))) {
// stop here if no SD card. but print a message
while (1) {
Serial.println("Unable to access the SD card");
delay(500);
}
}
*/
pinMode( BACKLIGHT_PIN, OUTPUT );
analogWrite( BACKLIGHT_PIN, 1023 );
tft.begin();
tft.setRotation( 3 );
tft.fillScreen(ILI9341_BLACK);
tft.setCursor(10, 1);
tft.setTextSize(2);
tft.setTextColor(ILI9341_WHITE);
tft.setFont(Arial_14);
tft.print("Floating point audio processing");
/****************************************************************************************
* init FIR filters
****************************************************************************************/
arm_fir_init_f32(&FIR_lowpass1_L, N_FIR1_COEFFS, FIR_lowpass1_coeffs, FIR_lowpass1_state_L,BUFFER_SIZE);
arm_fir_init_f32(&FIR_lowpass1_R, N_FIR1_COEFFS, FIR_lowpass1_coeffs, FIR_lowpass1_state_R,BUFFER_SIZE);
/****************************************************************************************
* init decimation and interpolation filters
****************************************************************************************/
if(arm_fir_decimate_init_f32(&FIR_dec1_L, N_DEC_COEFFS, DECIMATION_FACTOR, FIR_dec1_coeffs, FIR_decim_state_L, BUFFER_SIZE)) {
Serial.println("Init of decimation failed");
while(1);
}
if(arm_fir_interpolate_init_f32(&FIR_int1_L, DECIMATION_FACTOR, N_INT_COEFFS, FIR_int1_coeffs, FIR_interp_state_L, BUFFER_SIZE/DECIMATION_FACTOR)) {
Serial.println("Init of interpolation failed");
while(1);
}
if(arm_fir_decimate_init_f32(&FIR_dec1_R, N_DEC_COEFFS, DECIMATION_FACTOR, FIR_dec1_coeffs, FIR_decim_state_R, BUFFER_SIZE)) {
Serial.println("Init of decimation failed");
while(1);
}
if(arm_fir_interpolate_init_f32(&FIR_int1_R, DECIMATION_FACTOR, N_INT_COEFFS, FIR_int1_coeffs, FIR_interp_state_R, BUFFER_SIZE/DECIMATION_FACTOR)) {
Serial.println("Init of interpolation failed");
while(1);
}
/****************************************************************************************
* init complex FFT
****************************************************************************************/
switch (length_FFT) {
case 16:
S = &arm_cfft_sR_f32_len16;
break;
case 32:
S = &arm_cfft_sR_f32_len32;
break;
case 64:
S = &arm_cfft_sR_f32_len64;
break;
case 128:
S = &arm_cfft_sR_f32_len128;
break;
case 256:
S = &arm_cfft_sR_f32_len256;
break;
case 512:
S = &arm_cfft_sR_f32_len512;
break;
case 1024:
S = &arm_cfft_sR_f32_len1024;
break;
case 2048:
S = &arm_cfft_sR_f32_len2048;
break;
case 4096:
S = &arm_cfft_sR_f32_len4096;
break;
}
/****************************************************************************************
* begin to queue the audio from the audio library
****************************************************************************************/
delay(100);
Q_in_L.begin();
Q_in_R.begin();
} // END SETUP
int16_t *sp_L;
int16_t *sp_R;
void loop() {
elapsedMicros usec = 0;
/**********************************************************************************
* Get samples from queue buffers
**********************************************************************************/
// this is supposed to prevent overfilled queue buffers
if (Q_in_L.available() > 3 || Q_in_R.available() > 3) {
Q_in_L.clear();
Q_in_R.clear();
n_clear ++; // just for debugging to check how often this occurs
}
// is there at least one buffer in each channel available ?
if (Q_in_L.available() >= 1 && Q_in_R.available() >= 1)
{
sp_L = Q_in_L.readBuffer();
sp_R = Q_in_R.readBuffer();
// convert to float
arm_q15_to_float (sp_L, float_buffer_L, BUFFER_SIZE); // convert int_buffer to float 32bit
arm_q15_to_float (sp_L, float_buffer_R, BUFFER_SIZE); // convert int_buffer to float 32bit
Q_in_L.freeBuffer();
Q_in_R.freeBuffer();
/**********************************************************************************
* Put 128 floating point samples into FFT buffer and set flag FFT_state when it is filled
**********************************************************************************/
for(int i = 0; i < BUFFER_SIZE; i++)
{
if(FFT_state == false) // FFT buffer not yet filled
{
FFT_buffer [samp_ptr] = (float32_t)float_buffer_L[i]; // get floating point data for FFT for spectrum scope/waterfall display
samp_ptr++;
FFT_buffer [samp_ptr] = (float32_t)float_buffer_R[i];
samp_ptr++;
// On obtaining enough samples for spectrum scope/waterfall, update state machine, reset pointer and wait until we process what we have
if(samp_ptr > length_FFT)
{
samp_ptr = 0;
FFT_state = true; // FFT buffer filled with length_FFT samples
}
}
}
/**************************************************************************
* From here, all the 32 bit float audio processing can start
* ************************************************************************
/**************************************************************************
* This is the Weaver audio chain implemented by two quadrature oscillators
* and two IIR lowpass filters
* For more info see:
* http://csoundjournal.com/ezine/summer2000/processing/
* http://www.pa3ect.eu/start/weaver-derde-methode-ssb-visueel-uitgelegd/
* Summers, H. (2015): Weaver article library. – online: http://www.hanssummers.com/weaver/weaverlib [2015-12-06]
* Weaver, D. (1956): A Third Method of Generation of Single-Sideband Signals. – Proceedings of the IRE, Dec 1956.
* *************************************************************************/
// take the mono input signal and mix it with the I & Q outputs of a quadrature oscillator
// running at (sampling frequency / 16) = 2757.4Hz
freq_conv1(); // complex shift by fs/16
// I & Q are now filtered with a lowpass filter with Fstop = 2750Hz --> which will lead to a Bandwidth of 2 x 2750Hz = 5.5kHz
// The lowpass filters are identical IIR biquad filters with 5 stages = 10th order
arm_biquad_cascade_df1_f32 (&biquad_lowpass2_L, float_buffer_L,float_buffer_L_3, BUFFER_SIZE);
arm_biquad_cascade_df1_f32 (&biquad_lowpass2_R, float_buffer_R,float_buffer_R_3, BUFFER_SIZE);
// The second quadrature oscillator runs at (sampling frequency / 16) PLUS the desired offset frequency
// the two outputs of the oscillator (which are 90 degrees apart) are complex multiplied with I & Q
freq_conv2();
// We then take I & Q and simply add them together, the result is in float_buffer_R and is copied to float_buffer_L for mono output
arm_add_f32(float_buffer_R, float_buffer_L, float_buffer_R, BUFFER_SIZE); // add I & Q
arm_copy_f32 (float_buffer_R, float_buffer_L, BUFFER_SIZE); // copy to left audio chain --> Mono output of shifted signal
// THAT WAS THE WEAVER CHAIN !
// potential applications are:
// - Frequency shifting for passband tuning
// - SSB demodulation in Software defined radios without the need for phase shifting Hilbert transforms
// -
/*
// decimation by 4 --> 44118 / 4 = 11029.5 sps, means that a 4.2kHz decimation filter is fine [should have -80dB at 5.5kHz]
// decimation filter FIR 80 taps, lowpass 4.2kHz, Parks McClellan, Window off, Transition width 0.1, 80dB stopband attenuation
arm_fir_decimate_f32(&FIR_dec1_L, float_buffer_L, float_buffer_L_2, BUFFER_SIZE);
arm_fir_decimate_f32(&FIR_dec1_R, float_buffer_R, float_buffer_R_2, BUFFER_SIZE);
// test filter 1 stage
// arm_biquad_cascade_df1_f32 (&biquad_lowpass1_L, float_buffer_L_2,float_buffer_L_3, BUFFER_SIZE / DECIMATION_FACTOR);
// arm_biquad_cascade_df1_f32 (&biquad_lowpass1_R, float_buffer_R_2,float_buffer_R_3, BUFFER_SIZE / DECIMATION_FACTOR);
// test filter 5 to 8 stages
// arm_biquad_cascade_df1_f32 (&biquad_lowpass2_L, float_buffer_L_3,float_buffer_L_2, BUFFER_SIZE / DECIMATION_FACTOR);
// arm_biquad_cascade_df1_f32 (&biquad_lowpass2_R, float_buffer_R_3,float_buffer_R_2, BUFFER_SIZE / DECIMATION_FACTOR);
// test FIR filter 200 taps
// arm_fir_f32(&FIR_lowpass1_L,float_buffer_L_3, float_buffer_L_2, BUFFER_SIZE / DECIMATION_FACTOR);
// arm_fir_f32(&FIR_lowpass1_R,float_buffer_R_3, float_buffer_R_2, BUFFER_SIZE / DECIMATION_FACTOR);
// interpolation by 4
arm_fir_interpolate_f32(&FIR_int1_L, float_buffer_L_2, float_buffer_L, BUFFER_SIZE / DECIMATION_FACTOR);
arm_fir_interpolate_f32(&FIR_int1_R, float_buffer_R_2, float_buffer_R, BUFFER_SIZE / DECIMATION_FACTOR);
// this IIR filter works in 44118sps ! Quite expensive in terms of CPU power
// arm_biquad_cascade_df1_f32 (&biquad_lowpass2_L, float_buffer_L_3,float_buffer_L, BUFFER_SIZE);
// arm_biquad_cascade_df1_f32 (&biquad_lowpass2_R, float_buffer_R_3,float_buffer_R, BUFFER_SIZE);
// scaling after interpolation ? multiply volume by DECIMATION_FACTOR seems to be a reasonable figure (output = input volume)
arm_scale_f32 (float_buffer_R,(float32_t) DECIMATION_FACTOR,float_buffer_R, BUFFER_SIZE);
arm_scale_f32 (float_buffer_L,(float32_t) DECIMATION_FACTOR,float_buffer_L, BUFFER_SIZE);
*/
/**************************************************************************
* END of 32 bit float audio processing
* ************************************************************************
*/
sp_L = Q_out_L.getBuffer();
sp_R = Q_out_R.getBuffer();
arm_float_to_q15 (float_buffer_L, sp_L, BUFFER_SIZE);
arm_float_to_q15 (float_buffer_R, sp_R, BUFFER_SIZE);
Q_out_L.playBuffer(); // play it !
Q_out_R.playBuffer(); // play it !
/**********************************************************************************
* FFT
**********************************************************************************/
if (FFT_state) {
// arm_cfft_f32(S, FFT_buffer, 0, 1);
}
/**********************************************************************************
* PRINT ROUTINE FOR ELAPSED MICROSECONDS
**********************************************************************************/
sum = sum + usec;
idx_t++;
if (idx_t > 1000) {
tft.fillRect(240,50,90,20,ILI9341_BLACK);
tft.setCursor(240, 50);
mean = sum / idx_t;
tft.print (mean);
Serial.print (mean);
Serial.print (" microsec for 2 stereo blocks ");
Serial.println();
idx_t = 0;
sum = 0;
}
}
/**********************************************************************************
* PRINT ROUTINE FOR AUDIO LIBRARY PROCESSOR AND MEMORY USAGE
**********************************************************************************/
if (five_sec.check() == 1)
{
Serial.print("Proc = ");
Serial.print(AudioProcessorUsage());
Serial.print(" (");
Serial.print(AudioProcessorUsageMax());
Serial.print("), Mem = ");
Serial.print(AudioMemoryUsage());
Serial.print(" (");
Serial.print(AudioMemoryUsageMax());
Serial.println(")");
Serial.print("Cleared the audio buffer ");
Serial.print(n_clear); Serial.println (" times. ");
/* tft.fillRect(100,120,200,80,ILI9341_BLACK);
tft.setCursor(10, 120);
tft.setTextSize(2);
tft.setTextColor(ILI9341_WHITE);
tft.setFont(Arial_14);
tft.print ("Proc = ");
tft.setCursor(100, 120);
tft.print (AudioProcessorUsage());
tft.setCursor(180, 120);
tft.print (AudioProcessorUsageMax());
tft.setCursor(10, 150);
tft.print ("Mem = ");
tft.setCursor(100, 150);
tft.print (AudioMemoryUsage());
tft.setCursor(180, 150);
tft.print (AudioMemoryUsageMax());
*/
AudioProcessorUsageMaxReset();
AudioMemoryUsageMaxReset();
}
spectrum();
}
void spectrum() { // spectrum analyser code by rheslip - modified
if (myFFT.available()) {
int scale;
scale = 2;
for (int16_t x=2; x < 100; x+=2) {
int bar = (abs(myFFT.output[x]) * scale);
if (bar >180) bar=180;
// this is a very simple IIR filter to smooth the reaction of the bars
bar = 0.2 * bar + 0.8 * barm[x];
if (bar > peak[x]) peak[x]=bar;
// tft.drawFastVLine(x, 210-bar,bar, ILI9341_PURPLE);
tft.drawFastVLine(x*2+10, 210-bar,bar, ILI9341_PINK);
tft.drawFastVLine(x*2+10, 20, 210-bar-20, ILI9341_BLACK);
tft.drawPixel(x*2+10,209-peak[x], ILI9341_YELLOW);
if(peak[x]>0) peak[x]-=1;
barm[x] = bar;
}
} //end if
} // end void spectrum
/*************************************************************************************************
* FREQUENCY CONVERSION USING A SOFTWARE QUADRATURE OSCILLATOR
*
* THIS VERSION USES A PRECALCULATED COS AND SIN WAVE AND IS VERY FAST AND EFFICIENT
*
* MAJOR DRAWBACK: frequency conversion can only be done at sub-multiples of the sampling frequency
*
* large parts of the code taken from the mcHF code by Clint, KA7OEI, thank you!
* see here for more info on quadrature oscillators:
* Wheatley, M. (2011): CuteSDR Technical Manual Ver. 1.01. - http://sourceforge.net/projects/cutesdr/
* Lyons, R.G. (2011): Understanding Digital Processing. – Pearson, 3rd edition.
*************************************************************************************************/
/////////////////////////////////////////////////////////////////////
// Sine/cosine generation function
// Adjust rad_calc *= 32; p.e. ( 44100 / 128 * 32 = 11025 khz)
// this can only generate frequencies in sub-multiples of the SAMPLE_RATE !
/////////////////////////////////////////////////////////////////////
/**
*
*/
void freq_conv1()
{
uint i;
float32_t rad_calc;
if (!sine_flag1) { // have we already calculated the sine wave?
for (i = 0; i < BUFFER_SIZE; i++) { // No, let's do it!
rad_calc = (float32_t)i; // convert to float the current position within the buffer
rad_calc /= (BUFFER_SIZE); // make this a fraction
rad_calc *= (PI * 2); // convert to radians
rad_calc *= IF_FREQ1 * 128 / AUDIO_SAMPLE_RATE_EXACT; // multiply by number of cycles that we want within this block ( 44100 / 128 * 32 = 11025 khz)
//
Osc1_Q_buffer [i] = arm_cos_f32(rad_calc); // get sine and cosine values and store in pre-calculated array
Osc1_I_buffer [i] = arm_sin_f32(rad_calc); // they are in float32_t format
}
sine_flag1 = 1; // signal that once we have generated the quadrature sine waves, we shall not do it again
}
// Do frequency conversion using optimized ARM math functions [KA7OEI]
// there seems to be something wrong here: I is real! Q is imaginary
arm_mult_f32(float_buffer_L, Osc1_Q_buffer, c1_buffer, BUFFER_SIZE); // multiply products for converted I channel
arm_mult_f32(float_buffer_R, Osc1_I_buffer, d1_buffer, BUFFER_SIZE);
arm_mult_f32(float_buffer_L, Osc1_I_buffer, e1_buffer, BUFFER_SIZE);
arm_mult_f32(float_buffer_R, Osc1_Q_buffer, f1_buffer, BUFFER_SIZE); // multiply products for converted Q channel
if(!dir1) // Conversion is "above" on RX (LO needs to be set lower)
{
arm_add_f32(f1_buffer, e1_buffer, float_buffer_R, BUFFER_SIZE); // summation for I channel
arm_sub_f32(c1_buffer, d1_buffer, float_buffer_L, BUFFER_SIZE); // difference for Q channel
}
else // Conversion is "below" on RX (LO needs to be set higher)
{
arm_add_f32(c1_buffer, d1_buffer, float_buffer_L, BUFFER_SIZE); // summation for I channel
arm_sub_f32(f1_buffer, e1_buffer, float_buffer_R, BUFFER_SIZE); // difference for Q channel
}
} // end freq_conv1()
/*************************************************************************************************
* freq_conv2()
*
* FREQUENCY CONVERSION USING A SOFTWARE QUADRATURE OSCILLATOR
*
* THIS VERSION calculates the COS AND SIN WAVE on the fly AND IS SLOW AND INEFFICIENT
*
* MAJOR ADVANTAGE: frequency conversion can be done for any frequency !
*
* large parts of the code taken from the mcHF code by Clint, KA7OEI, thank you!
* see here for more info on quadrature oscillators:
* Wheatley, M. (2011): CuteSDR Technical Manual Ver. 1.01. - http://sourceforge.net/projects/cutesdr/
* Lyons, R.G. (2011): Understanding Digital Processing. – Pearson, 3rd edition.
*************************************************************************************************/
// INPUT I = float_buffer_L_3
// INPUT Q = float_buffer_R_3
// OUTPUT I = float_buffer_L
// OUTPUT Q = float_buffer_R
/*
void set_freq_conv2(float32_t NCO_FREQ) {
// float32_t NCO_FREQ = AUDIO_SAMPLE_RATE_EXACT / 16; // + 20;
float32_t NCO_INC = 2 * PI * NCO_FREQ / AUDIO_SAMPLE_RATE_EXACT;
float32_t OSC_COS = cos (NCO_INC);
float32_t OSC_SIN = sin (NCO_INC);
float32_t Osc_Vect_Q = 1.0;
float32_t Osc_Vect_I = 0.0;
float32_t Osc_Gain = 0.0;
float32_t Osc_Q = 0.0;
float32_t Osc_I = 0.0;
float32_t i_temp = 0.0;
float32_t q_temp = 0.0;
}
*/
void freq_conv2()
{
uint i;
for(i = 0; i < BUFFER_SIZE; i++) {
// generate local oscillator on-the-fly: This takes a lot of processor time!
Osc_Q = (Osc_Vect_Q * OSC_COS) - (Osc_Vect_I * OSC_SIN); // Q channel of oscillator
Osc_I = (Osc_Vect_I * OSC_COS) + (Osc_Vect_Q * OSC_SIN); // I channel of oscillator
Osc_Gain = 1.95 - ((Osc_Vect_Q * Osc_Vect_Q) + (Osc_Vect_I * Osc_Vect_I)); // Amplitude control of oscillator
// rotate vectors while maintaining constant oscillator amplitude
Osc_Vect_Q = Osc_Gain * Osc_Q;
Osc_Vect_I = Osc_Gain * Osc_I;
//
// do actual frequency conversion
float_buffer_L[i] = (float_buffer_L_3[i] * Osc_Q) + (float_buffer_R_3[i] * Osc_I); // multiply I/Q data by sine/cosine data to do translation
float_buffer_R[i] = (float_buffer_R_3[i] * Osc_Q) - (float_buffer_L_3[i] * Osc_I);
//
}
} // end freq_conv2()
Filter_coeffs.h
Code:
// pass-thru coefficients
float32_t biquad_lowpass1_coeffs[5] = {1,0,0,0,0};
// lowpass elliptic 2.75kHz, IIR biquad 5 stages = 10 pole
// fs 44118Hz
// IIR Filter designer Iowa Hills
// a1 and a2 negated
// order of coefficients: b0, b1, b2, -a1, -a2
float32_t biquad_lowpass2_coeffs [25]= {
0.171265552838014729,
-0.306003398567613216,
0.171265552838014729,
1.650009424162684810,
-0.686537131271100942,
0.282254516360715746,
-0.496801858839390431,
0.282254516360715746,
1.692471910721256470,
-0.760179084603297528,
0.327140207187691878,
-0.548411686425261835,
0.327140207187691878,
1.746570461822357200,
-0.852439189772479011,
0.215776470640781259,
-0.297830800920305616,
0.215776470640781259,
1.791885744998123590,
-0.925607885359380322,
0.063661755193962916,
0.021017226351199236,
0.063661755193962916,
1.829338531467033620,
-0.977679268206158580
};
/*
// lowpass elliptic 4.0kHz, IIR biquad 6 stages = 12 pole
// fs 11027Hz
// IIR Filter designer Iowa Hills
// a1 and a2 negated
// order of coefficients: b0, b1, b2, -a1, -a2
float32_t biquad_lowpass2_coeffs [30]= {
0.390955637301552750,
0.556640362257148746,
0.390955637301552750,
-0.283754748892188324,
-0.054796887968065991,
0.539269449233711784,
0.786587988954648765,
0.539269449233711784,
-0.561224440037602634,
-0.303902447384469809,
0.691524236619742827,
1.062809822558603120,
0.691524236619742827,
-0.866733199372850538,
-0.579125096425238239,
0.776938719898112917,
1.296485971846480290,
0.776938719898112917,
-1.078306368085572900,
-0.772057043557133338,
0.805342431114359658,
1.482472831743180340,
0.805342431114359658,
-1.202727577982673820,
-0.890430115989226167,
0.814418978226675527,
1.611693182699644120,
0.814418978226675527,
-1.273176548353518900,
-0.967354590799476166
};
*/
// FIR 200 taps, Raised Cosine 0.940
// Fc = 3.000kHz, 75dB stopband
// fs 44118Hz
// just to test if this works,
// a 200 tap FIR is ridicously big and lots of calculation work for the processor
// --> it works !
float32_t FIR_lowpass1_coeffs[200] =
{ 94.97870007611338390E-9,
-1.009924021425136600E-6,
-3.955509751038522200E-6,
-8.167421170036924140E-6,
-12.23385614931934380E-6,
-14.14794389888045070E-6,
-11.83601267362507410E-6,
-3.874707993346245160E-6,
9.784409190338859470E-6,
27.31190087730911390E-6,
44.93736028076724410E-6,
57.46570395053866780E-6,
59.37446254570318160E-6,
46.32396853632215540E-6,
16.75260847587800940E-6,
-26.87322881082111080E-6,
-77.58971029118251290E-6,
-124.5675198379420860E-6,
-154.9396363165045050E-6,
-156.6193065484440300E-6,
-121.5946328111847800E-6,
-48.95032523151913040E-6,
53.22076923312532420E-6,
167.7426457861406560E-6,
270.5249246949471740E-6,
334.8987951214899680E-6,
337.4870793426350130E-6,
264.4733746227138910E-6,
116.8225367007839280E-6,
-87.00810636093507360E-6,
-311.9857878058807610E-6,
-511.4939298316222110E-6,
-635.9492759051854590E-6,
-643.6567497319431370E-6,
-511.7505262783362240E-6,
-244.6910091391752930E-6,
122.0563006080398620E-6,
524.7970382117916870E-6,
881.2433438597264510E-6,
0.001105897264525711,
0.001128584366403902,
912.4295155936849820E-6,
467.1426447303733770E-6,
-146.0617101058849410E-6,
-820.2329826975067140E-6,
-0.001419395225860928,
-0.001804044124385529,
-0.001861021146414526,
-0.001531785067380846,
-832.6535950832496840E-6,
138.4931805272558170E-6,
0.001212301918764672,
0.002174825940792150,
0.002807530450404730,
0.002933442142689424,
0.002460262925219348,
0.001410796432743954,
-67.40198731776092700E-6,
-0.001717749520116288,
-0.003215444385842726,
-0.004227679123662836,
-0.004482787665470263,
-0.003834751661105775,
-0.002308921085773049,
-117.1613956177577340E-6,
0.002364261142304478,
0.004653686825192925,
0.006251853354833492,
0.006743484329240796,
0.005893195184744424,
0.003715295847640134,
500.2318515875483630E-6,
-0.003212572302564382,
-0.006715879276971413,
-0.009259324885670151,
-0.010197256146698041,
-0.009132578886926748,
-0.006027440073467311,
-0.001255112944564374,
0.004423632349767859,
0.009963382868447519,
0.014204955965723772,
0.016088984764790849,
0.014876395277668464,
0.010335310034903742,
0.002856460155015252,
-0.006531028294970668,
-0.016255904142096857,
-0.024406432741385542,
-0.029011490508947704,
-0.028359336686767097,
-0.021302866903510787,
-0.007499202509221305,
0.012460426961983004,
0.037058421514990315,
0.064002438051746699,
0.090503604392753637,
0.113642557166307487,
0.130768970678281582,
0.139873389683470517,
0.139873389683470517,
0.130768970678281582,
0.113642557166307487,
0.090503604392753637,
0.064002438051746699,
0.037058421514990315,
0.012460426961983004,
-0.007499202509221305,
-0.021302866903510787,
-0.028359336686767097,
-0.029011490508947704,
-0.024406432741385542,
-0.016255904142096857,
-0.006531028294970668,
0.002856460155015252,
0.010335310034903742,
0.014876395277668464,
0.016088984764790849,
0.014204955965723772,
0.009963382868447519,
0.004423632349767859,
-0.001255112944564374,
-0.006027440073467311,
-0.009132578886926748,
-0.010197256146698041,
-0.009259324885670151,
-0.006715879276971413,
-0.003212572302564382,
500.2318515875483630E-6,
0.003715295847640134,
0.005893195184744424,
0.006743484329240796,
0.006251853354833492,
0.004653686825192925,
0.002364261142304478,
-117.1613956177577340E-6,
-0.002308921085773049,
-0.003834751661105775,
-0.004482787665470263,
-0.004227679123662836,
-0.003215444385842726,
-0.001717749520116288,
-67.40198731776092700E-6,
0.001410796432743954,
0.002460262925219348,
0.002933442142689424,
0.002807530450404730,
0.002174825940792150,
0.001212301918764672,
138.4931805272558170E-6,
-832.6535950832496840E-6,
-0.001531785067380846,
-0.001861021146414526,
-0.001804044124385529,
-0.001419395225860928,
-820.2329826975067140E-6,
-146.0617101058849410E-6,
467.1426447303733770E-6,
912.4295155936849820E-6,
0.001128584366403902,
0.001105897264525711,
881.2433438597264510E-6,
524.7970382117916870E-6,
122.0563006080398620E-6,
-244.6910091391752930E-6,
-511.7505262783362240E-6,
-643.6567497319431370E-6,
-635.9492759051854590E-6,
-511.4939298316222110E-6,
-311.9857878058807610E-6,
-87.00810636093507360E-6,
116.8225367007839280E-6,
264.4733746227138910E-6,
337.4870793426350130E-6,
334.8987951214899680E-6,
270.5249246949471740E-6,
167.7426457861406560E-6,
53.22076923312532420E-6,
-48.95032523151913040E-6,
-121.5946328111847800E-6,
-156.6193065484440300E-6,
-154.9396363165045050E-6,
-124.5675198379420860E-6,
-77.58971029118251290E-6,
-26.87322881082111080E-6,
16.75260847587800940E-6,
46.32396853632215540E-6,
59.37446254570318160E-6,
57.46570395053866780E-6,
44.93736028076724410E-6,
27.31190087730911390E-6,
9.784409190338859470E-6,
-3.874707993346245160E-6,
-11.83601267362507410E-6,
-14.14794389888045070E-6,
-12.23385614931934380E-6,
-8.167421170036924140E-6,
-3.955509751038522200E-6,
-1.009924021425136600E-6,
94.97870007611338390E-9
};
// decimation filter FIR 80 taps, lowpass 4.2kHz, Parks McClellan, Window off, Transition width 0.1, 80dB stopband attenuation
float32_t FIR_dec1_coeffs[80] = {
-142.9442739751772820E-6,
-252.2211831462206530E-6,
-347.5945623945052030E-6,
-321.0839601413753140E-6,
-97.22225175627986000E-6,
326.3117239250896090E-6,
836.1439112274239280E-6,
0.001209027523486800,
0.001184222537420025,
590.9136933192837660E-6,
-520.3527006113926060E-6,
-0.001812698722149656,
-0.002724775642735784,
-0.002673582490820696,
-0.001340728074926084,
0.001078011666779184,
0.003797965250984124,
0.005638404522781618,
0.005468470132752291,
0.002761334048239419,
-0.001980549872799683,
-0.007167492106103557,
-0.010557955026781203,
-0.010112507043340404,
-0.004978767523812532,
0.003808453095602158,
0.013295208276611466,
0.019416922705541947,
0.018489926838076525,
0.008904119903017756,
-0.007643552646951099,
-0.025966841359731561,
-0.038527883680910445,
-0.037684191047641612,
-0.018373491864758760,
0.019689814633387617,
0.071157875324433559,
0.125978118210888251,
0.171947745227487625,
0.198145667117833019,
0.198145667117833019,
0.171947745227487625,
0.125978118210888251,
0.071157875324433559,
0.019689814633387617,
-0.018373491864758760,
-0.037684191047641612,
-0.038527883680910445,
-0.025966841359731561,
-0.007643552646951099,
0.008904119903017756,
0.018489926838076525,
0.019416922705541947,
0.013295208276611466,
0.003808453095602158,
-0.004978767523812532,
-0.010112507043340404,
-0.010557955026781203,
-0.007167492106103557,
-0.001980549872799683,
0.002761334048239419,
0.005468470132752291,
0.005638404522781618,
0.003797965250984124,
0.001078011666779184,
-0.001340728074926084,
-0.002673582490820696,
-0.002724775642735784,
-0.001812698722149656,
-520.3527006113926060E-6,
590.9136933192837660E-6,
0.001184222537420025,
0.001209027523486800,
836.1439112274239280E-6,
326.3117239250896090E-6,
-97.22225175627986000E-6,
-321.0839601413753140E-6,
-347.5945623945052030E-6,
-252.2211831462206530E-6,
-142.9442739751772820E-6 };
float32_t FIR_int1_coeffs[80] = {
-142.9442739751772820E-6,
-252.2211831462206530E-6,
-347.5945623945052030E-6,
-321.0839601413753140E-6,
-97.22225175627986000E-6,
326.3117239250896090E-6,
836.1439112274239280E-6,
0.001209027523486800,
0.001184222537420025,
590.9136933192837660E-6,
-520.3527006113926060E-6,
-0.001812698722149656,
-0.002724775642735784,
-0.002673582490820696,
-0.001340728074926084,
0.001078011666779184,
0.003797965250984124,
0.005638404522781618,
0.005468470132752291,
0.002761334048239419,
-0.001980549872799683,
-0.007167492106103557,
-0.010557955026781203,
-0.010112507043340404,
-0.004978767523812532,
0.003808453095602158,
0.013295208276611466,
0.019416922705541947,
0.018489926838076525,
0.008904119903017756,
-0.007643552646951099,
-0.025966841359731561,
-0.038527883680910445,
-0.037684191047641612,
-0.018373491864758760,
0.019689814633387617,
0.071157875324433559,
0.125978118210888251,
0.171947745227487625,
0.198145667117833019,
0.198145667117833019,
0.171947745227487625,
0.125978118210888251,
0.071157875324433559,
0.019689814633387617,
-0.018373491864758760,
-0.037684191047641612,
-0.038527883680910445,
-0.025966841359731561,
-0.007643552646951099,
0.008904119903017756,
0.018489926838076525,
0.019416922705541947,
0.013295208276611466,
0.003808453095602158,
-0.004978767523812532,
-0.010112507043340404,
-0.010557955026781203,
-0.007167492106103557,
-0.001980549872799683,
0.002761334048239419,
0.005468470132752291,
0.005638404522781618,
0.003797965250984124,
0.001078011666779184,
-0.001340728074926084,
-0.002673582490820696,
-0.002724775642735784,
-0.001812698722149656,
-520.3527006113926060E-6,
590.9136933192837660E-6,
0.001184222537420025,
0.001209027523486800,
836.1439112274239280E-6,
326.3117239250896090E-6,
-97.22225175627986000E-6,
-321.0839601413753140E-6,
-347.5945623945052030E-6,
-252.2211831462206530E-6,
-142.9442739751772820E-6 };
/*
float32_t FIR_int1_coeffs[4] = {
0.210904123894329720,
0.301141108924676826,
0.301141108924676826,
0.210904123894329720 };
*/
/*
// brickwall CW filter
// bandpass, elliptic, Fc=700Hz, BW=300Hz, -80dB at 522Hz, -80dB at 932Hz
// 8 stages IIR
float32_t biquad_lowpass2_coeffs [40]= {
0.364001063989044915,
-0.696473959523549957,
0.364001063989044971,
1.788480069318890250,
-0.925376106125455067,
0.361186709947353413,
-0.620959940244166786,
0.361186709947353468,
1.751855271152122700,
-0.918221357851320086,
0.433058689492207660,
-0.831580060059418424,
0.433058689492207660,
1.831455406178845950,
-0.949172844759517020,
0.426669654838185863,
-0.723801208927776329,
0.426669654838185863,
1.738051696214021780,
-0.935169435187165887,
0.377153150378401081,
-0.731453681768477582,
0.377153150378401136,
1.864734872952909410,
-0.973104720052835992,
0.373007288989870733,
-0.599811328494536000,
0.373007288989870733,
1.743134620930931970,
-0.962407852528818108,
0.167662570128892602,
-0.332516665519374976,
0.167662570128892630,
1.886699807792520560,
-0.991809863590435659,
0.167030532646816859,
-0.144078663443187288,
0.167030532646816859,
1.757154994466575190,
-0.988071038589723227
};
*/