信号变换成傅里叶级数的迭代算法

Question

我想向你们展示一种迭代算法，把信号转换成傅里叶级数。当我运行单元测试时，我获得了以下结果：

• 一次迭代:1秒，通过；
• 十次迭代:3秒，通过。

我使用VisualStudioCommunity2015环境。我正在用C#和MathNet NuGet包编写这个库。

我如何优化这个算法，你对这个库有什么看法？

单元测试：

[TestMethod]
    public void OneIterationTest()
    {
        double[] points = new double[628];
        FourierSeries fourier = new FourierSeries((x) => Math.Sign(Math.Sin(x)), 1);
        for (int current = 0; current < 628; current++)
        {
            points[current] = fourier.GetHarmonicWavePoint(current);
        }
    }

    [TestMethod]
    public void TenIterationTest()
    {
        double[] points = new double[628];
        FourierSeries fourier = new FourierSeries((x) => Math.Sign(Math.Sin(x)), 10);
        for (int current = 0; current < 628; current++)
        {
            points[current] = fourier.GetHarmonicWavePoint(current);
        }
    }

这就是我的实现：

using MathNet.Numerics.Integration;
using System;

namespace SignalProcessingLibrary
{
    /**
     * <summary>
     * The assumption of following class is distribution signals to Fourier Series.
     * First, we calculate the coefficient a_0 by integrate function in -T/2 to T/2 interval
     * and multiply it by 2/T. Next, we apply iterative algorithm to superposition of harmonic waves.
     * At the beginning this algorithm, we calculate coefficient a_n and b_n and sum it using
     * math formula. We substitute to n index of the loop and loop is performed until reach N.
     * </summary>
     */
    public class FourierSeries
    {
        /**
         * <summary>
         * Initalizes fields by default values.
         * </summary>
         */
        public FourierSeries() : this((x) => Math.Sign(Math.Sin(x)), 10, 2 * Math.PI)
        {
        }

        /**
         * <summary>
         * Sets function and iterations. This leaves the default value of _period.
         * </summary>
         * <param name="function">Analyzed function</param>
         * <param name="iterations">Degree of approximation</param>
         */
        public FourierSeries(Func<double, double> function, double iterations) : this(function, iterations, 2 * Math.PI)
        {
        }

        /**
         * <summary>
         * Sets all fields.
         * </summary>
         * <param name="function">Analyzed function</param>
         * <param name="iterations">Degree of approximation</param>
         * <param name="period">Period of analyzed function</param>
         */
        public FourierSeries(Func<double, double> function, double iterations, double period)
        {
            Function = function;
            Iterations = iterations;
            Period = period;
        }

        /**
         * <summary>
         * Implements a_n coefficient.
         * </summary>
         * <param name="n">Number of iteration</param>
         * <returns>a_n coefficient</returns>
         */
        private double GetA(double n)
        {
            return (2 / Period) *
                GaussLegendreRule.Integrate(
                (x) => (Function.Invoke(x) * Math.Cos((2 * n * Math.PI * x) / Period)),
                -Period / 2,
                Period / 2,
                1024);
        }

        /**
         * <summary>
         * Implements b_n coefficient.
         * </summary>
         * <param name="n">Number of iteration</param>
         * <returns>b_n coefficient</returns>
         */
        private double GetB(double n)
        {
            return (2 / Period) *
                GaussLegendreRule.Integrate(
                (x) => (Function.Invoke(x) * Math.Sin((2 * n * Math.PI * x) / Period)),
                -Period / 2,
                Period / 2,
                1024);
        }

        /**
         * <summary>
         * Implements a_0 coefficient.
         * </summary>
         * <returns>a_0 coefficient</returns>
         */
        private double GetA0()
        {
            return (2 / Period) *
                GaussLegendreRule.Integrate(
                    (x) => (Function.Invoke(x)),
                    -Period / 2,
                    Period / 2,
                    1024);
        }

        /**
         * <summary>
         * This does superposition of harmonics waves.
         * </summary>
         * <param name="x">Desired point</param>
         * <returns>Superposition of harmonics waves</returns>
         */
        public double GetHarmonicWavePoint(double x)
        {
            double sum = this.GetA0() / 2;
            for (int i = 0; i < Iterations; i++)
            {
                sum = sum +
                    (GetA(i) * Math.Cos((2 * i * Math.PI * x) / Period) +
                    GetB(i) * Math.Sin((2 * i * Math.PI * x) / Period));
            }
            return sum;
        }

        public Func<double, double> Function { get; set; }
        public double Iterations { get; set; }
        public double Period { get; set; }
    }
}

Answer 1

可以计算一次多次计算的所有值。我不确定这是否对演出有任何重大影响.最耗时的调用可能是GaussLegendreRule.Integrate，但是值得一试。

这需要将传递给构造函数的迭代和周期变为只读。但无论如何，这样做是明智的。

检查传递给(public)构造函数的参数的值对该类是否有效也是很好的做法。

这将将代码更改为以下内容：

public class FourierSeries
{
    private readonly double _iterations;
    private readonly double _period;
    private readonly double _periodDividedByTwo;
    private readonly double _twoDividedByPeriode;
    private readonly Func<double, double> _function;

    private const double TAU = Math.PI * 2;

    public FourierSeries() : this((x) => Math.Sign(Math.Sin(x)), 10, TAU)
    {}

    public FourierSeries(Func<double, double> function, double iterations) : this(function, iterations, TAU)
    {}

    public FourierSeries(Func<double, double> function, double iterations, double period)
    {
        if (function == null) throw new ArgumentNullException("function");
        if (iterations < 1) throw new ArgumentException("iterations must be greater than 0")
        if (period < 1) throw new ArgumentException("period must be greater than 0")

        _function = function;
        _iterations = iterations;
        _period = period;
        _periodDividedByTwo = period / 2;
        _twoDividedByPeriode = 2 / period;
    }

    private double GetA(double n)
    {
        var factor = TAU * n;
        return _twoDividedByPeriode *
            GaussLegendreRule.Integrate(
            (x) => (_function.Invoke(x) * Math.Cos((factor * x) / _period)),
            - _periodDividedByTwo,
            _periodDividedByTwo,
            1024);
    }
    private double GetB(double n)
    {
        var factor = TAU * n;
        return _twoDividedByPeriode *
            GaussLegendreRule.Integrate(
            (x) => (_function.Invoke(x) * Math.Sin((factor * x) / _period)),
            -_periodDividedByTwo,
            _periodDividedByTwo,
            1024);
    }

    private double GetA0()
    {
        return _twoDividedByPeriode *
            GaussLegendreRule.Integrate(
                (x) => (_function.Invoke(x)),
                -_periodDividedByTwo,
                _periodDividedByTwo,
                1024);
    }

    public double GetHarmonicWavePoint(double x)
    {
        double sum = this.GetA0() / 2;
        for (int i = 0; i < _iterations; i++)
        {
            sum += GetA(i) * Math.Cos((2 * i * Math.PI * x) / _period) +
                   GetB(i) * Math.Sin((2 * i * Math.PI * x) / _period);
        }
        return sum;
    }
}

我认为GaussLegendreRule.Integrate使用了一个迭代集成算法，这是相当昂贵的。因此，另一个优化可能是在可能的情况下以分析的方式集成功能--但这取决于您的用例.