python生成不对称正弦波与傅里叶级数分解

lilili~

已于 2022-03-23 16:15:51 修改

阅读量4.7k

点赞数 2

分类专栏： python 文章标签： python 机器学习

于 2022-03-15 19:22:13 首次发布

本文链接：https://blog.csdn.net/weixin_46536094/article/details/123509404

版权

python 专栏收录该内容

3 篇文章 1 订阅

订阅专栏

不对称正弦波介绍

不对称正弦波图像如下
在这里插入图片描述
参数方程可以表示为
$\begin{array}{l} x=t-k \sin t \\ y=\sin t \end{array}$
t为参数方程的参数，k表示不对称的程度

python生成

求解思路

目标是画出一个x与y的函数图像，因此先通过x求解非线性方程得到t，而后通过t求y,其中求解非线性方程使用scipy.optimize.fsolve

python代码

#导入包
import numpy as np
from scipy.optimize import fsolve, curve_fit
import matplotlib.pyplot as plt
#中间函数，通过x求解t
def feixianxing(t):
    global buffer  # buffer就是x的取值，在asin函数中修改buffer，反解参数t
    k = -.4  # 参数k
    return t - k * np.sin(t) - buffer
def asin(x):
    funy = []
    global buffer
    for m in x:
        buffer = m
        result = fsolve(feixianxing, 0)
        buffery = np.sin(result[0])
        funy.append(buffery)
    return funy

以上函数定义完成后，可以设置x的取值范围，并且生成y

x = np.linspace(0, 10., 100)
y = asin(x)
plt.plot(x, y, 'r-', label='original data')
plt.show()

图像如下：
在这里插入图片描述

傅里叶分解

由于这个函数是周期奇函数，所以可以用傅里叶级数对其进行分解，使用的函数为scipy.optimize.curve_fit,代码如下：

# 模板函数，因为是奇函数，所以没有cos的傅里叶项，period是周期，可以自己定义
def fourier(x, *a):
    global period
    tau = period / 2
    f = a[0] * np.sin(np.pi / tau * x)
    for deg in range(1, len(a)):
        f += a[deg] * np.sin((deg + 1) * np.pi / tau * x)
    return f
# 使用curce_fit函数
period = 2 * np.pi
x = np.linspace(0, 10., 100)
y = asin(x)
#x，y是拟合的数据集，[1.0]*3是傅里叶的项数，函数中的意思是估计拟合参数的初值是三个1
popt, pcov = curve_fit(fourier, x, y, [1.0] * 3)
#画出函数并比对
plt.plot(x, y, 'r-', label='original data')
plt.plot(x, fourier(x, *popt), "g--", label='fit curve')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

图像如下，拟合程度可以通过调节项数来改变。
在这里插入图片描述

完整代码

import numpy as np
from scipy.optimize import fsolve, curve_fit
import matplotlib.pyplot as plt


# 中间函数，通过x求解t
def feixianxing(t):
    global buffer  # buffer就是x的取值，在asin函数中修改buffer，反解参数t
    k = -.4  # 参数k
    return t - k * np.sin(t) - buffer


def asin(x):
    funy = []
    global buffer
    for m in x:
        buffer = m
        result = fsolve(feixianxing, 0)
        buffery = np.sin(result[0])
        funy.append(buffery)
    return funy


def fourier(x, *a):
    global period
    tau = period / 2
    f = a[0] * np.sin(np.pi / tau * x)
    for deg in range(1, len(a)):
        f += a[deg] * np.sin((deg + 1) * np.pi / tau * x)
    return f


period = 2 * np.pi
x = np.linspace(0, 10., 100)
y = asin(x)
popt, pcov = curve_fit(fourier, x, y, [1.0] * 3)
plt.plot(x, y, 'r-', label='original data')
plt.plot(x, fourier(x, *popt), "g--", label='fit curve')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()