Oversampling

[Stefan Westerfeld]

   Hi!

I've implmeneted some code for per-module oversampling. It is based on
FIR filters which are designed from a windowed sinc function.
Oversampling ratios from 2 to 16 are supported.

To test it, I've implemented Bse::Rectify, a plugin which should -
without oversampling - generate extreme aliasing.

I am attaching the code here, and an example .bse file. When running the
bse file, use a fft scope to monitor the rectify output, while playing
with the oversample ratio setting.

Comments are welcome.

   Cu... Stefan
--
Stefan Westerfeld, Hamburg/Germany, http://space.twc.de/~stefan

/* BseNoise - BSE FM Operator -*-mode: c++;-*-
 * Copyright (C) 2006 Stefan Westerfeld
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General
 * Public License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place, Suite 330,
 * Boston, MA 02111-1307, USA.
 */
#include <bse/bsecxxmodule.idl>

namespace Bse {

class Rectify : Effect {
  Info    authors     = "Stefan Westerfeld";
  Info    license     = _("GNU Lesser General Public License");
  Info    category    = _("/Distortion/Rectify");
  Info    blurb       = _("Outputs -1 for signals < 0 and +1 for signals > 0");
  IStream audio_in    = (_("Audio In"), _("Audio Input"));
  OStream audio_out   = (_("Audio Out"), _("Audio Output"));

  group _("General") {
    Int   oversample  = (_("Oversampling"), _("The amount of internal oversampling. If BSE is running with a rate of 48000 Hz, "
                                              "an oversampling of 2 means to perform internal computations with 2*48000 Hz = 96000 Hz. "
                                              "\n\n"
                                              "Higher oversampling means better quality (less aliasing), but also more CPU usage."),
                         1, 1, 16, 1, STANDARD) ;
    Bool  bypass      = Bool (_("Bypass Rectification"), _("Bypasses Rectification (but still does oversampling"), FALSE, STANDARD);

  };
};
};

---

/* BseRectify - generates rectangle signals out of normal ones
 * Copyright (C) 2006 Stefan Westerfeld
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General
 * Public License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place, Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#include "bserectify.genidl.hh"
#include <bse/gslfft.h>
#include <bse/bsemathsignal.h>
#include <vector>
#include <complex> /* debugging only */

using namespace std;

