C# : Cublic Spline Interpolation

From:
Cubic Interpolating Splines
and
Algorithms for Cubic Spline Interpolation

private double[] CubicSpline(double[] x, double[] y, int xi_length)
        {
            int n = x.Length;
            double[] h = new double[n];
            double[] b = new double[n];
            for (int i = 0; i < n-1; i++)
            {
                h[i] = x[i + 1] - x[i];
                b[i] = (y[i + 1] - y[i]) / h[i];
            }

            double[] u = new double[n];
            double[] v = new double[n];
            u[1] = 2*(h[0] + h[1]);
            v[1] = 6*(b[1] - b[0]);
            for (int i = 2; i < n-1; i++)
            {
                u[i] = 2*(h[i-1] + h[i]) - (h[i-1]*h[i-1])/u[i-1];
                v[i] = 6*(b[i] - b[i-1]) - (h[i-1]*v[i-1])/u[i-1];
            }

            double[] z = new double[n];
            z[n-1] = 0;
            for (int i = n-2; i > 0; i--)
                z[i] = (v[i] - h[i]*z[i+1]) / u[i];
            z[0] = 0;

            double[] S = new double[xi_length];
            int j = 0;
            for (int i = 0; i < S.Length; i++)
            {
                if (i >= x[j+1] && j < x.Length - 2)
                    j++;
                double va = y[j];
                double vb = -(h[j]/6)*z[j+1] - (h[j]/3)*z[j] + (y[j+1]-y[j])/h[j];
                double vc = z[j]/2;
                double vd = (z[j+1]-z[j])/(6*h[j]);
                S[i] = va + (i-x[j])*(vb+(i-x[j])*(vc+(i-x[j])*vd));
            }

            return S;
        }

C# : So and Chan Method - R Peak Detection

Learned From:
Chiu CC, Lin TH, Liau BY : Using correlation coefficient in ECG waveform for arrhythmia detection. Biomedical Engineering - Applications, Basis & Communications 2005; 17: 147-152.

Based On:
So HH and Chan KL : Development of QRS detection method for real-time ambulatory cardiac monitor. Engineering in Medicine and Biology Society, 1997. Proceedings of the 19th Annual International Conference of the IEEE 1997; 1: 289-292.

// We need some initial data...
const double THRESHOLD_PARAM = 8;
const double FILTER_PARAMETER = 16;
const int SAMPLE_RATE = 250;

public double[] SoAndChan(double[] voltages)
{

    // initial maxi should be the max slope of the first 250 points.
    double initial_maxi = -2 * voltages[0] - voltages[1] + voltages[3] + 2 * voltages[4];
    for (int i = 2; i < SAMPLE_RATE; i++)
    {
        double slope = -2 * voltages[i - 2] - voltages[i - 1] + voltages[i + 1] + 2 * voltages[i + 2];
        if (slope > initial_maxi)
            initial_maxi = slope;
    }

    // Since we don't know how many R peaks we'll have, we'll use an ArrayList
    ArrayList rTime = new ArrayList();

    // set initial maxi
    double maxi = initial_maxi;
    bool first_satisfy = false;
    bool second_satisfy = false;
    int onset_point = 0;
    int R_point = 0;
    bool rFound = false;
    // I want a way to plot all the r dots that are found...
    int[] rExist = new int[voltages.Length];
    // First two voltages should be ignored because we need rom length
    for (int i = 2; i < voltages.Length - 2; i++)
    {

        // Last two voltages should be ignored too
        if (!first_satisfy || !second_satisfy)
        {
            // Get Slope:
            double slope = -2 * voltages[i - 2] - voltages[i - 1] + voltages[i + 1] + 2 * voltages[i + 2];

            // Get slope threshold
            double slope_threshold = (THRESHOLD_PARAM / 16) * maxi;

            // We need two consecutive datas that satisfy slope > slope_threshold
            if (slope > slope_threshold)
            {
                if (!first_satisfy)
                {
                    first_satisfy = true;
                    onset_point = i;
                }
                else
                {
                    if (!second_satisfy)
                    {
                        second_satisfy = true;
                    }
                }
            }
        }
        // We found the ONSET already, now we find the R point
        else
        {
            
                if (voltages[i] < voltages[i - 1])
                {
                    rTime.Add(i - 1);
                    R_point = i - 1;

                    // Since we have the R, we should reset
                    first_satisfy = false;
                    second_satisfy = false;

                    // and update maxi
                    double first_maxi = voltages[R_point] - voltages[onset_point];
                    maxi = ((first_maxi - maxi) / FILTER_PARAMETER) + maxi;
                }
        }
    }

    double[] results = new double[rTime.Count];

    // Now we convert the ArrayList to an array and return it
    for (int i = 0; i < rTime.Count; i++)
    {
        results[i] = rTime[i];
    }

    return results;
}

C# : Discrete Fourier Transform - Magnitude

public double[] dft(double[] x)
{
    double[] real = new double[x.Length];
    double[] imaginary = new double[x.Length];

    double[] magnitude = new double[x.Length];

    for (int i = 0; i < x.Length; i++)
    {
        real[i] = 0;
        imaginary[i] = 0;

        for (int j = 0; j < x.Length; j++)
        {
            real[i] += x[j] * Math.Cos(-2 * Math.PI * i * j / x.Length);
            imaginary[i] += x[j] * Math.Sin(-2 * Math.PI * i * j / x.Length);
        }

        magnitude[i] = Math.Sqrt(real[i]*real[i] + imaginary[i]*imaginary[i]);
    }

    return magnitude;
}