Hilbert.cs

using System;
using System.Collections.Generic;
using System.Text;

// Copyright (c) 2006, 2007 by Hugh Pyle, inguzaudio.com

namespace DSPUtil
{
    /// <summary>
    /// FIR filter with a Dirac pulse at its center
    /// </summary>
    public class Dirac : SoundObj
    {
        private int _length;
        
        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="length">Length (preferably odd)</param>
        public Dirac(int length)
        {
            if (length % 2 == 0)
            {
                length++;
            }
            _length = length;
            base.NumChannels = 1;
        }

        /// <summary>Length of the filter</summary>
        public override int Iterations
        {
            get { return _length; }
        }

        /// <summary>
        /// Number of channels: always 1
        /// </summary>
        public override ushort NumChannels
        {
            get
            {
                return 1;
            }
            set
            {
                if (value != 1)
                {
                    throw new ArgumentOutOfRangeException();
                }
            }
        }

        /// <summary>
        /// Get an iterator for samples of the filter
        /// </summary>
        public override IEnumerator<ISample> Samples
        {
            get
            {
                int mid = (int)((_length - 1) / 2);
                for (int j = 0; j < _length; j++)
                {
                    double val = 0;
                    int n = j - mid;
                    if (n == 0)
                    {
                        val = 1;
                    }
                    yield return new Sample(val);
                }
            }
        }
    }

    /// <summary>
    /// FIR filter approximating a Hilbert transform; not windowed
    /// </summary>
    public class Hilbert : SoundObj
    {
        private int _length;

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="length">Length (preferably odd)</param>
        public Hilbert(int length)
        {
            if (length % 2 == 0)
            {
                length++;
            }
            _length = length;
            base.NumChannels = 1;
        }

        /// <summary>Length of the filter</summary>
        public override int Iterations
        {
            get { return _length; }
        }

        /// <summary>
        /// Number of channels: always 1
        /// </summary>
        public override ushort NumChannels
        {
            get
            {
                return 1;
            }
            set
            {
                if (value != 1)
                {
                    throw new ArgumentOutOfRangeException();
                }
            }
        }

        /// <summary>
        /// Get an iterator for samples of the filter
        /// </summary>
        public override IEnumerator<ISample> Samples
        {
            get
            {
                int mid = (int)((_length - 1) / 2);
                for (int j = 0; j < _length; j++)
                {
                    int n = j - mid;
                    double val = 0;
                    if (n % 2 != 0)
                    {
                        val = 2 / (Math.PI * n);
                    }
                    yield return new Sample(val);
                }
            }
        }
    }


    /// <summary>
    /// Given a signal input
    /// Return the (positive only) magnitude Hilbert envelope
    /// </summary>
    public class HilbertEnvelope : SoundObj
    {
        private int _length;
        private ISoundObj _i;
        private ISoundObj _r;

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="length">Length for Hilbert FIR (preferably odd)</param>
        public HilbertEnvelope(int length)
        {
            if (length % 2 == 0)
            {
                length++;
            }
            _length = length;
            base.NumChannels = 1;

            // Imaginary portion is generated by a Hilbert transform of the real input
            _i = new FastConvolver(new Hilbert(length));

            // Real portion is a delayed version of the real input
            _r = new FastConvolver(new Dirac(length));
        }

        public override IEnumerator<ISample> Samples
        {
            get
            {
                if (_input == null)
                {
                    yield break;
                }

                _i.Input = _input;
                _r.Input = _input;

                ushort nc = _input.NumChannels;

                IEnumerator<ISample> ienum = _i.Samples;
                IEnumerator<ISample> renum = _r.Samples;

                bool imore = ienum.MoveNext();
                bool rmore = renum.MoveNext();
                while (imore && rmore)
                {
                    ISample icurr = ienum.Current;
                    ISample rcurr = renum.Current;
                    ISample ret;
                    if (nc == 2)
                    {
                        ret = new Sample2();
                    }
                    else
                    {
                        ret = new Sample(nc);
                    }
                    for (ushort c = 0; c < nc; c++)
                    {
                        ret[c] = new Complex(rcurr[c], icurr[c]).Magnitude;
                    }
                    yield return ret;

                    imore = ienum.MoveNext();
                    rmore = renum.MoveNext();
                }
            }
        }
    }

    /// <summary>
    /// Given a signal X
    /// produce a filter returning (a + jb)X
    /// where 'a' is the in-phase multiplier, 'b' is the quadrature multiplier, and 'j' means a 90 degree phase shift
    /// </summary>
    public class PhaseMultiplier : SoundObj
    {
        private double _a;
        private double _b;
        private int _length;
        private double _f;

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="phase">in-phase and quadrature multipliers</param>
        /// <param name="length">Length (preferably odd)</param>
        public PhaseMultiplier(Complex phase, int length)
        {
            _init(phase, length, 0, 0);
        }

        public PhaseMultiplier(Complex phase, int length, uint sampleRate)
        {
            _init(phase, length, sampleRate, 0);
        }

        public PhaseMultiplier(Complex phase, int length, uint sampleRate, double cornerFreq)
        {
            _init(phase, length, sampleRate, cornerFreq);
        }

        private void _init(Complex phase, int length, uint sampleRate, double cornerFreq)
        {
            if (length % 2 == 0)
            {
                length++;
            }
            _a = phase.Re; // inphase scale
            _b = phase.Im; // quadrature scale
            _length = length;
            SampleRate = sampleRate;
            _f = cornerFreq;
            base.NumChannels = 1;
        }

        public double PhaseRe
        {
            get
            {
                return _a;
            }
            set
            {
                _a = value;
            }
        }
        public double PhaseIm
        {
            get
            {
                return _b;
            }
            set
            {
                _b = value;
            }
        }

        /// <summary>Length of the filter</summary>
        public override int Iterations
        {
            get { return _length; }
        }

        /// <summary>
        /// Number of channels: always 1
        /// </summary>
        public override ushort NumChannels
        {
            get
            {
                return 1;
            }
            set
            {
                if (value != 1)
                {
                    throw new ArgumentOutOfRangeException();
                }
            }
        }

        public override IEnumerator<ISample> Samples
        {
            get
            {
                int mid = (int)((_length - 1) / 2);

                // Make a Hilbert transformation filter (for the 'j' component)
                Hilbert hilbert = new Hilbert(_length);

                // And a window for the filter
                BlackmanHarris window = new BlackmanHarris(mid, mid);

                if (_f > 0)
                {
                    FilterProfile fp = new FilterProfile();
                    fp.Add(new FreqGain(0, 0));
                    fp.Add(new FreqGain(_f / 2, 0));
                    fp.Add(new FreqGain(_f * 2, -200));
                    FilterImpulse fi = new FilterImpulse(_length, fp, FilterInterpolation.COSINE, _sr);
                    FastConvolver fc = new FastConvolver(hilbert, fi);
                    SoundBuffer sb = new SoundBuffer(fc);
                    window.Input = sb.Subset(_length / 2, _length);
                }
                else
                {
                    window.Input = hilbert;
                }

                // Return the samples
                int j=0;
                foreach(ISample sample in window)
                {
                    int n = j - mid;
                    if(n==0)
                    {
                        // Set the in-phase component
                        sample[0] = _a;
                    }
                    else
                    {
                        // Scale the quadrature component
                        sample[0] = sample[0] * _b;
                    }
                    yield return sample;
                    j++;
                }
            }
        }
    }
}