Pitch shift, frequency shifting, vocoder and granular effect for bird song listening

DD4WH

Well-known member
Many people suffer from high frequency hearing loss, especially people above the age of 50. Lang Elliott had the idea of lowering the pitch of bird songs with a special device in order to help bird enthusiasts with high frequency hearing loss to hear those higher-pitched bird songs [https://hearbirdsagain.org/]. These higher-pitched songs would normally be outside of their hearing ranges. His idea motivated me to try this on the Teensy 4.0. The plan is to experiment with different pitch shifting algorithms. I want to evaluate their specific suitability for pitch shifting of bird songs. [Bird songs differ a lot from the normal audio application range, i.e. speech or music !]

My scan of the Teensy forum revealed several existing libraries/algorithms for altering the pitch of audio sounds:


  1. Granular effect - in the Teensy Audio library
  2. Heterodyning, i.e. multiplying the audio with a sine wave of given frequency - implementable with the Audio Library
  3. Vocoder / pitch shifting in the frequency domain - sophisticated algorithm by Stephan Bernsee, implemented by Duff2013 for the Teensy
  4. Formant shifting in the frequency domain - sophisticated algorithm by Chip Audette, implemented in the Tympan library and later by Bob Larkin in the OpenAudio_Arduino_lib
  5. Frequency shifting by Quadrature Modulation - implemented by Bob Larkin in the OpenAudio library, originally intended for Software Defined Radio

So, after finding those implementations, I wanted to experiment with them in realtime by sending bird song audio from my smartphone to the Teensy and switch between those algorithms and tweak their parametres in order to find the best sounding solution for altering the pitch of bird songs by ear.

I use a Teensy 4.0 with an ADC PCM1808 and DAC5102a, but it should also run with a Teensy 4 and the Teensy audio board by using the Stereo line inputs.

The following sketch contains a real-time implementation of all the above mentioned methods, so I can switch and tweak them. Audio comes from the 3.5mm audio plug of my smartphone into the ADC and I can hear the result in my headphones. I implemented all methods in Stereo. Five buttons are used to toggle ON/OFF an input 2.5kHz highpass filter, to switch between the algorithms and to increase/decrease the amount of pitch shift of the specific method.

Code:
/*  
 *  PITCH SHIFTING, FORMANT SHIFTING AND FREQUENCY SHIFTING OF BIRD SOUNDS
 *  (c) Frank DD4WH
 *  
 *  Experiments with the Teensy 4.0 & ADC PCM1808 & DAC PCM5102
 *    
 *  to aid people with hearing loss to hear birds (again)
 *  
 *  inspired & motivated by Lang Elliott
 *  https://hearbirdsagain.org/  
 *  
 *  using the following code libraries:
 *  
 *  Teensy Audio Library: Paul Stoffregen
 *  
 *  Pitch Shifting using FFT: Stephan Bernsee & duff2013
 * 
 *  OpenAudioLibrary: Chip Audette & Bob Larkin 
 * 
 *  experimental platform
 *  all licenses remain those of the original code snippets that I used,
 *  mostly MIT license, but also some use GNU GPLv3 license
 *  
 *  July, 27th 2022
 */
 
/*
 *  Pitch Shifting:
 *  1.) advanced algorithm (smbPitchShift) modified by duff2013
 * 
 *  Heterodyning:
 *  2.) simple mixing of Audio signal with sine oscillator
 *  
 *  Granular Effect:
 *  3.) from Audio Library
 *  
 *  Frequency Shifting: 
 *  4.) Quadrature Modulator: designed for Software Defined Radios by Bob Larkin
 *  
 *  Formant shifting:
 *  5.) Tympan library by Chip Audette, version taken from OpenAudio lib
*/

/*
 * Only for Teensy 4.0 or Teensy 4.1 !!! This does not work with Teensy 3.2/3.5/3.6
 * Use the latest Arduino 1.8.19 and Teensyduino version 1.57
 * 
 * In order for this code to compile correctly, please install these libraries into your Arduino folder: 
 * 
 * 1.) https://github.com/duff2013/AudioEffectPhaseVocoder
 *     you have to correct line 26 of this file: effect_phaseVocoder.cpp --> should be: "#include "effect_phaseVocoder.h"" 
 * 2.) https://github.com/chipaudette/OpenAudio_ArduinoLibrary
 *     this is the source for Formant Shift and Quadrature modulation-based frequency shift
 * 3.) copy the file "hilbert251A.h" from folder OpenAudio_ArduinoLibrary/examples/FineFreqShift_OA to your sketch folder   
 * 
 */

/*  
Audio MEM INT16   Cur/Peak: 20/24
Audio MEM Float32 Cur/Peak: 24/26
CPU INT16 Cur/Peak: 60.39%/62.68%   
CPU F32   Cur/Peak: 15.60%/16.17% 
 *  
Memory Usage on Teensy 4.0:
   FLASH: code:98420, data:81032, headers:8960   free for files:1843204
   RAM1: variables:227104, code:95816, padding:2488   free for local variables:198880
   RAM2: variables:23808  free for malloc/new:500480

 */

#include "AudioStream_F32.h"
#include "OpenAudio_ArduinoLibrary.h"
#include "AudioEffectFormantShiftFD_OA_F32.h"  //the local file holding your custom function
#include <Audio.h>
#include <Bounce.h>
#include "effect_phaseVocoder.h"
//#include <Wire.h>
//#include <SPI.h>
//#include <SD.h>
//#include <SerialFlash.h>
// Filter for AudioFilter90Deg_F32 hilbert1
#include "hilbert251A.h"

// this is only needed for F32 OpenAudio lib !
//set the sample rate and block size
const float sample_rate_Hz = 44100.f; ; // other frequencies in the table in AudioOutputI2S_F32 for T3.x only
const int audio_block_samples = 128;     //for freq domain processing, a power of 2: 16, 32, 64, 128
AudioSettings_F32 audio_settings(sample_rate_Hz, audio_block_samples);

// GUItool: begin automatically generated code
AudioInputI2S            i2s1;           //xy=153,323
AudioAmplifier           ampL;           //xy=327,398
AudioAmplifier           ampR;           //xy=335,254
AudioFilterBiquad        highpassL;        //xy=479,350
AudioFilterBiquad        highpassR;        //xy=480,296
AudioMixer4              HPswitchR;         //xy=639,241
AudioMixer4              HPswitchL;         //xy=639,421
AudioSynthWaveformSine   sine1;          //xy=645,325
AudioEffectPhaseVocoder  VocoderR;      //xy=855,117
AudioEffectPhaseVocoder  VocoderL;      //xy=858,174
AudioEffectMultiply      multiplyR;      //xy=859,287
AudioEffectMultiply      multiplyL;      //xy=859,363
AudioEffectGranular      granularR;      //xy=504,155
AudioEffectGranular      granularL;      //xy=504,155
AudioMixer4              switchR;         //xy=1042,202
AudioMixer4              switchL;         //xy=1046,362
AudioFilterBiquad        lowpassR;        //xy=1195,230
AudioFilterBiquad        lowpassL;        //xy=1199,331
AudioOutputI2S           i2s2;           //xy=1338,285

// Formant Shifter
AudioConvert_I16toF32    convertInR;  // Convert to float
AudioConvert_I16toF32    convertInL;  // Convert to float
AudioEffectFormantShiftFD_OA_F32 formantShiftR(audio_settings); // the frequency-domain processing block
AudioEffectFormantShiftFD_OA_F32 formantShiftL(audio_settings); // the frequency-domain processing block
AudioEffectGain_F32      formantGainR; //Applies digital gain to audio data.
AudioEffectGain_F32      formantGainL; //Applies digital gain to audio data.
AudioConvert_F32toI16    convertOutR;
AudioConvert_F32toI16    convertOutL;

AudioMixer4_F32          MixerR;  
AudioMixer4_F32          MixerL;

AudioConnection          patchCord31(HPswitchR, 0, convertInR, 0);
AudioConnection          patchCord32(HPswitchL, 0, convertInL, 0);
AudioConnection_F32      patchCord33(convertInR, 0, formantShiftR, 0);
AudioConnection_F32      patchCord34(convertInL, 0, formantShiftL, 0);
AudioConnection_F32      patchCord35(formantShiftR, 0, MixerR, 2);
AudioConnection_F32      patchCord36(formantShiftL, 0, MixerL, 2);

AudioConnection_F32      patchCord39(MixerR, 0, convertOutR, 0);
AudioConnection_F32      patchCord40(MixerL, 0, convertOutL, 0);
AudioConnection          patchCord37(convertOutR, 0, switchR, 2);
AudioConnection          patchCord38(convertOutL, 0, switchL, 2);


// Frequency Mixer
RadioIQMixer_F32         IQMixerR;
RadioIQMixer_F32         IQMixerL;
AudioFilter90Deg_F32     HilbertR;
AudioFilter90Deg_F32     HilbertL;
AudioConnection_F32      patchCord51(convertInR, 0, IQMixerR, 0);
AudioConnection_F32      patchCord52(convertInR, 0, IQMixerR, 1);
AudioConnection_F32      patchCord53(convertInL, 0, IQMixerL, 0);
AudioConnection_F32      patchCord54(convertInL, 0, IQMixerL, 1);
AudioConnection_F32      patchCord55(IQMixerR, 0, HilbertR, 0);
AudioConnection_F32      patchCord56(IQMixerR, 1, HilbertR, 1);
AudioConnection_F32      patchCord57(IQMixerL, 0, HilbertL, 0);
AudioConnection_F32      patchCord58(IQMixerL, 1, HilbertL, 1);
AudioConnection_F32      patchCord59(HilbertR, 0, MixerR, 0);
AudioConnection_F32      patchCord60(HilbertR, 1, MixerR, 1);
AudioConnection_F32      patchCord61(HilbertL, 0, MixerL, 0);
AudioConnection_F32      patchCord62(HilbertL, 1, MixerL, 1);


AudioConnection          patchCord1(i2s1, 0, ampR, 0);
AudioConnection          patchCord2(i2s1, 1, ampL, 0);
AudioConnection          patchCord3(ampL, highpassL);
AudioConnection          patchCord4(ampL, 0, HPswitchL, 0);
AudioConnection          patchCord5(ampR, highpassR);
AudioConnection          patchCord6(ampR, 0, HPswitchR, 0);
AudioConnection          patchCord7(highpassL, 0, HPswitchL, 1);
AudioConnection          patchCord8(highpassR, 0, HPswitchR, 1);
AudioConnection          patchCord9(HPswitchR, VocoderR);
AudioConnection          patchCord10(HPswitchR, 0, multiplyR, 0);
AudioConnection          patchCord11(HPswitchL, VocoderL);
AudioConnection          patchCord12(HPswitchL, 0, multiplyL, 1);
AudioConnection          patchCord13(sine1, 0, multiplyR, 1);
AudioConnection          patchCord14(sine1, 0, multiplyL, 0);
AudioConnection          patchCord15(VocoderR, 0, switchR, 0);
AudioConnection          patchCord16(VocoderL, 0, switchL, 0);
AudioConnection          patchCord17(multiplyR, 0, switchR, 1);
AudioConnection          patchCord18(multiplyL, 0, switchL, 1);
AudioConnection          patchCord19(switchR, lowpassR);
AudioConnection          patchCord20(switchL, lowpassL);

AudioConnection          patchCord70(HPswitchR, granularR);
AudioConnection          patchCord71(HPswitchL, granularL);
AudioConnection          patchCord72(granularR, 0, switchR, 3);
AudioConnection          patchCord73(granularL, 0, switchL, 3);

AudioConnection          patchCord21(lowpassR, 0, i2s2, 0);
AudioConnection          patchCord22(lowpassL, 0, i2s2, 1);
AudioControlSGTL5000     sgtl5000_1;     //xy=154,532
// GUItool: end automatically generated code


elapsedMillis count;

int shift = -12; // 12 semitones == one octave downshift
int heterodyne_freq = 3600;
float formant_scaling = 0.5f; // 1.0f == no shift
int formant_overlap_factor = 4;
float quadrature_freq = 2500.0f;
float delayms = 1.0f;
int highPass_freq = 2500;
int lowPass_freq = 3500;
uint8_t HPfilter_ON = 0; // 0 = OFF, 1 = ON
#define GRANULAR_MEMORY_SIZE 12800  // enough for 290 ms at 44.1 kHz
int16_t granularMemoryR[GRANULAR_MEMORY_SIZE];
int16_t granularMemoryL[GRANULAR_MEMORY_SIZE];
float granular_ratio = 0.5f;
float granular_ms = 30.0f;

#define MODE_MIN          0
#define MODE_HETERODYNE   0
#define MODE_VOCODER      1
#define MODE_FORMANT      2
#define MODE_QUADRATURE   3
#define MODE_GRANULAR     4
#define MODE_MAX          4

int mode = MODE_HETERODYNE;

#define HETERODYNE_STANDARD_FREQ  2500

Bounce button0 = Bounce(26, 20); // UP              Decrease the pitch of the sound
Bounce button1 = Bounce(24, 20); // DOWN            Increase the pitch of the sound
Bounce button2 = Bounce(28, 20); // Klick middle    Restore shift parameter to STANDARD value (= zero shift)  
Bounce button3 = Bounce(32, 20); // RIGHT           Toggle different "modes"= shifting methods
Bounce button4 = Bounce(30, 20); // LEFT            Toggle HighPass Filter ON/OFF

void setup() {
  // put your setup code here, to run once:
  Serial.begin(9600);
  AudioMemory(30);
  AudioMemory_F32(200, audio_settings);
  pinMode(26, INPUT_PULLUP);
  pinMode(24, INPUT_PULLUP);
  pinMode(28, INPUT_PULLUP);
  pinMode(32, INPUT_PULLUP);
  pinMode(30, INPUT_PULLUP);

  sgtl5000_1.enable();
  sgtl5000_1.volume(0.9);

  ampR.gain(2.5f);
  ampL.gain(2.5f);

    // Configure the frequency-domain algorithm for Formant Shift
  //int overlap_factor = 4;  //set to 4 or 8 or either 75% overlap (4x) or 87.5% overlap (8x)
  int N_FFT = audio_block_samples * formant_overlap_factor;
  formantShiftR.setup(audio_settings, N_FFT); //do after AudioMemory_F32();
  formantShiftL.setup(audio_settings, N_FFT); //do after AudioMemory_F32();

  Set_Formant(formant_scaling);

  Set_Vocoder(shift);

  Set_Heterodyne(heterodyne_freq, lowPass_freq, highPass_freq);

  Set_Quadrature(quadrature_freq);
  HilbertR.begin(hilbert251A, 251);  // Set the Hilbert transform FIR filter
  HilbertL.begin(hilbert251A, 251);  // Set the Hilbert transform FIR filter

  Set_HighPass (HPfilter_ON, highPass_freq);
  
  Switch_Mode(mode);

  // the Granular effect requires memory to operate
  granularR.begin(granularMemoryR, GRANULAR_MEMORY_SIZE);
  granularL.begin(granularMemoryL, GRANULAR_MEMORY_SIZE);

  Set_Granular(granular_ratio, granular_ms);
}

void loop() {
  // put your main code here, to run repeatedly:
  button0.update();
  button1.update();
  button2.update();
  button3.update();
  button4.update();
  
  if(count > 2000)
  {
    print_Info();
    count = 0;
  }

    if(mode == MODE_VOCODER)
    {
        if (button0.fallingEdge()) 
        {
          shift += 2;
          Set_Vocoder(shift);
        }
        if (button1.fallingEdge()) 
        {
          shift -= 2;
          Set_Vocoder(shift);
        }
        if (button2.fallingEdge()) 
        {
          shift = 0;
          Set_Vocoder(shift);
        }
    }
    else if(mode == MODE_HETERODYNE)
    {
        if (button0.fallingEdge()) 
        {
          heterodyne_freq -= 200;
          Set_Heterodyne(heterodyne_freq, lowPass_freq, highPass_freq);
        }
        if (button1.fallingEdge()) 
        {
          heterodyne_freq += 200;
          Set_Heterodyne(heterodyne_freq, lowPass_freq, highPass_freq);
        }
        if (button2.fallingEdge()) 
        {
          heterodyne_freq = HETERODYNE_STANDARD_FREQ;
          Set_Heterodyne(heterodyne_freq, lowPass_freq, highPass_freq);
        }
    }
    else if(mode == MODE_FORMANT)
    {
        if (button0.fallingEdge()) 
        {
          formant_scaling -= 0.05f;
          Set_Formant(formant_scaling);
        }
        if (button1.fallingEdge()) 
        {
          formant_scaling += 0.05f;
          Set_Formant(formant_scaling);
        }
        if (button2.fallingEdge()) 
        {
          formant_scaling = 1.0f;
          Set_Formant(formant_scaling);
        }
    }
    else if(mode == MODE_QUADRATURE)
    {
        if (button0.fallingEdge()) 
        {
          quadrature_freq -= 200.0f;
          Set_Quadrature(quadrature_freq);
        }
        if (button1.fallingEdge()) 
        {
          quadrature_freq += 200.0f;
          Set_Quadrature(quadrature_freq);
        }
        if (button2.fallingEdge()) 
        {
          quadrature_freq = 2500.0f;
          Set_Quadrature(quadrature_freq);
        }
    }
    else if(mode == MODE_GRANULAR)
    {
        if (button0.fallingEdge()) 
        {
          granular_ratio -= 0.05f;
          Set_Granular(granular_ratio, granular_ms);
        }
        if (button1.fallingEdge()) 
        {
          granular_ratio += 0.05f;
          Set_Granular(granular_ratio, granular_ms);
        }
        if (button2.fallingEdge()) 
        {
          granular_ratio = 0.5f;
          Set_Granular(granular_ratio, granular_ms);
        }
    }    
    // this is for mode changing    
        if (button3.fallingEdge()) 
        {
          mode += 1;
          mode = clamp_turnover(mode,MODE_MIN,MODE_MAX);
          Switch_Mode(mode);
        }
        if (button4.fallingEdge()) 
        {
          if(HPfilter_ON == 1)
          {
            HPfilter_ON = 0;
          }
          else
          {
            HPfilter_ON = 1;
          }
          Set_HighPass(HPfilter_ON, highPass_freq);
        }
}

void Switch_Mode(int mode)
{
  if(mode == MODE_VOCODER)
  {
    Serial.println("MODE Vocoder was selected");
    ampR.gain(3.5f);
    ampL.gain(3.5f);
    // switch on Vocoder, switch off everything else
    switchR.gain(0,1.0f);
    switchR.gain(1,0.0f);
    switchR.gain(2,0.0f);
    switchR.gain(3,0.0f);
    switchL.gain(0,1.0f);
    switchL.gain(1,0.0f);
    switchL.gain(2,0.0f);
    switchL.gain(3,0.0f);
    // lowpass 
    lowpassR.setLowpass(0,lowPass_freq,0.54);
    lowpassR.setLowpass(1,lowPass_freq,1.3);
    lowpassR.setLowpass(2,lowPass_freq,0.54);
    lowpassR.setLowpass(3,lowPass_freq,1.3);
    lowpassL.setLowpass(0,lowPass_freq,0.54);
    lowpassL.setLowpass(1,lowPass_freq,1.3);
    lowpassL.setLowpass(2,lowPass_freq,0.54);
    lowpassL.setLowpass(3,lowPass_freq,1.3);
        
  }
  else if(mode == MODE_HETERODYNE)
  {
    sine1.amplitude(1.0f);
    sine1.frequency(heterodyne_freq);
    // switch on Heterodyne, switch off everything else
    switchR.gain(0,0.0f);
    switchR.gain(1,1.3f);
    switchR.gain(2,0.0f);
    switchR.gain(3,0.0f);
    switchL.gain(0,0.0f);
    switchL.gain(1,1.3f);
    switchL.gain(2,0.0f);
    switchL.gain(3,0.0f);
    // lowpass
    lowpassR.setLowpass(0,lowPass_freq,0.54);
    lowpassR.setLowpass(1,lowPass_freq,1.3);
    lowpassR.setLowpass(2,lowPass_freq,0.54);
    lowpassR.setLowpass(3,lowPass_freq,1.3);
    lowpassL.setLowpass(0,lowPass_freq,0.54);
    lowpassL.setLowpass(1,lowPass_freq,1.3);
    lowpassL.setLowpass(2,lowPass_freq,0.54);
    lowpassL.setLowpass(3,lowPass_freq,1.3);
    Serial.println("MODE Heterodyne was selected");
  }
  else if(mode == MODE_FORMANT)
  {
    // switch on Formant, switch off everything else
    switchR.gain(0,0.0f);
    switchR.gain(1,0.0f);
    switchR.gain(2,1.0f);
    switchR.gain(3,0.0f);
    switchL.gain(0,0.0f);
    switchL.gain(1,0.0f);
    switchL.gain(2,1.0f);
    switchL.gain(3,0.0f);
    // lowpass
    lowpassR.setLowpass(0,lowPass_freq,0.54);
    lowpassR.setLowpass(1,lowPass_freq,1.3);
    lowpassR.setLowpass(2,lowPass_freq,0.54);
    lowpassR.setLowpass(3,lowPass_freq,1.3);
    lowpassL.setLowpass(0,lowPass_freq,0.54);
    lowpassL.setLowpass(1,lowPass_freq,1.3);
    lowpassL.setLowpass(2,lowPass_freq,0.54);
    lowpassL.setLowpass(3,lowPass_freq,1.3);
    MixerR.gain(0,0.0f);
    MixerR.gain(1,0.0f);
    MixerR.gain(2,1.0f);
    MixerR.gain(3,0.0f);
    MixerL.gain(0,0.0f);
    MixerL.gain(1,0.0f);
    MixerL.gain(2,1.0f);
    MixerL.gain(3,0.0f);
    Serial.println("MODE Formant Shift was selected");  
  }
  else if(mode == MODE_QUADRATURE)
  {
    // switch on Formant or Freq shifting, switch off everything else
    switchR.gain(0,0.0f);
    switchR.gain(1,0.0f);
    switchR.gain(2,1.0f);
    switchR.gain(3,0.0f);
    switchL.gain(0,0.0f);
    switchL.gain(1,0.0f);
    switchL.gain(2,1.0f);
    switchL.gain(3,0.0f);
    // lowpass
    lowpassR.setLowpass(0,lowPass_freq,0.54);
    lowpassR.setLowpass(1,lowPass_freq,1.3);
    lowpassR.setLowpass(2,lowPass_freq,0.54);
    lowpassR.setLowpass(3,lowPass_freq,1.3);
    lowpassL.setLowpass(0,lowPass_freq,0.54);
    lowpassL.setLowpass(1,lowPass_freq,1.3);
    lowpassL.setLowpass(2,lowPass_freq,0.54);
    lowpassL.setLowpass(3,lowPass_freq,1.3);
    MixerR.gain(0,0.5f);
    MixerR.gain(1,0.5f);
    MixerR.gain(2,0.0f);
    MixerR.gain(3,0.0f);
    MixerL.gain(0,0.5f);
    MixerL.gain(1,0.5f);
    MixerL.gain(2,0.0f);
    MixerL.gain(3,0.0f);
    Serial.println("MODE Quadrature was selected");  
  }
  else if(mode == MODE_GRANULAR)
  {
    // switch on Formant or Freq shifting, switch off everything else
    switchR.gain(0,0.0f);
    switchR.gain(1,0.0f);
    switchR.gain(2,0.0f);
    switchR.gain(3,1.0f);
    switchL.gain(0,0.0f);
    switchL.gain(1,0.0f);
    switchL.gain(2,0.0f);
    switchL.gain(3,1.0f);
    // lowpass
    lowpassR.setLowpass(0,lowPass_freq,0.54);
    lowpassR.setLowpass(1,lowPass_freq,1.3);
    lowpassR.setLowpass(2,lowPass_freq,0.54);
    lowpassR.setLowpass(3,lowPass_freq,1.3);
    lowpassL.setLowpass(0,lowPass_freq,0.54);
    lowpassL.setLowpass(1,lowPass_freq,1.3);
    lowpassL.setLowpass(2,lowPass_freq,0.54);
    lowpassL.setLowpass(3,lowPass_freq,1.3);
    Serial.println("MODE Granular was selected");  
  }
  else
  {
    Serial.println("ERROR: no sensible mode selected !");
  }
} // END "Switch_Mode"

void Set_Vocoder(int shift)
{
  VocoderR.setPitchShift(shift);
  VocoderL.setPitchShift(shift);  
  Serial.print(" PITCH SHIFT IN SEMITONES: "); Serial.println(shift);
  Serial.println();
}

void Set_Heterodyne(int freq, int lowpass_freq, int highpass_freq)
{
  heterodyne_freq = clamp(freq,0,10000);
  sine1.frequency(heterodyne_freq);
  Serial.print(" HETERODYNE SHIFT IN HZ: "); Serial.println(heterodyne_freq);
  Serial.println();
}

void Set_Formant(float scale_factor)
{
//  formantShiftR.setScaleFactor(1.587401f); // up-shift !  1.0 is no formant shifting.
//  formantShiftL.setScaleFactor(1.587401f); // up-shift !  1.0 is no formant shifting.
  formant_scaling = clamp_f(scale_factor,0.1f,2.0f);
  formantShiftR.setScaleFactor(formant_scaling); //1.0 is no formant shifting.
  formantShiftL.setScaleFactor(formant_scaling); //1.0 is no formant shifting.
}

void Set_Quadrature(float frequencyLO)
{
  IQMixerR.frequency(frequencyLO);   // Frequency shift, Hz
  IQMixerL.frequency(frequencyLO);   // Frequency shift, Hz
}

void Set_Granular(float ratio, float ms)
{
    granular_ratio = clamp_f(ratio,0.2f, 3.0f);
    granular_ms = clamp_f(ms, 25, 2500);
    granularR.beginPitchShift(granular_ms); // ??? maybe 25 to 2500 ?
    granularR.setSpeed(ratio); // 0.5 to 2.0
    granularL.beginPitchShift(granular_ms); // ??? maybe 25 to 2500 ?
    granularL.setSpeed(ratio); // 0.5 to 2.0
    Serial.print(" GRANULAR RATIO is: "); Serial.print(granular_ratio);
    Serial.println();
}

void Set_HighPass (uint8_t on, int freq)
{
    if(on == 1)
    {
      HPswitchR.gain(0,0.0f);  
      HPswitchR.gain(1,1.0f);
      HPswitchL.gain(0,0.0f);
      HPswitchL.gain(1,1.0f);
    }
    else
    {
      HPswitchR.gain(0,1.0f);  
      HPswitchR.gain(1,0.0f);
      HPswitchL.gain(0,1.0f);
      HPswitchL.gain(1,0.0f);
    }

    HPswitchL.gain(2,0.0f);
    HPswitchL.gain(3,0.0f);
    HPswitchR.gain(2,0.0f);
    HPswitchR.gain(3,0.0f);

    highpassR.setHighpass(0, freq, 0.54);
    highpassR.setHighpass(1, freq, 1.3);
    highpassR.setHighpass(2, freq, 0.54);
    highpassR.setHighpass(3, freq, 1.3);
    highpassL.setHighpass(0, freq, 0.54);
    highpassL.setHighpass(1, freq, 1.3);
    highpassL.setHighpass(2, freq, 0.54);
    highpassL.setHighpass(3, freq, 1.3);  
}

int clamp_turnover(int var, int Min, int Max)
{
  if(var > Max) var = Min;
  if(var < Min) var = Max;
  return var;
}

int clamp(int var, int Min, int Max)
{
  if(var > Max) var = Max;
  if(var < Min) var = Min;
  return var;
}

float clamp_f(float var, float Min, float Max)
{
  if(var > Max) var = Max;
  if(var < Min) var = Min;
  return var;
}

void print_Info(void)
{
    Serial.println();
    Serial.print("We are in ");
    //Serial.println(mode);
    if(mode == MODE_VOCODER) 
    {
      Serial.println(" VOCODOER mode");
      Serial.print(" PITCH SHIFT IN SEMITONES: "); Serial.println(shift);
   }
    else if (mode == MODE_HETERODYNE) 
    {
      Serial.println(" HETERODYNE mode");
      Serial.print(" HETERODYNE FREQUENCY is: "); Serial.print(heterodyne_freq); Serial.println(" Hz"); 
    }
    else if(mode == MODE_QUADRATURE)
    {
      Serial.println("QUADRATURE/FREQUENCY SHIFT mode");
      Serial.print(" FREQUENCY SHIFT is: "); Serial.print(quadrature_freq); Serial.println(" Hz"); 
    }
    else if(mode == MODE_FORMANT)
    {
      Serial.println("FORMANT SHIFT mode");
      Serial.print(" Formant scaling factor is: "); Serial.println(formant_scaling); 
    }    
    else if(mode == MODE_GRANULAR)
    {
      Serial.println("GRANULAR mode");
      Serial.print(" Granular ratio is: "); Serial.println(granular_ratio); 
    }    
    else
    {
      Serial.println(" no sensible mode, something wrong here . . .");
    }
    Serial.print("Audio Input HighPass Filter is ");
    if(HPfilter_ON == 1) Serial.println("ON");
    else Serial.println("OFF"); 
    Serial.print("Audio MEM INT16   Cur/Peak: ");
    Serial.print(AudioMemoryUsage());
    Serial.print("/");
    Serial.println(AudioMemoryUsageMax());
    Serial.print("Audio MEM Float32 Cur/Peak: ");
    Serial.print(AudioMemoryUsage_F32());
    Serial.print("/");
    Serial.println(AudioMemoryUsageMax_F32());    
    
    Serial.print("CPU INT16 Cur/Peak: "); 
    Serial.print (AudioProcessorUsage());
    Serial.print("%/"); 
    Serial.print (AudioProcessorUsageMax());
    Serial.println("%,   ");
    Serial.print("CPU F32   Cur/Peak: ");
    Serial.print(audio_settings.processorUsage());
    Serial.print("%/");
    Serial.print(audio_settings.processorUsageMax());
    Serial.print("%,   ");
    Serial.println();

}

This is the setup (please note that all modules in the bottom right corner are float 32-bit processing modules; everything right and below of the "convertInR")

pitch shifter DD4WH.jpg

If you want to run the sketch, read carefully through the comments in the header: you need to download and copy some files into your local Arduino folder. It runs on Teensy 4 only (and has 75% CPU load).

Now I have the following questions:

* Are there any other methods for pitch shifting that could be useful for bird sounds and that I am missing here ?
* Can you recommend useful settings for the methods above ? I am new to all of those algorithms and they all seem to have their pros and cons.
* Also, it would be great if you could point me to any errors in the implementations of the different methods


EDIT: the current state of the code can also be found here
https://github.com/DD4WH/BirdSongPitchShifter
 
Last edited:
With the help of Lang Elliott and Harold Mills I was able to implement their "simple" time domain overlap-add (OLA) pitch shift algorithm.

It works surprisingly well for pitch shifting bird songs downwards one, one-and-a-half or two octaves (simply change the variable "shift" and recompile). It preserves the harmonics, thus it is a true pitch shift and not a frequency shift method. It does not work well for music, though.

The following sketch should work out-of-the-box without any additional libraries.

To test it out, you can use an audio board OR an audio ADC and DAC OR substitute the I2S input and output by USB objects and play bird song audio from xenocanto into it (e.g. https://xeno-canto.org/explore?query=waldbaumläufer) and listen to the down-pitched songs.

It has about 1.5% MCU load on the Teensy 4.0.

By drawing a sketch of the processing, I try to explain, how this algorithm works: https://user-images.githubuserconte...3110-f01d8397-0838-47c0-8373-3df8eebc1835.png
https://github.com/DD4WH/BirdSongPitchShifter

Code:
/*
 * Pitch shifting in the time domain
 * 
 * implements pitch shifting in the time domain and is based on an idea by Lang Elliott & Herb Susmann for the "Hear birds again"-project,
 * specifically for the now deprecated SongFinder pitch shifter units. 
 * The algorithm can shift the audio down by a factor of two (one octave), three (1.5 octaves) or four (two octaves). 
 * Find a graph showing the implementation for the case of downshifting by 4 here: https://github.com/DD4WH/BirdSongPitchShifter#readme
 * 
 * Many thanks go to Harold Mills & Lang Elliott for explaining this algorithm to me and answering my questions ! :-) 
 * https://hearbirdsagain.org/
 * 
 * (c) Frank DD4WH 2022-10-05
 * 
 * uses Teensy 4.0 and external ADC / DAC connected on perf board with ground plane
 * 
 *  PCM5102A DAC module
    VCC = Vin
    3.3v = NC
    GND = GND
    FLT = GND
    SCL = 23 / MCLK via series 100 Ohm
    BCK = BCLK (21)
    DIN = TX (7)
    LCK = LRCLK (20)
    FMT = GND
    XMT = 3.3V (HIGH)
    
    PCM1808 ADC module:    
    FMT = GND
    MD1 = GND
    MD0 = GND
    GND = GND
    3.3V = 3.3V --> ADC needs both: 5V AND 3V3
    5V = VIN
    BCK = BCLK (21) via series 100 Ohm
    OUT = RX (8)
    LRC = LRCLK (20) via series 100 Ohm
    SCK = MCLK (23) via series 100 Ohm
    GND = GND
    3.3V = 3.3V
 *  
 * MIT license
 */

#include <Audio.h>
#include <Wire.h>
#include <SPI.h>
#include <arm_math.h>
#include <arm_const_structs.h> // in the Teensy 4.0 audio library, the ARM CMSIS DSP lib is already a newer version
#include <utility/imxrt_hw.h>

#define DEBUG
#define PIH 1.5707963267948966192313216916398f
#define FOURPI  (2.0 * TWO_PI)
#define SIXPI   (3.0 * TWO_PI)
int32_t sum;
int idx_t = 0;
float32_t mean;
int16_t *sp_L;
int16_t *sp_R;
double SAMPLE_RATE = 44100;
float32_t audio_gain = 4.0f; 

// this audio comes from the codec by I2S
AudioInputI2S            i2s_in;
AudioRecordQueue         Q_in_L;
AudioRecordQueue         Q_in_R;
AudioMixer4              mixleft;
AudioMixer4              mixright;
AudioPlayQueue           Q_out_L;
AudioPlayQueue           Q_out_R;
AudioOutputI2S           i2s_out;
//AudioOutputUSB           i2s_out;

AudioConnection          patchCord1(i2s_in, 0, Q_in_L, 0);
AudioConnection          patchCord2(i2s_in, 1, Q_in_R, 0);

AudioConnection          patchCord3(Q_out_L, 0, mixleft, 0);
AudioConnection          patchCord4(Q_out_R, 0, mixright, 0);
AudioConnection          patchCord9(mixleft, 0,  i2s_out, 1);
AudioConnection          patchCord10(mixright, 0, i2s_out, 0);

int shift = 4; 
#define BLOCK_SIZE 128
const int N_BLOCKS = 6; //9; // 6 blocks á 128 samples  = 768 samples. No. of samples has to be dividable by 2 AND by 3 AND by 4 // 6 blocks of 128 samples == 768 samples = 17.4ms
const int WINDOW_LENGTH = N_BLOCKS * BLOCK_SIZE;
const int WINDOW_LENGTH_D_2 = WINDOW_LENGTH / 2;
const int IN_BUFFER_SIZE = WINDOW_LENGTH;
#define HOPSIZE_4 WINDOW_LENGTH/4
#define HOPSIZE_3 WINDOW_LENGTH/3
float32_t in_buffer_L[4][IN_BUFFER_SIZE];
float32_t in_buffer_R[4][IN_BUFFER_SIZE];
float32_t out_buffer_L[IN_BUFFER_SIZE];
float32_t out_buffer_R[IN_BUFFER_SIZE];
float32_t add_buffer_L[IN_BUFFER_SIZE / 2];
float32_t add_buffer_R[IN_BUFFER_SIZE / 2];
float32_t window[N_BLOCKS * BLOCK_SIZE];
int buffer_idx = 0;

const float32_t n_att = 90.0; // desired stopband attenuation
const int num_taps = 36; // can be divided by 2, 3 and 4 
// interpolation-by-N
// num_taps has to be a multiple integer of interpolation factor L
// pState is of length (numTaps/L)+blockSize-1 words where blockSize is the number of input samples processed by each call
arm_fir_interpolate_instance_f32 interpolation_R;
float32_t DMAMEM interpolation_R_state [num_taps / 2 + N_BLOCKS * BLOCK_SIZE / 2]; 
arm_fir_interpolate_instance_f32 interpolation_L;
float32_t DMAMEM interpolation_L_state [num_taps / 2 + N_BLOCKS * BLOCK_SIZE / 2]; 
float32_t DMAMEM interpolation_coeffs[num_taps];
int hop0 = 0;
int hop1 = 1;
int hop2 = 2;
int hop3 = 3;

void setup() {
  Serial.begin(115200);
  delay(100);

  AudioMemory(60); // must be high enough to deliver 16 stereo blocks = 32 * 128 samples at any point in time! 
  delay(100);

  /****************************************************************************************
     Audio Shield Setup
  ****************************************************************************************/
  mixleft.gain(0, 1.0);
  mixright.gain(0, 1.0);

  /****************************************************************************************
     begin to queue the audio from the audio library
  ****************************************************************************************/
//  setI2SFreq(SAMPLE_RATE);
  
  Serial.println("DD4WH time domain pitch shifter following Hear Birds Again by Harolds Mills");

  for(unsigned idx=0; idx < WINDOW_LENGTH; idx++)
  { // von Hann window
     window[idx] = 0.5f * (1.0f - cosf(TWO_PI * (float)idx / ((float)(WINDOW_LENGTH - 1))));  
    // Blackman-Nuttall
    //window[idx] = 0.3635819f - 0.4891775f*cosf(2.0*M_PI*(float)idx/((float)(WINDOW_LENGTH-1))) + 0.1365995*cosf(FOURPI*(float)idx/((float)(WINDOW_LENGTH-1))) - 0.0106411f*cosf(SIXPI*(float)idx/((float)(WINDOW_LENGTH-1)));  
  }

  // Interpolation filter
  // the interpolation filter is AFTER the upsampling, so it has to be in the target sample rate!
  calc_FIR_coeffs (interpolation_coeffs, num_taps, (float32_t)4000, n_att, 0, 0.0, SAMPLE_RATE);
  if (arm_fir_interpolate_init_f32(&interpolation_R, (uint8_t)shift, num_taps, interpolation_coeffs, interpolation_R_state, BLOCK_SIZE * N_BLOCKS / (uint32_t)shift)) 
  {
    Serial.println("Init of interpolation failed");
    while(1);
  }
  if (arm_fir_interpolate_init_f32(&interpolation_L, (uint8_t)shift, num_taps , interpolation_coeffs, interpolation_L_state, BLOCK_SIZE * N_BLOCKS / (uint32_t)shift)) 
  {
    Serial.println("Init of interpolation failed");
    while(1);
  }
  delay(100);
  Q_in_L.begin();
  Q_in_R.begin();

} // END OF SETUP


void loop() {
  elapsedMicros usec = 0;
  // are there at least N_BLOCKS buffers in each channel available ?
    if (Q_in_L.available() > N_BLOCKS + 0 && Q_in_R.available() > N_BLOCKS + 0)
    {
      // get audio samples from the audio  buffers and convert them to float
      for (unsigned i = 0; i < N_BLOCKS; i++)
      {
        sp_L = Q_in_L.readBuffer();
        sp_R = Q_in_R.readBuffer();

        // convert to float one buffer_size
        // float_buffer samples are now standardized from > -1.0 to < 1.0
        arm_q15_to_float (sp_L, &in_buffer_L[buffer_idx][BLOCK_SIZE * i], BLOCK_SIZE); // convert int_buffer to float 32bit
        arm_q15_to_float (sp_R, &in_buffer_R[buffer_idx][BLOCK_SIZE * i], BLOCK_SIZE); // convert int_buffer to float 32bit
        Q_in_L.freeBuffer();
        Q_in_R.freeBuffer();
      }
 
      /**********************************************************************************
          Time Domain pitch shift algorithm "OLA" by Harald Mills "Hear birds again" 
       **********************************************************************************/

      /**********************************************************************************
          1 Windowing 
       **********************************************************************************/

//             we apply the second half of the window to the first half of the input buffer
//             works for N==2 
        if(shift == 2)
        {
            for (unsigned i = 0; i < WINDOW_LENGTH_D_2; i++)
            {
              in_buffer_L[0][i] = in_buffer_L[0][i] * window[i + WINDOW_LENGTH_D_2];
              in_buffer_R[0][i] = in_buffer_R[0][i] * window[i + WINDOW_LENGTH_D_2];
            }
//             we apply the first half of the window to the second half of the input buffer
    
            for (unsigned i = 0; i < WINDOW_LENGTH_D_2; i++)
            {
              in_buffer_L[0][i + WINDOW_LENGTH_D_2] = in_buffer_L[0][i + WINDOW_LENGTH_D_2] * window[i];
              in_buffer_R[0][i + WINDOW_LENGTH_D_2] = in_buffer_R[0][i + WINDOW_LENGTH_D_2] * window[i];
            }
        }
        else if(shift == 3)
        {
           for (unsigned i = 0; i < WINDOW_LENGTH; i++)
           {
              in_buffer_L[buffer_idx][i] = in_buffer_L[buffer_idx][i] * window[i];
              in_buffer_R[buffer_idx][i] = in_buffer_R[buffer_idx][i] * window[i];
           }          
        }
        else if(shift == 4)
        {
           for (unsigned i = 0; i < WINDOW_LENGTH; i++)
           {
              in_buffer_L[buffer_idx][i] = in_buffer_L[buffer_idx][i] * window[i];
              in_buffer_R[buffer_idx][i] = in_buffer_R[buffer_idx][i] * window[i];
           }
        } // end shift == 4

    
      /**********************************************************************************
          2 Overlap & Add 
       **********************************************************************************/
        if(shift == 2)
        {
            for (unsigned i = 0; i < WINDOW_LENGTH_D_2; i++)
            {
              add_buffer_L[i] = in_buffer_L[0][i] + in_buffer_L[0][i + WINDOW_LENGTH_D_2];
              add_buffer_R[i] = in_buffer_R[0][i] + in_buffer_R[0][i + WINDOW_LENGTH_D_2];
            }
        }
        else if(shift == 3)
        {   // index of in_buffer    [0][x]  [1][x]  [2][x]  
            if(buffer_idx==2)       {hop0=2, hop1=1, hop2=0;}
            else if(buffer_idx==1)  {hop0=1, hop1=0, hop2=2;}
            else if(buffer_idx==0)  {hop0=0, hop1=2, hop2=1;}
            hop0=hop0*HOPSIZE_3;
            hop1=hop1*HOPSIZE_3;
            hop2=hop2*HOPSIZE_3;
            for (unsigned i = 0; i < HOPSIZE_3; i++)
            {
                add_buffer_L[i] = in_buffer_L[0][i + hop0] + in_buffer_L[1][i + hop1] + in_buffer_L[2][i + hop2]; 
                add_buffer_R[i] = in_buffer_R[0][i + hop0] + in_buffer_R[1][i + hop1] + in_buffer_R[2][i + hop2]; 
            }

            buffer_idx = buffer_idx + 1;  // increment
            if (buffer_idx >=3) {buffer_idx = 0;} // flip-over           
        }
        else if(shift == 4)
        {   // index of in_buffer    [0][x]  [1][x]  [2][x]  [3][x]
            if(buffer_idx == 3)     {hop0=3, hop1=2, hop2=1, hop3=0;}
            else if(buffer_idx==2)  {hop0=2, hop1=1, hop2=0, hop3=3;}
            else if(buffer_idx==1)  {hop0=1, hop1=0, hop2=3, hop3=2;}
            else if(buffer_idx==0)  {hop0=0, hop1=3, hop2=2, hop3=1;}
            hop0=hop0*HOPSIZE_4;
            hop1=hop1*HOPSIZE_4;
            hop2=hop2*HOPSIZE_4;
            hop3=hop3*HOPSIZE_4;
            for (unsigned i = 0; i < HOPSIZE_4; i++)
            {
                add_buffer_L[i] = in_buffer_L[0][i + hop0] + in_buffer_L[1][i + hop1] + in_buffer_L[2][i + hop2] + in_buffer_L[3][i + hop3]; 
                add_buffer_R[i] = in_buffer_R[0][i + hop0] + in_buffer_R[1][i + hop1] + in_buffer_R[2][i + hop2] + in_buffer_R[3][i + hop3]; 
            }

            buffer_idx = buffer_idx + 1;  // increment
            if (buffer_idx >=4) {buffer_idx = 0;} // flip-over  
        }

      /**********************************************************************************
          3 Interpolate 
       **********************************************************************************/

      // receives 512 samples and makes 1024 samples out of it
      // interpolation-by-2
      // interpolation-in-place does not work
      // blocksize is BEFORE zero stuffing
      // recycle in_buffer as temporary buffer
      if(shift == 2)
      {
          arm_fir_interpolate_f32(&interpolation_L, add_buffer_L, in_buffer_L[0], BLOCK_SIZE * N_BLOCKS / (uint32_t)(shift));
          arm_fir_interpolate_f32(&interpolation_R, add_buffer_R, in_buffer_R[0], BLOCK_SIZE * N_BLOCKS / (uint32_t)(shift));
          // do some scaling / gain application after the interpolation
          arm_scale_f32(in_buffer_L[0], audio_gain, out_buffer_L, BLOCK_SIZE * N_BLOCKS);
          arm_scale_f32(in_buffer_R[0], audio_gain, out_buffer_R, BLOCK_SIZE * N_BLOCKS);

      }
      else if(shift == 3)
      {
          arm_fir_interpolate_f32(&interpolation_L, add_buffer_L, in_buffer_L[buffer_idx], BLOCK_SIZE * N_BLOCKS / (uint32_t)(shift));
          arm_fir_interpolate_f32(&interpolation_R, add_buffer_R, in_buffer_R[buffer_idx], BLOCK_SIZE * N_BLOCKS / (uint32_t)(shift));
          // do some scaling / gain application after the interpolation
          arm_scale_f32(in_buffer_L[buffer_idx], audio_gain, out_buffer_L, BLOCK_SIZE * N_BLOCKS);
          arm_scale_f32(in_buffer_R[buffer_idx], audio_gain, out_buffer_R, BLOCK_SIZE * N_BLOCKS);
      }
      else if(shift == 4)
      {
          arm_fir_interpolate_f32(&interpolation_L, add_buffer_L, in_buffer_L[buffer_idx], BLOCK_SIZE * N_BLOCKS / (uint32_t)(shift));
          arm_fir_interpolate_f32(&interpolation_R, add_buffer_R, in_buffer_R[buffer_idx], BLOCK_SIZE * N_BLOCKS / (uint32_t)(shift));
          // do some scaling / gain application after the interpolation
          arm_scale_f32(in_buffer_L[buffer_idx], audio_gain, out_buffer_L, BLOCK_SIZE * N_BLOCKS);
          arm_scale_f32(in_buffer_R[buffer_idx], audio_gain, out_buffer_R, BLOCK_SIZE * N_BLOCKS);
      }

      /**********************************************************************************
           
       **********************************************************************************/

      /**********************************************************************************
          Just copy input into output buffer for testing
       **********************************************************************************/
/*        for (unsigned i = 0; i < WINDOW_LENGTH; i++)
        {
          out_buffer_L[i] = in_buffer_L[i] * audio_gain;
          out_buffer_R[i] = in_buffer_R[i] * audio_gain;
        } */

       /**********************************************************************
          CONVERT TO INTEGER AND PLAY AUDIO - Push audio into I2S audio chain
       **********************************************************************/
      for (unsigned i = 0; i < N_BLOCKS; i++)
        {
          sp_L = Q_out_L.getBuffer();    
          sp_R = Q_out_R.getBuffer();
          arm_float_to_q15 (&out_buffer_L[BLOCK_SIZE * i], sp_L, BLOCK_SIZE); 
          arm_float_to_q15 (&out_buffer_R[BLOCK_SIZE * i], sp_R, BLOCK_SIZE);
          Q_out_L.playBuffer(); // play it !  
          Q_out_R.playBuffer(); // play it !
        }

       /**********************************************************************************
          PRINT ROUTINE FOR ELAPSED MICROSECONDS
       **********************************************************************************/
#ifdef DEBUG
      sum = sum + usec;
      idx_t++;
      if (idx_t > 40) {
        mean = sum / idx_t;
        if (mean / 29.00 / N_BLOCKS * SAMPLE_RATE / AUDIO_SAMPLE_RATE_EXACT < 100.0)
        {
          Serial.print("processor load:  ");
          Serial.print (mean / 29.00 / N_BLOCKS * SAMPLE_RATE / AUDIO_SAMPLE_RATE_EXACT);
          Serial.println("%");
        }
        else
        {
          Serial.println("100%");
        }
        Serial.print (mean);
        Serial.print (" microsec for ");
        Serial.print (N_BLOCKS);
        Serial.print ("  stereo blocks    ");
        //Serial.print("FFT-length = "); Serial.print(FFT_length);
        //Serial.print(";   FIR filter length = "); Serial.println(m_NumTaps);

        idx_t = 0;
        sum = 0;
      }
#endif
    } // end of audio process loop
    
      /**********************************************************************************
          Add button check etc. here
       **********************************************************************************/

} // end loop

void calc_FIR_coeffs (float * coeffs_I, int numCoeffs, float32_t fc, float32_t Astop, int type, float dfc, float Fsamprate)
// pointer to coefficients variable, no. of coefficients to calculate, frequency where it happens, stopband attenuation in dB,
// filter type, half-filter bandwidth (only for bandpass and notch)
{ // modified by WMXZ and DD4WH after
  // Wheatley, M. (2011): CuteSDR Technical Manual. www.metronix.com, pages 118 - 120, FIR with Kaiser-Bessel Window
  // assess required number of coefficients by
  //     numCoeffs = (Astop - 8.0) / (2.285 * TPI * normFtrans);
  // selecting high-pass, numCoeffs is forced to an even number for better frequency response

  float32_t Beta;
  float32_t izb;
  float fcf = fc;
  int nc = numCoeffs;
  fc = fc / Fsamprate;
  dfc = dfc / Fsamprate;
  // calculate Kaiser-Bessel window shape factor beta from stop-band attenuation
  if (Astop < 20.96)
    Beta = 0.0;
  else if (Astop >= 50.0)
    Beta = 0.1102 * (Astop - 8.71);
  else
    Beta = 0.5842 * powf((Astop - 20.96), 0.4) + 0.07886 * (Astop - 20.96);

  for (int i = 0; i < numCoeffs; i++) //zero pad entire coefficient buffer, important for variables from DMAMEM
  {
    coeffs_I[i] = 0.0;
  }

  izb = Izero (Beta);
  if (type == 0) // low pass filter
    //     {  fcf = fc;
  { fcf = fc * 2.0;
    nc =  numCoeffs;
  }
  else if (type == 1) // high-pass filter
  { fcf = -fc;
    nc =  2 * (numCoeffs / 2);
  }
  else if ((type == 2) || (type == 3)) // band-pass filter
  {
    fcf = dfc;
    nc =  2 * (numCoeffs / 2); // maybe not needed
  }
  else if (type == 4) // Hilbert transform
  {
    nc =  2 * (numCoeffs / 2);
    // clear coefficients
    for (int ii = 0; ii < 2 * (nc - 1); ii++) coeffs_I[ii] = 0;
    // set real delay
    coeffs_I[nc] = 1;

    // set imaginary Hilbert coefficients
    for (int ii = 1; ii < (nc + 1); ii += 2)
    {
      if (2 * ii == nc) continue;
      float x = (float)(2 * ii - nc) / (float)nc;
      float w = Izero(Beta * sqrtf(1.0f - x * x)) / izb; // Kaiser window
      coeffs_I[2 * ii + 1] = 1.0f / (PIH * (float)(ii - nc / 2)) * w ;
    }
    return;
  }

  for (int ii = - nc, jj = 0; ii < nc; ii += 2, jj++)
  {
    float x = (float)ii / (float)nc;
    float w = Izero(Beta * sqrtf(1.0f - x * x)) / izb; // Kaiser window
    coeffs_I[jj] = fcf * m_sinc(ii, fcf) * w;

  }

  if (type == 1)
  {
    coeffs_I[nc / 2] += 1;
  }
  else if (type == 2)
  {
    for (int jj = 0; jj < nc + 1; jj++) coeffs_I[jj] *= 2.0f * cosf(PIH * (2 * jj - nc) * fc);
  }
  else if (type == 3)
  {
    for (int jj = 0; jj < nc + 1; jj++) coeffs_I[jj] *= -2.0f * cosf(PIH * (2 * jj - nc) * fc);
    coeffs_I[nc / 2] += 1;
  }

} // END calc_FIR_coeffs

float32_t Izero (float32_t x)
{
  float32_t x2 = x / 2.0;
  float32_t summe = 1.0;
  float32_t ds = 1.0;
  float32_t di = 1.0;
  float32_t errorlimit = 1e-9;
  float32_t tmp;
  do
  {
    tmp = x2 / di;
    tmp *= tmp;
    ds *= tmp;
    summe += ds;
    di += 1.0;
  }   while (ds >= errorlimit * summe);
  return (summe);
}  // END Izero

float m_sinc(int m, float fc)
{ // fc is f_cut/(Fsamp/2)
  // m is between -M and M step 2
  //
  float x = m * PIH;
  if (m == 0)
    return 1.0f;
  else
    return sinf(x * fc) / (fc * x);
}
 
This looks really interesting. I have some desires to make a vocoder for voice input, and also pitch shifting sounds of bats and other animals that have content well above even good human hearing. Thanks for posting! I'm forking your GitHub now.
 
Back
Top