namespace Bse {
namespace {

static double
sinc (double x)
{
  return sin (x * M_PI) / (x * M_PI);
}

static void
design_filter (int oversample_ratio, int taps_per_sample, double *filter /* length: taps_per_sample * (oversample_ratio - 1) */)
{
  g_return_if_fail ((taps_per_sample & 1) == 0);

  int o = 0;
  for (int offset = 1; offset < oversample_ratio; offset++)
    {
      gdouble window_width = taps_per_sample / 2.0 * oversample_ratio;
      for (int i = 0; i < taps_per_sample; i++)
        {
          gdouble pos = oversample_ratio * i - window_width + offset;
          filter[o++] = bse_window_blackman (pos / window_width) * sinc (pos / oversample_ratio);
          //printf ("filter[%d] = %f\n", o-1, filter[o-1]);
        }
    }

  g_assert (o == taps_per_sample * (oversample_ratio - 1));

#if 0 // debugging
  gdouble sinc_filter[oversample_ratio*taps_per_sample];
  for (int i = 0; i < taps_per_sample; i++)
    {
      for (int j = 0; j < oversample_ratio-1; j++)
        sinc_filter[i * oversample_ratio+j] = filter[i+j*taps_per_sample];
      sinc_filter[i * oversample_ratio + oversample_ratio - 1] = 0.0;
    }
  sinc_filter[taps_per_sample / 2 * oversample_ratio - 1] = 1.0;
  gfloat freq;
  for (freq = 0; freq < PI; freq += 0.001)
    {
      std::complex<double> z = exp (complex<double> (0, freq));
      std::complex<double> r = 0;

      guint i;
      for (i = 0; i < oversample_ratio*taps_per_sample; i++)
        {
          r /= z;
          r += sinc_filter[i];
          printf ("sinc_filter[%d] = %f\n", i, sinc_filter[i]);
        }
      printf ("%f %f FTRANS\n", freq, abs (r));
    }
  gdouble linear_filter[3] = { 0.5, 1, 0.5 };
  for (freq = 0; freq < PI; freq += 0.001)
    {
      std::complex<double> z = exp (complex<double> (0, freq));
      std::complex<double> r = 0;

      guint i;
      for (i = 0; i < 3; i++)
        {
          r /= z;
          r += linear_filter[i];
        }
      printf ("%f %f FTLIN\n", freq, abs (r));
    }
  gdouble downsampler[3] = { 0.5, 0.5 };
  gfloat freq;
  for (freq = 0; freq < PI; freq += 0.001)
    {
      std::complex<double> z = exp (complex<double> (0, freq));
      std::complex<double> r = 0;

      guint i;
      for (i = 0; i < 2; i++)
        {
          r /= z;
          r += downsampler[i];
        }
      printf ("%f %f FTDOWN\n", freq, abs (r));
    }
#endif
}

class OversampleUp {
  const IStream& istream;
  guint ratio;
  static const int HMASK = 511; /* must be > 16 * TAPS_PER_SAMPLE, must be bitmask ==> 2^k-1 */
  static const int TAPS_PER_SAMPLE = 20;
  vector<gfloat> history;
  vector<double> filter;
  guint history_write_pos;
public:
  OversampleUp (const IStream& istream, guint ratio)
    : istream (istream),
      ratio (ratio),
      history (HMASK + 1),
      history_write_pos (0),
      filter (TAPS_PER_SAMPLE * (ratio - 1))
  {
    g_return_if_fail (ratio >= 1 && ratio <= 16);
    design_filter (ratio, TAPS_PER_SAMPLE, &filter[0]);
  }
  gfloat hist (int index) const
  {
    return history[(history_write_pos + index) & HMASK];
  }
  const gfloat *upsample (gfloat *buffer, unsigned int n_buffer_values)
  {
    if (ratio == 1)
      return istream.values;

    if (ratio == 2)
      {
        const gfloat *input = istream.values;
        for (int i = 0, j = 0; i < n_buffer_values; i++)
          {
            if (i % 2) /* apply filter */
              {
                //buffer[i] = (hist(-2) + hist(-1)) * 0.5; // linear
                buffer[i] = 0;
                for (int tap = 0; tap < TAPS_PER_SAMPLE; tap++)
                  buffer[i] += hist(tap-TAPS_PER_SAMPLE) * filter[tap];
              }
            else /* copy data */
              {
                buffer[i] = hist(-TAPS_PER_SAMPLE/2);
                //buffer[i] = hist(-1);
                history[history_write_pos++ & HMASK] = input[j++];
              }
          }
      }
    else
      {
        const gfloat *input = istream.values;
        for (int i = 0, j = 0; i < n_buffer_values; i++)
          {
            if (i % ratio) /* apply filter */
              {
                /*
                 * FIXME: why is fstart computation soo difficult?
                 * I thought I had it tweaked to the point where filter could be
                 * simply read from start to end, instead of "quasi-reverse" adressing...?
                 */
                int fstart = (ratio - (i % ratio) - 1) * TAPS_PER_SAMPLE;
                buffer[i] = 0;
                for (int tap = 0; tap < TAPS_PER_SAMPLE; tap++)
                  buffer[i] += hist(tap-TAPS_PER_SAMPLE) * filter[tap + fstart];
              }
            else /* copy data */
              {
                buffer[i] = hist(-TAPS_PER_SAMPLE/2);
                history[history_write_pos++ & HMASK] = input[j++];
              }
          }
      }
    return buffer;
  }
};

class OversampleDown {
  const OStream& ostream;
  guint ratio;
  gfloat *buffer;
  unsigned int n_buffer_values;
  static const int HMASK = 511; /* must be > 16 * TAPS_PER_SAMPLE, must be bitmask ==> 2^k-1 */
  static const int TAPS_PER_SAMPLE = 20;
  vector<gfloat> history;
  vector<double> filter;
  guint history_write_pos;
public:
  OversampleDown (const OStream& ostream, guint ratio)
    : ostream (ostream),
      ratio (ratio),
      history (HMASK + 1),
      history_write_pos (0),
      filter (TAPS_PER_SAMPLE * (ratio - 1))
  {
    g_return_if_fail (ratio >= 1 && ratio <= 16);
    design_filter (ratio, TAPS_PER_SAMPLE, &filter[0]);
    for (guint i = 0; i < filter.size(); i++)
      filter[i] /= ratio;
  }
  gfloat hist (int index) const
  {
    return history[(history_write_pos + index) & HMASK];
  }
  gfloat *downsample_begin (gfloat *buffer, unsigned int n_buffer_values)
  {
    if (ratio == 1)
      return ostream.values;

    this->buffer = buffer;
    this->n_buffer_values = n_buffer_values;

    return buffer;
  }
  void downsample_end()
  {
    if (ratio == 1)
      return;
    if (ratio == 2)
      {
        gfloat *output = ostream.values;
        for (int i = 0; i < n_buffer_values; i += 2)
          {
            history[history_write_pos++ & HMASK] = buffer[i];
            history[history_write_pos++ & HMASK] = buffer[i+1];

            double out = 0;
            int pos = -2 * TAPS_PER_SAMPLE + 1;
            for (int tap = 0; tap < TAPS_PER_SAMPLE; tap++, pos += 2)
              out += hist (pos) * filter[tap];
            out += hist (-TAPS_PER_SAMPLE) * 0.5;

            //out = hist (-2) + hist (-1); // linear
            output[i / 2] = out;
          }
      }
    else
      {
      gfloat *output = ostream.values;
        for (int i = 0; i < n_buffer_values; i += ratio)
          {
            for (int j = 0; j < ratio; j++)
              history[history_write_pos++ & HMASK] = buffer[i+j];

            double out = 0;
            int pos = -ratio * TAPS_PER_SAMPLE + 1; // slow! swap loops!
            for (int tap = 0; tap < TAPS_PER_SAMPLE; tap++, pos += ratio)
              {
                for (int j = 0; j < (ratio-1); j++)
                  {
                    /*
                     * FIXME: if fstart computation is "the wrong direction" above,
                     * why is it easy here? maybe we don't do the right thing?
                     */
                    int fstart = j * TAPS_PER_SAMPLE;
                    out += hist (pos+j) * filter[tap+fstart];
                  }
              }
            out += hist (-TAPS_PER_SAMPLE / 2 * ratio) / ratio;

            // out = (hist (-3) + hist (-2) + hist (-1)) / 3.0; // linear
            output[i / ratio] = out;
          }
}
  }
};

}; /* anon namespace */

class Rectify : public RectifyBase {
  /* actual computation */
  class Module : public SynthesisModule {
    int              oversample;
    bool             bypass;
    OversampleUp    *oversample_audio_in;
    OversampleDown  *oversample_audio_out;
  public:
    void
    config (RectifyProperties *properties)
    {
      set_oversample (properties->oversample);
      bypass = properties->bypass;
    }
    void
    reset()
    {
    }
    Module()
    {
      oversample_audio_in = 0;
      oversample_audio_out = 0;
    }
    ~Module()
    {
      set_oversample (0); /* delete oversample modules */
    }
    void
    set_oversample (int ratio)
    {
      oversample = ratio;
      if (oversample_audio_in)
        delete oversample_audio_in;
      if (oversample_audio_out);
        delete oversample_audio_out;
      if (ratio)
        {
          /* FIXME: malloc in RT thread */
          oversample_audio_in = new OversampleUp (istream (ICHANNEL_AUDIO_IN), ratio);
          oversample_audio_out = new OversampleDown (ostream (OCHANNEL_AUDIO_OUT), ratio);
        }
    }

