Python 利用插值法实现数值重现——blog9

数字重现是指在对同样的数据或实验进行多次重复时，得到相似或相同的结果。

目录

基本介绍

实现步骤

设置模拟时长和采样点数

生成信号

采样

傅立叶变换

可视化输出

运行结果

讨论与思考

---

想象一下你正在研究一个新的药物，你希望知道它是否能够治疗某种疾病。为了测试这个药物的有效性，你会进行多次实验。如果每次实验的结果都很相似，即药物在不同实验中都有相似的治疗效果，那就是数字重现。这意味着你的实验是可靠的并且结果是可信的。

为什么数字重现很重要呢？因为它能够验证科学研究的可靠性。如果某个研究只是在单次实验中得到了某些结果，那么这个结果的可信度就不高。但是，如果这个结果在多次重复的实验中反复出现，那么我们就可以相信这个结果是稳定和可靠的。这就是为什么科学研究经常要求进行多次实验并进行数字重现。

数字重现还有助于排除偶然性的误差。有时候，一个实验的结果可能受到一些偶然的因素影响，例如实验条件、测量误差等。通过进行多次实验，我们可以看到这些偶然因素的影响是否被消除，以确保结果的可靠性。

基本介绍

本文的这个实验主要利用插值法进行抽样信号的复原，对于插值法，不了解的小伙伴可以上网查查，百度上给出的定义是：

插值法，也称为“内插法”，是一种利用函数在某个区间中已知的若干点的函数值，构造一个特定的函数，以便在该区间的其他点上估计函数值的近似方法。这个特定的函数通常是多项式，因此被称为多项式插值。插值法可以分为几种类型：

1.线性插值。使用两个相邻数据点之间的直线进行插值，这种方法简单直观，但在数据点之间可能会产生不连续或不光滑的曲线。

2.多项式插值。使用多项式函数来拟合数据点之间的曲线，包括拉格朗日插值和牛顿插值。拉格朗日插值构造一个多项式函数，使得该函数在已知数据点上与原始函数完全一致，保证插值曲线通过已知数据点，但在高次插值中可能会出现振荡现象。牛顿插值则是通过构造一个递推的多项式，同样保证插值曲线通过已知数据点，且在高次插值中不容易出现振荡现象。

3.样条插值。这是一种更平滑的插值方法，通过使用分段连续的低次多项式来逼近数据，从而得到一条平滑的曲线。最常见的样条插值方法是三次样条插值，它在给定的数据点上构造一组三次多项式来实现数据的插值，可以更好地逼近实际曲线，并且可以通过调整边界条件来控制曲线的性质。

实现步骤

不过python作为面向对象编程语言，使用起来倒也不需要管这么多，scipy有自带的interpolate方法，轻松实现你想要的插值方法，以下是具体实现步骤：

设置模拟时长和采样点数

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import time
frequency = 2
duty_cycle = 0.4
pulse_width = 0.3
T = 10/frequency   # 设置模拟时长
N = 10000          # 设置采样点数
f = 1/T

生成信号

def creat_sig():
    x = np.linspace(0, 10, num=11, endpoint=True)
    y = np.sin(2*np.pi*f*x) + np.random.normal(loc=0, scale=0.3, size=11)
    return x, y

采样

def sampling(x, y):
    f_linear = interp1d(x, y, kind='linear')
    x_new = np.linspace(0, 10, num=101, endpoint=True)
    y_linear = f_linear(x_new)
    f_cubic = interp1d(x, y, kind='cubic')
    y_cubic = f_cubic(x_new)
    return y_cubic, y_linear

傅立叶变换

start = time.time()

x, y = creat_sig()
y1, y2 = sampling(x, y)
x_ = np.linspace(0, 10, 101)

freqs = np.fft.fftfreq(N, T/N)   # 计算频率轴

fft_result_1 = np.fft.fft(y1)/N      # 进行FFT变换，归一化处理
amplitudes_1 = 2 * np.abs(fft_result_1[:N//2])  # 获取振幅信息（仅需使用正半轴数据）

fft_result_2 = np.fft.fft(y2)/N      # 进行FFT变换，归一化处理
amplitudes_2 = 2 * np.abs(fft_result_1[:N//2])  # 获取振幅信息（仅需使用正半轴数据）

可视化输出

plt.subplot(511)
plt.title("original signal")
plt.plot(x, y)
x_0 = np.linspace(0, 10, 5000)
plt.subplot(512)
plt.title("linear interpolate")
plt.plot(x_, y2, 'r')
plt.subplot(513)
plt.title("cubic interpolate")
plt.plot(x_, y1, 'b')

plt.subplot(514)
plt.title("fft of linear one")
plt.plot(x_, amplitudes_2)

plt.subplot(515)
plt.title("fft of cubic one")
plt.plot(x_, amplitudes_1)

plt.show()

end = time.time()

print(f"总耗时：{end - start}秒")

运行结果

讨论与思考

import numpy as np
import matplotlib.pyplot as plt

t1 = np.linspace(0, 2, 2000)
y1 = 10 * np.sin(2 * np.pi * 505 * t1) + np.sin(2 * np.pi *504 * t1) + np.sin(2 * np.pi *506 * t1)

t2 = np.linspace(0, 2, 20000)
y2 = 10 * np.sin(2 * np.pi * 505 * t2) + np.sin(2 * np.pi *504 * t2) + np.sin(2 * np.pi *506 * t2)

plt.subplot(211)
plt.plot(t1, y1)

plt.subplot(212)
plt.plot(t2, y2)

plt.show()

k3 = 5
k4 = 10
k5 = 11
k6 = 20

t3 = t2[1:len(t2) + 1:k3]
y3 = y2[1:len(y2) + 1:k3]

t4 = t2[1:len(t2) + 1:k4]
y4 = y2[1:len(t2) + 1:k4]

t5 = t2[1:len(t2) + 1:k5]
y5 = y2[1:len(t2) + 1:k5]

t6 = t2[1:len(t2) + 1:k6]
y6 = y2[1:len(t2) + 1:k6]


plt.subplot(411)
plt.plot(t3, y3)
plt.title("interval:5")

plt.subplot(412)
plt.plot(t4, y4)
plt.title("interval:10")

plt.subplot(413)
plt.plot(t5, y5)
plt.title("interval:11")

plt.subplot(414)
plt.plot(t6, y6)
plt.title("interval:20")

plt.show()

要求：运行以上代码，分析其结果含义，欢迎在评论区留言交流。

觉得有帮助的小伙伴还请点个关注

后续会持续分享 免费、高质量的高校相关以及Python学习文章

（拒绝AI水文章）