added resampling

Closes #2

go.mod

@@ -2,5 +2,5 @@ module zikichombo.org/dsp
 
 require (
 	zikichombo.org/codec v0.0.1-alpha.3 // indirect
-	zikichombo.org/sound v0.1.1-alpha.2
+	zikichombo.org/sound v0.1.2-alpha.2
 )

go.sum

@@ -2,3 +2,5 @@ zikichombo.org/codec v0.0.1-alpha.3 h1:g09S6DOXyO56Q3rpwHL813dwEDeNIUoFhxCOk/8sk
 zikichombo.org/codec v0.0.1-alpha.3/go.mod h1:JQb94Z/twHniAfp0ydPKk8ibq61UoBMriWo7iQopcR8=
 zikichombo.org/sound v0.1.1-alpha.2 h1:x2bLgomOt5nKdX5feIkRn1+cthP8YLXLqa62cBmFfVs=
 zikichombo.org/sound v0.1.1-alpha.2/go.mod h1:akiZR7uLjsosCVBmP0t0hqsOgZ7iydYzjz92wFK2pLQ=
+zikichombo.org/sound v0.1.2-alpha.2 h1:x6hKOgVBKH85M7JlOmjNEJ1mPogRYiLdkWAN4QqX3OU=
+zikichombo.org/sound v0.1.2-alpha.2/go.mod h1:akiZR7uLjsosCVBmP0t0hqsOgZ7iydYzjz92wFK2pLQ=

resample/ct.go

@@ -2,16 +2,21 @@ package resample
 
 import (
 	"errors"
+	"io"
 	"math"
 
 	"zikichombo.org/dsp/wfn"
 	"zikichombo.org/sound"
+	"zikichombo.org/sound/cil"
+	"zikichombo.org/sound/freq"
 )
 
 const (
 	shiftSize = 64
 )
 
