Single sideband amplitude modulation shifts the frequency either upward or downward, but not both directions (as in AM). This technique is based on the trigonometric identity we used earlier. Shifting a cosine up simply takes half of a complex multiplication.

cos(\alpha+\beta)=cos(\alpha)cos(\beta)−sin(\alpha)sin(\beta)

Alpha and beta represent the frequencies of the signal sources. Each of these frequencies requires a cosine and sine oscillator pair. By subtracting the sin() product from the cos() product, a frequency shifted (alpha+beta) oscillator is synthesized. If the products are summed, a downward (alpha-beta) oscillator is synthesized. As cosine and sine are separated by 90 degrees, it is simple enough to create a cosine/sine pair from a sinusoidal oscillator – simply add 90 degrees to the phase of the oscillator to create the second of the pair.

However, single sideband AM is more often used to frequency shift sampled signals. To do this, it is useful to consider the sample as a sum of multiple sinusoidal oscillators. To perform SSB AM, the original sample as well as a 90 degree phase shifted sample is needed. Phase shifting all of the sinusoidal components of a sample 90 degrees can be done with a type of allpass filter called a Hilbert transformer. If well designed, a Hilbert transformer will shift most frequencies 90 degrees, and will maintain this phase shift until close to 0Hz and the Nyquist frequency. Many implementations have two outputs which are 90 degrees apart, but not in phase with the original sample.

Calculating coefficients for a Hilbert transformer is beyond the scope of this class, but one could read articles by Scott Wardle or Dan Harris, Edgar Berdahl and Jonathan Abel for more information in this area.

Implementation

There are several open source Hilbert transformers in popular audio software. One is in Miller Puckette’s Pure Data, and is implemented as 2 pairs of biquad filters (a biquad filter was used in the earlier second order allpass filter example). It is based on a 4X synthesizer patch by Emmanuel Favreau. The other is in the Csound source code, designed by Sean Costello and based on data from Bernie Hutchins.

The difference equations for the Favreau Hilbert transformer are:

\scriptsize x1[n]=x3[n]=in\\x2[n]=y1[n]=0.94657x1[n]-1.94632x1[n-1]+x1[n-2]+1.94632y1[n-1]-0.94657y1[n-2]\\ outsin = y2[n] = 0.06338x2[n]-0.83774x2[n-1]+x2[n-2]+0.83774y2[n-1]-0.06338y2[n-2]\\x4[n]=y3[n]=-0.260502x3[n]+0.02569x3[n-1]+x3[n-2]-0.02569y3[n-1]+0.260502y3[n-2]\\ outcos = y4[n] = 0.870686x4[n]-1.8685x4[n-1]+x4[n-2]+01.8685y4[n-1]-0.870686y4[n-2]

Notice that this specifies two pairs of allpass filters in which each pair is in series. The C code follows directly from these equations.

// the samplerate in Hz
float sampleRate = 44100;
// initial phase
float phase = 0.0;
long tableSize = 8192;
long tableMask = 8191;
float sineTab[8192];
float x1[3], x2[3], x3[3], x4[3];
float y1[3], y2[3], y3[3], y4[3];

void InitTables(void)
{
  // create sineTab from sin()
  for(i = 0; i < tableSize; i++)
  {
    sineTab[i] = sin((6.283185307179586 * i)/tableSize);
   }
}


// 4 biquad filters
void Hilbert(float input, float *sinout, float *cosout)
{
  x1[0] = x3[0] = input;
  x2[0] = y1[0] = 0.94657*x1[0] − 1.94632*x1[1] + x1[2]
    + 1.94632*y1[1] − 0.94657*y1[2];
  y1[2] = y1[1];
  y1[1] = y1[0];
  x1[2] = x1[1];
  x1[1] = x1[0];
  *sinout = y2[0] = 0.06338*x2[0] − 0.83774x2[1] + x2[2] 
    + 0.83774*y2[1] − 0.06338*y2[2];
  y2[2] = y2[1];
  y2[1] = y2[0];
  x2[2] = x2[1];
  x2[1] = x2[0];

  x4[0] = y3[0] = -0.260502*x3[0] + 0.02569*x3[1] + x3[2] 
    - 0.02569*y3[1] + 0.260502*y3[2];
  y3[2] = y3[1];
  y3[1] = y3[0];
  x3[2] = x3[1];
  x3[1] = x3[0];
  *cosout = y4[0] = 0.870686*x4[0] − 1.8685*x4[1] + x4[2] 
    + 1.8685*y4[1] − 0.870686y4[2];
  y4[2] = y4[1];
  y4[1] = y4[0];
  x4[2] = x4[1];
  x4[1] = x4[0];
}
void FreqShift(float *input, float *output, float frequency, float ratio, float ringGain, float sineGain, float squareGain, long samples)
{
    long sample;
    float phaseInc;
    float sinOsc, cosOsc, sinSamp, cosSamp;
    // calculate for each sample in a block
    for(sample = 0; sample < samples; sample++)
    {
        // get the phase increment for this sample
        phaseInc = frequency/sampleRate;

        // calculate the waveforms
        sinOsc = sineTab[(long)(phase * tableSize)];
        cosOsc = sineTab[((long)(phase+0.25 * tableSize))%tableMask];
        Hilbert(*output+sample, &sinSamp, &cosSamp);

        *(output+sample) = sinSamp * sinOsc - cosSamp * cosOsc; 

        // increment the phases
        phase = phase + phaseInc;
        if(phase >= 1.0)
            phase = phase - 1.0;
  }
}