    void
    process (unsigned int n_values)
    {
      if (!istream (ICHANNEL_AUDIO_IN).connected)
        {
          ostream_set (OCHANNEL_AUDIO_OUT, const_values (0));
          return;
        }

      n_values *= oversample;

      gfloat audio_in_buffer[n_values], audio_out_buffer[n_values];
      const gfloat *audio_in = oversample_audio_in->upsample (audio_in_buffer, n_values);
      gfloat *audio_out = oversample_audio_out->downsample_begin (audio_out_buffer, n_values);

      if (bypass)
        {
          copy (audio_in, audio_in + n_values, audio_out);
        }
      else
        {
          for (unsigned int i = 0; i < n_values; i++)
            audio_out[i] = (audio_in[i] < 0) ? -1 : 1;
        }

      oversample_audio_out->downsample_end();
    }
  };
public:
  /* implement creation and config methods for synthesis Module */
  BSE_EFFECT_INTEGRATE_MODULE (Rectify, Module, RectifyProperties);
};

BSE_CXX_DEFINE_EXPORTS();
BSE_CXX_REGISTER_EFFECT (Rectify);

} // Bse

#ifdef TEST_STANDALONE
int
main()
{
  // test upsampling
    {
      const int R = 7;
      const int N = 1024;
      float in[N];
      for (int i = 0; i < N; i++)
        in[i] = sin (i * 2 * M_PI * 440.0 / 44100.0);

      Bse::IStream is;
      is.values = in;
      float out[R * N];
      Bse::OversampleUp os (is, R);
      os.upsample (out, R * N);
      for (int i = 0; i < R * N; i++)
        printf ("%d %f UP\n", i, out[i]);
    }

  // test downsampling
    {
      float high_buffer[2048];
      float out[1024];
      Bse::OStream os;
      os.values = out;
      Bse::OversampleDown ds (os, 2);

      float *high = ds.downsample_begin (high_buffer, 2048);
      for (int i = 0; i < 2048; i++)
        high[i] = sin (i * 2 * M_PI * 440.0 / 88200.0) // desired signal
                + sin (i * 2 * M_PI * 30000.0 / 88200.0)*0.1; // undesired signal

      ds.downsample_end();

      for (int i = 0; i < 1024; i++)
        printf ("%d %f DOWN\n", i, out[i]);
    }
  return 0;
}
#endif

