fouriertransform.cs

心电信号处理库文件

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

namespace DSProcessing
{
    /// <summary>
    /// FFT, a part of DSProcessing library.
    /// 
    /// based on: AForge.NET
    ///      adress: http://code.google.com/p/aforge/
    ///      author: Andrew Kirillov
    ///      
    /// modifyed under GNU GPL by Jan Sova
    /// mailto: [e-mail]
    /// a part of DSProcessing library
    /// -------------------------------------------------------------------------
    ///
    /// DSProcessing - C#/C++ library of signal processing, speech processing,
    ///              and communications classes and functions
    ///
    /// Copyright (C) 2007-2008
    ///
    /// This program is free software; you can redistribute it and/or modify
    /// it under the terms of the GNU General Public License as published by
    /// the Free Software Foundation; either version 2 of the License, or
    /// (at your option) any later version.
    ///
    /// This program 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 General Public License for more details.
    ///
    /// You should have received a copy of the GNU General Public License
    /// along with this program; if not, write to the Free Software
    /// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
    ///
    ///-------------------------------------------------------------------------
    /// </summary>
    /// 
    /// <remarks>The class implements one dimensional and two dimensional
    /// Discrete and Fast Fourier Transformation.</remarks>
    /// 
    public class FourierTransform
    {
        /// <summary>
        /// Tools class.
        /// </summary>
        /// 
        public class Tools
        {
            /// <summary>
            /// Calculates power of 2
            /// </summary>
            /// 
            /// <param name="power">Power</param>
            /// 
            /// <returns>Returns specified power of 2 in the case if power is in the range of
            /// [0, 30]. Otherwise returns 0.</returns>
            /// 
            public static int Pow2(int power)
            {
                return ((power >= 0) && (power <= 30)) ? (1 << power) : 0;
            }

            /// <summary>
            /// Checks if the specified integer is power of 2
            /// </summary>
            /// 
            /// <param name="x">Integer number to check</param>
            /// 
            /// <returns>Returns <b>true</b> if the specified number is power of 2.
            /// Otherwise returns <b>false</b>.</returns>
            /// 
            public static bool IsPowerOf2(int x)
            {
                return (x & (x - 1)) == 0;
            }

            /// <summary>
            /// Get base of binary logarithm
            /// </summary>
            /// 
            /// <param name="x">Source integer number</param>
            /// 
            /// <returns>Power of the number (base of binary logarithm).</returns>
            /// 
            public static int Log2(int x)
            {
                if (x <= 65536)
                {
                    if (x <= 256)
                    {
                        if (x <= 16)
                        {
                            if (x <= 4)
                            {
                                if (x <= 2)
                                {
                                    if (x <= 1)
                                        return 0;
                                    return 1;
                                }
                                return 2;
                            }
                            if (x <= 8)
                                return 3;
                            return 4;
                        }
                        if (x <= 64)
                        {
                            if (x <= 32)
                                return 5;
                            return 6;
                        }
                        if (x <= 128)
                            return 7;
                        return 8;
                    }
                    if (x <= 4096)
                    {
                        if (x <= 1024)
                        {
                            if (x <= 512)
                                return 9;
                            return 10;
                        }
                        if (x <= 2048)
                            return 11;
                        return 12;
                    }
                    if (x <= 16384)
                    {
                        if (x <= 8192)
                            return 13;
                        return 14;
                    }
                    if (x <= 32768)
                        return 15;
                    return 16;
                }

                if (x <= 16777216)
                {
                    if (x <= 1048576)
                    {
                        if (x <= 262144)
                        {
                            if (x <= 131072)
                                return 17;
                            return 18;
                        }
                        if (x <= 524288)
                            return 19;
                        return 20;
                    }
                    if (x <= 4194304)
                    {
                        if (x <= 2097152)
                            return 21;
                        return 22;
                    }
                    if (x <= 8388608)
                        return 23;
                    return 24;
                }
                if (x <= 268435456)
                {
                    if (x <= 67108864)
                    {
                        if (x <= 33554432)
                            return 25;
                        return 26;
                    }
                    if (x <= 134217728)
                        return 27;
                    return 28;
                }
                if (x <= 1073741824)
                {
                    if (x <= 536870912)
                        return 29;
                    return 30;
                }
                return 31;
            }
        }
        /// <summary>
        /// Fourier transformation direction
        /// </summary>
        public enum Direction
        {
            /// <summary>
            /// Forward direction of Fourier transformation
            /// </summary>
            Forward = 1,

            /// <summary>
            /// Backward direction of Fourier transformation
            /// </summary>
            Backward = -1
        };

        /// <summary>
        /// One dimensional Discrete Fourier Transform
        /// </summary>
        /// 
        /// <param name="data">Data to transform</param>
        /// <param name="direction">Transformation direction</param>
        /// 
        public static void DFT(Complex[] data, Direction direction)
        {
            int n = data.Length;
            double arg, cos, sin;
            Complex[] dst = new Complex[n];

            // for each destination element
            for (int i = 0; i < n; i++)
            {
                dst[i] = Complex.Zero;

                arg = -(int)direction * 2.0 * System.Math.PI * (double)i / (double)n;

                // sum source elements
                for (int j = 0; j < n; j++)
                {
                    cos = System.Math.Cos(j * arg);
                    sin = System.Math.Sin(j * arg);

                    dst[i].Re += (data[j].Re * cos - data[j].Im * sin);
                    dst[i].Im += (data[j].Re * sin + data[j].Im * cos);
                }
            }

            // copy elements
            if (direction == Direction.Forward)
            {
                // devide also for forward transform
                for (int i = 0; i < n; i++)
                {
                    data[i].Re = dst[i].Re / n;
                    data[i].Im = dst[i].Im / n;
                }
            }
            else
            {
                for (int i = 0; i < n; i++)
                {
                    data[i].Re = dst[i].Re;
                    data[i].Im = dst[i].Im;
                }
            }
        }
        /// <summary>
        /// Two dimensional Discrete Fourier Transform
        /// </summary>
        /// 
        /// <param name="data">Data to transform</param>
        /// <param name="direction">Transformation direction</param>
        /// 
        public static void DFT2(Complex[,] data, Direction direction)
        {
            int n = data.GetLength(0);	// rows
            int m = data.GetLength(1);	// columns
            double arg, cos, sin;
            Complex[] dst = new Complex[System.Math.Max(n, m)];

            // process rows
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < m; j++)
                {
                    dst[j] = Complex.Zero;

                    arg = -(int)direction * 2.0 * System.Math.PI * (double)j / (double)m;