+// Type C holds state for giving a continuous time
+// representation of a sound.Source.
 type C struct {
 	src     sound.Source
 	shift   int
@@ -24,6 +29,123 @@ type C struct {
 	itper   Itper
 }
 
+// SampleRateConverter provides an interface to a dynamic resample rate
+// conversion.
+type SampleRateConverter interface {
+	// Convert is called by the resampling methods in this package to determine
+	// the result frequency in a conversion.
+	//
+	// It is called once for every output sample except the first sample, which
+	// is taken to be at the same point in time as the first input sample.
+	//
+	// The return value should provide the ratio of the output rate to the input
+	// rate.  It is assumed the input rate is fixed and determined by calling
+	// context.
+	Convert() float64
+}
+
+// ConstSampleRateStretcher is a type which implements SampleRateConverter
+// based on a float64 constant sample rate conversion ratio.
+type constSampleRateConverter float64
+
+// Stretch implements SampleRateStretcher
+func (c constSampleRateConverter) Convert() float64 {
+	return float64(c)
+}
+
+// DynResampler is used to dynamically resample a source.
+// It does not implement sound.Source, since the sample rate is
+// fixed.
+type DynResampler struct {
+	ct    *C
+	src   sound.Source
+	conv  SampleRateConverter
+	lasti float64
+	buf   []float64
+}
+
+// NewDynResampler creates a new Dynamic Resampler from a continuous
+// time representation and a sample rate converter.
+func NewDynResampler(c *C, conv SampleRateConverter) *DynResampler {
+	return &DynResampler{ct: c, conv: conv, lasti: 0.0, buf: make([]float64, c.Channels())}
+}
+
+// DynResampler returns the number of channels.
+func (r *DynResampler) Channels() int {
+	return r.ct.src.Channels()
+}
+
+// Close implements sound.Close
+func (r *DynResampler) Close() error {
+	return r.ct.Close()
+}
+
+// Receive is as in sound.Source.Receive.
+func (r *DynResampler) Receive(d []float64) (int, error) {
+	nC := r.ct.Channels()
+	if len(d)%nC != 0 {
+		return 0, sound.ErrChannelAlignment
+	}
+	nF := len(d) / nC
+	for f := 0; f < nF; f++ {
+		if err := r.ct.FrameAt(r.buf, r.lasti); err != nil {
+			if err == io.EOF {
+				cil.Compact(d, nC, f)
+				return f, nil
+			}
+			return 0, err
+		}
+		r.lasti += r.conv.Convert()
+		for c := range r.buf {
+			d[c*nF+f] = r.buf[c]
+		}
+	}
+	return nF, nil
+}
+
+type constResampler struct {
+	*DynResampler
+	outRate freq.T
+}
+
+// SampleRate returns the output sample rate of c.
+func (c *constResampler) SampleRate() freq.T {
+	return c.outRate
+}
+
+// Resample takes a sound.Source src, a desired samplerate r, and
+// an interpolator itp.
+//
+// If itp is nil, it will default to a high quality interpolator
+// (order 10 Blackman windowed sinc interpolation).
+//
+// Resample returns a sound.Source whose SampleRate() is equal to
+// r.
+//
+// After a call to Resample, either the Receive method of src
+// should not be called, or the Receive method of the result
+// should not be called.  Clearly, the former is the usual use case.
+func Resample(src sound.Source, r freq.T, itp Itper) sound.Source {
+	sr := src.SampleRate()
+	if sr == r {
+		return src
+	}
+	tr := float64(sr) / float64(r)
+	conv := constSampleRateConverter(tr)
+	ct := NewC(src, itp)
+	dyn := NewDynResampler(ct, conv)
+	return &constResampler{DynResampler: dyn, outRate: r}
+}
+
+// NewC creates a new continuous time representation
+// of the source src using an interpolator itp.
+//
+// if itp is nil, then a default interpolator of high
+// quality will be used (order 10 Blackman windowed sinc interpolation).
+//
+// NewC calls src.Receive in this process, so src.Receive
+// should not be called if the resulting continuous time
+// interface is used.
 func NewC(src sound.Source, itp Itper) *C {
 	order := 10
 	if itp != nil {
@@ -58,7 +180,7 @@ var errMultiChanAt = errors.New("ErrMultiChanAt")
 // At returns a continuous time interpolated sample at index i.
 // It is the equivalent of
 //
-//  var buf [1]float64
+//	var buf [1]float64
 //	if err := c.FrameAt(buf[:], i); err != nil {
 //		return 0.0, err
 //	}
@@ -68,7 +190,6 @@ func (c *C) At(i float64) (float64, error) {
 	if c.Channels() != 1 {
 		return 0.0, errMultiChanAt
 	}
-	//return c.at(i)
 	var buf [1]float64
 	if err := c.FrameAt(buf[:], i); err != nil {
 		return 0.0, err
@@ -82,7 +203,13 @@ func (c *C) Channels() int {
 	return c.src.Channels()
 }
 
-// FrameAt returns a continuous time interpolated frame at index i.
+// Close closes the underlying source and returns the resulting error.
+func (c *C) Close() error {
+	return c.src.Close()
+}
+
+// FrameAt places a continuous time interpolated frame at index i in
+// dst.
 //
 // FrameAt should be called with i increasing monotonically to guarantee that c
 // does not need to go back in time arbitrarily in its underlying source.
@@ -96,11 +223,11 @@ func (c *C) Channels() int {
 // At the edges, where insufficient or no neighbors are available, the
 // interpolation is truncated symmetrically.
 //
-// FrameAt returns sound.ChannelAlignmentError if len(dst) != c.Channels().
+// FrameAt returns sound.ErrChannelAlignment if len(dst) != c.Channels().
 func (c *C) FrameAt(dst []float64, i float64) error {
 	nC := c.Channels()
 	if len(dst) != nC {
-		return sound.ChannelAlignmentError
+		return sound.ErrChannelAlignment
 	}
 	jf, jr := math.Modf(i)
 	j := int(jf)

resample/ct_test.go

@@ -121,3 +121,34 @@ func TestCLanczos(t *testing.T) {
 		}
 	}
 }
+
+func TestResampleConst(t *testing.T) {
+	fa := 800 * freq.Hertz
+	src := gen.Sin(fa)
+	r := 48000 * freq.Hertz
+	rez := Resample(src, r, nil)
+	rps := fa.RadsPerAt(rez.SampleRate())
+	rads := 0.0
+	N := 1000
+	C := 4
+	d := make([]float64, N)
+	err := 0.0
+	for c := 0; c < C; c++ {
+		n, e := rez.Receive(d)
+		if e != nil {
+			t.Fatal(e)
+		}
+		if n != N {
+			t.Fatalf("got %d frames not %d\n", n, N)
+		}
+		for i := range d {
+			ref := math.Sin(rads)
+			rads += rps
+			got := d[i]
+			err += math.Abs(got - ref)
+		}
+	}
+	if err > 0.001*float64(N)*float64(C) {
+		t.Errorf("resample error too large %f\n", err)
+	}
+}

resample/doc.go

@@ -1,6 +1,24 @@
-// Copyright 2017 The IriFrance Audio Authors. All rights reserved.  Use of
+// Copyright 2017 The ZikiChombo Authors. All rights reserved.  Use of
 // this source code is governed by a license that can be found in the License
 // file.
 
 // Package resample implements resampling/changes in resolution.
+//
+// Package resample uses interpolation for resampling which provides easy
+// control over the quality/cost tradeoff and can produce very high quality
+// resampling.  Other resampling methods may be more appropriate for a given
+// calling context, package resample doesn't yet provide other mechanisms.
+//
+// When resampling audio, any decrease in sample rate from rate S to a rate R
+// must be applied to a signal which does not contain frequencies at or above
+// R/2, or aliasing will produce strange results.
+//
+// This is often achieved by first applying a low pass filter and then
+// resampling.  As ZikiChombo does not yet have filter design support, we
+// recommend in the meantime simply taking a moving average of the signal with
+// a window size W = ceil(S/R) before decreasing the sample rate if you do not
+// have access to or knowledge about low pass filtering design.
+//
+// BUG(wsc) the shift size, effecting interpolation order limits and
+// latency of implementations is constant (64 frames).
 package resample /* import "zikichombo.org/dsp/resample" */