---

; BseProject

(bse-version "0.6.6")

(container-child "BseWaveRepo::Wave-Repository"
  (modification-time "2006-02-17 07:26:13")
  (creation-time "2006-02-17 07:26:13"))
(container-child "BseCSynth::Unnamed"
  (auto-activate #t)
  (modification-time "2006-02-17 07:26:18")
  (creation-time "2006-02-17 07:26:18")
  (container-child "BseStandardOsc::StandardOsc-1"
    (pulse-mod-perc 0)
    (pulse-width 50)
    (fm-n-octaves 1)
    (exponential-fm #f)
    (fm-perc 0)
    (base-freq 440)
    (wave-form bse-standard-osc-saw-fall)
    (pos-y 1)
    (pos-x -3))
  (container-child "BseAmplifier::Amplifier-1"
    (master-volume 1)
    (base-level 2.7000000000000011)
    (ostrength 100)
    (ctrl-exp #f)
    (ctrl-mul #t)
    (clevel2 100)
    (clevel1 100)
    (alevel2 100)
    (alevel1 100)
    (pos-y 1)
    (pos-x 1)
    (source-input "audio-in1" (link 1 "Rectify-1") "audio-out"))
  (container-child "BseRectify::Rectify-1"
    (oversample 1)
    (pos-y 1)
    (pos-x -1)
    (source-input "audio-in" (link 1 "StandardOsc-1") "audio-out"))
  (container-child "BsePcmOutput::PcmOutput-1"
    (pos-y 1.03)
    (pos-x 3.15)
    (source-input "left-audio-in" (link 1 "Amplifier-1") "audio-out")
    (source-input "right-audio-in" (link 1 "Amplifier-1") "audio-out")))

_______________________________________________
beast mailing list
beast@...
http://mail.gnome.org/mailman/listinfo/beast

[Tim Janik]

On Tue, 21 Feb 2006, Stefan Westerfeld wrote:

>   Hi!
>
> I've implmeneted some code for per-module oversampling. It is based on
> FIR filters which are designed from a windowed sinc function.
> Oversampling ratios from 2 to 16 are supported.
>
> To test it, I've implemented Bse::Rectify, a plugin which should -
> without oversampling - generate extreme aliasing.
>
> I am attaching the code here, and an example .bse file. When running the
> bse file, use a fft scope to monitor the rectify output, while playing
> with the oversample ratio setting.

hm, thinking about the general test casing you developed here,
could you turn this into an audio unit test based on bsefeature
extract/compare?

>
> Comments are welcome.
>
>   Cu... Stefan
> --
> Stefan Westerfeld, Hamburg/Germany, http://space.twc.de/~stefan
>

---
ciaoTJ
_______________________________________________
beast mailing list
beast@...
http://mail.gnome.org/mailman/listinfo/beast

[Tim Janik]

On Tue, 2 May 2006, Stefan Westerfeld wrote:

>   Hi!
>
> On Tue, May 02, 2006 at 12:11:02AM +0200, Tim Janik wrote:
>> On Tue, 21 Feb 2006, Stefan Westerfeld wrote:
>>> I've implmeneted some code for per-module oversampling. It is based on
>>> FIR filters which are designed from a windowed sinc function.
>>> Oversampling ratios from 2 to 16 are supported.
>>>
>>> To test it, I've implemented Bse::Rectify, a plugin which should -
>>> without oversampling - generate extreme aliasing.
>>>
>>> I am attaching the code here, and an example .bse file. When running the
>>> bse file, use a fft scope to monitor the rectify output, while playing
>>> with the oversample ratio setting.
>>
>> hm, thinking about the general test casing you developed here,
>> could you turn this into an audio unit test based on bsefeature
>> extract/compare?
>
> I don't exactly know what you mean by that. If you mean that bsefextract
> bsefcompare should somehow take the role of a human being looking at the
> spectrum and saying "well, that looks like it doesn't have aliasing", I
> don't think thats feasible. That relies too much on listening and
> experience to be automated.
>
> But of course what could be done is to oversample a plugin properly (say
> 16 times), and generate a average spectrum from that. After you listened
> to it and convinced yourself it sounds "right", you could then compare
> it to the average spectrum of a non-oversampled run automatically, and
> see how this differs.
>
> Hoever, this method may still fail, because the non-oversampled version
> may not produce output in "inaudible" frequencies (>18kHz), or less
> output in high frequencies in general, so they may still differ.
>
> So I think in the end, when it comes to aliasing you always have to
> check manually in some way or other.

i was really thinking of something simpler:

> Of your bsefextract/compare can be used to detect regressions after you
> once have validated that some version is correct.

yeah, that's what i was thinking of. i.e. a simple test that figures:
a) the signal is really resampled, i.e. != the original
b) the resampled signal matches the features extracted by a previous run.

so we'd figure if some change breaks the resampling logic in principle
(e.g. the move to SSE or similar).

>   Cu... Stefan

---
ciaoTJ
_______________________________________________
beast mailing list
beast@...
http://mail.gnome.org/mailman/listinfo/beast

[Stefan Westerfeld]

   Hi!

On Tue, May 02, 2006 at 12:11:02AM +0200, Tim Janik wrote:

> On Tue, 21 Feb 2006, Stefan Westerfeld wrote:
> >I've implmeneted some code for per-module oversampling. It is based on
> >FIR filters which are designed from a windowed sinc function.
> >Oversampling ratios from 2 to 16 are supported.
> >
> >To test it, I've implemented Bse::Rectify, a plugin which should -
> >without oversampling - generate extreme aliasing.
> >
> >I am attaching the code here, and an example .bse file. When running the
> >bse file, use a fft scope to monitor the rectify output, while playing
> >with the oversample ratio setting.
>
> hm, thinking about the general test casing you developed here,
> could you turn this into an audio unit test based on bsefeature
> extract/compare?

I don't exactly know what you mean by that. If you mean that bsefextract
bsefcompare should somehow take the role of a human being looking at the
spectrum and saying "well, that looks like it doesn't have aliasing", I
don't think thats feasible. That relies too much on listening and
experience to be automated.

But of course what could be done is to oversample a plugin properly (say
16 times), and generate a average spectrum from that. After you listened
to it and convinced yourself it sounds "right", you could then compare
it to the average spectrum of a non-oversampled run automatically, and
see how this differs.

Hoever, this method may still fail, because the non-oversampled version
may not produce output in "inaudible" frequencies (>18kHz), or less
output in high frequencies in general, so they may still differ.

So I think in the end, when it comes to aliasing you always have to
check manually in some way or other.

Of your bsefextract/compare can be used to detect regressions after you
once have validated that some version is correct.

   Cu... Stefan
--
Stefan Westerfeld, Hamburg/Germany, http://space.twc.de/~stefan
_______________________________________________
beast mailing list
beast@...
http://mail.gnome.org/mailman/listinfo/beast