超级会员免费看
科学计算与可视化：Python 中的数值分析与信号处理
在科学计算和可视化领域，Python 提供了丰富的工具和库，本文将深入探讨其中的一些重要功能，包括非线性函数的线性回归、多项式逼近、样条插值、非线性方程求解、特殊函数以及信号处理等方面。
非线性函数的线性回归
在某些情况下，即使要拟合的函数不是线性的，仍然可以进行线性回归。以下是一个拟合指数数据的示例代码：
from pylab import *
# number of data points
N = 100
start = 0
end = 2
A = rand()+0.5
B = rand()
# our linear line will be:
# y = B*exp(A*x) = exp(A*x + log(B))
x = linspace(start, end, N)
y = exp(A*x+B)
y += randn(N)/5
# linear regression
p = polyfit(x, log(y), 1)
figure()
title(r'Linear regression with polyfit(), $y=Be^{Ax}$')
plot(x, y, 'o',
label='Measured data; A=%.2f, B=%.2f' % (A, exp(B)))
plot(x, exp(polyval(p, x)), '-',
label='Linear regression; A=%.2f, B=%.2f' % (p[0], exp(p[1])))
legend(loc='best')
show()
在这个示例中，我们通过调用
polyfit()
函数进行回归。这里传递了
x
和
log(y)
作为参数，实现了对
log(y)
的线性回归，也就是对
y
的指数回归。
多项式逼近函数
polyfit()
还可以用于通过插值来逼近函数。以逼近
sin(x)
函数为例，具体步骤如下：
1.
选择插值点
：从 0 到 $\frac{\pi}{2}$ 选择一组点作为插值点，这里选择了 0、30、45、60 和 90 度对应的弧度值，即
[0, pi/6, pi/4, pi/3, pi/2]
。
2.
计算正弦值
：使用
sqrt()
函数计算这些点的正弦值。
values = [0, pi/6, pi/4, pi/3, pi/2]
sines = sqrt(arange(5))/2
-
进行插值
：使用
polyfit()函数计算插值多项式的系数。
p = polyfit(values, sines, len(values)-1)
-
绘制误差图
：比较我们实现的
sin(x)和 Python 内置的sin(x)函数之间的差异。
figure()
x = linspace(0, pi/2, 100)
plot(x, polyval(p, x)-sin(x), label='error', lw=3)
grid()
ylabel('polyval(p, x)-sin(x)')
xlabel('x')
title('Error approximating sin(x) using polyfit()')
xlim(0, pi/2)
show()
结果显示，绝对误差小于 0.003，说明逼近效果较为准确。
样条插值
scipy.interpolate
模块提供了额外的插值函数，其中
spline(xp, yp, x)
函数可以实现分段多项式插值，以获得平滑的结果。以下是一个比较分段线性插值和样条插值的示例：
from scipy.interpolate import spline
from pylab import *
xp = linspace(0, 1, 6)
yp = sqrt(1-xp**2)
xi = linspace(0, 1, 100)
yi = interp(xi, xp, yp)
ys = spline(xp, yp, xi)
figure()
hold(True)
plot(xi, yi, '--', label='piecewise linear', lw=2)
plot(xi, ys, '-', label='spline', lw=2)
legend(loc='best')
grid()
title(r'Spline interpolation of $y=\sqrt{1-x^2}$')
xlabel('x')
ylabel('y')
axis('scaled')
axis([0, 1.2, 0, 1.2])
show()
从结果可以看出，样条插值看起来更加“平滑”。
非线性方程求解
scipy.optimize
模块提供了一些工具来求解非线性方程，这里介绍三个函数：
fsolve(f, x0)
、
bisect(f, a, b)
和
newton(f, x0)
。以计算 $\sqrt{3}$ 为例，我们构造函数 $f(x) = x^2 - 3$，然后使用这些函数求解方程 $f(x) = 0$。
def f(x):
"""Returns x**2-3"""
return x**2-3
from scipy import optimize
print(optimize.fsolve(f, 1))
print(optimize.newton(f, 1))
print(optimize.bisect(f, 1, 2))
虽然在计算 $\sqrt{3}$ 这种简单情况下，这些函数都能提供准确的结果，但它们的算法计算量较大。在大多数情况下，可以通过传递适当的参数来控制结果的精度。
特殊函数
scipy.special
模块提供了许多在高等数学和物理学中出现的特殊函数，包括贝塞尔函数、艾里函数、伽马函数、误差函数以及特殊多项式（如勒让德多项式、切比雪夫多项式）等。使用示例如下：
from scipy import special
print(special.chebyt(2))
如果需要了解更多信息，可以使用
help(special)
命令。
信号处理
信号处理是一个广泛的领域，涉及处理随时间变化的值。
scipy.signal
模块提供了一些功能，下面介绍一些基本的信号检测算法和设计滤波器的函数。
常用函数
- find(cond) ：查找数组中满足条件的元素的索引。
from pylab import *
squares = arange(10)**2
I = find(squares<50)
print(squares[I])
-
nonzero(seq)
：返回序列中的非零元素。与
find(seq!=0)类似，但nonzero()返回一个二维矩阵，而find()返回扁平化序列的索引。
A = eye(3)
I1 = find(A!=0)
I2 = nonzero(A)
A[I2] = 3
print(A)
-
where(cond, x, y)
：根据条件
cond从x和y中选择元素。
up = arange(10)
down = arange(10, 0, -1)
highest = where(up > down, up, down)
print(highest)
- select(cond, vals, default=0) ：允许设置多个条件，根据条件选择相应的元素。
up = arange(10)
ramp = select([up < 4, up > 7], [4, 7], up)
print(ramp)
简单信号检测示例
在噪声中检测信号在许多应用中都起着重要作用。以下是一个简单的示例，首先生成一个包含噪声的信号，然后进行信号检测。
生成信号 ：
from pylab import *
from scipy import signal
# parameters controlling the signal
n = 100
t = arange(n)
y = zeros(n)
num_pulses = 3
pw = 11
amp = 20
for i in range(num_pulses):
loc = floor(rand()*(n-pw+1))
y[loc:loc+pw] = signal.triang(pw)*amp
# add some noise
y += randn(n)
figure()
title('Signal and noise')
xlabel('t')
ylabel('y')
plot(t, y)
show()
信号检测 ：
# detect signals
thr = amp/2
I = find(y > thr)
# plot signal with noise plus detection
figure()
hold(True)
plot(t, y, 'b', label='signal with noise')
plot(t[I], y[I], 'ro', label='detections')
plot([0, n], [thr, thr], 'g--')
# annotate the threshold
text(2, thr+.2, 'Threshold', va='bottom')
title('Simple signal detection in noise')
legend(loc='best')
show()
峰值检测
在许多情况下，我们更关注高于阈值的最高点。以下是一个峰值检测的示例：
# peak detections
J = find(diff(I) > 1)
for K in split(I, J+1):
ytag = y[K]
peak = find(ytag==max(ytag))
plot(peak+K[0], ytag[peak], 'sg', ms=7)
波形函数
scipy.signal
模块还提供了一些波形函数，如
sawtooth()
、
square()
、
gausspulse()
和
chirp()
。以下是生成锯齿波和方波的示例：
from pylab import *
from scipy import signal
cycles = 10
t = arange(0, 2*pi*cycles, pi/10)
waveforms = ['sawtooth', 'square']
figure()
for i, waveform in enumerate(waveforms):
subplot(2, 2, i+1)
exec('y = signal.' + waveform + '(t)')
plot(t, y)
title(waveform)
axis([0, 2*pi*cycles, -1.1, 1.1])
show()
这些波形函数与之前使用的三角窗函数的区别在于它们是重复的，而
triang()
函数只生成一个窗口。
gausspulse()
和
chirp()
函数则更为专业化，具体使用方法可以参考交互式帮助。
综上所述，Python 中的这些工具和函数为科学计算和信号处理提供了强大的支持，通过合理运用它们，可以解决许多实际问题。
科学计算与可视化：Python 中的数值分析与信号处理
信号处理函数的详细分析
为了更清晰地理解前面提到的信号处理函数，我们可以通过一个表格来对比它们的功能和使用场景：
| 函数名 | 功能 | 使用场景 | 示例代码 |
| — | — | — | — |
|
find(cond)
| 查找数组中满足条件的元素的索引 | 当需要筛选出数组中符合特定条件的元素位置时使用 |
squares = arange(10)**2; I = find(squares<50); print(squares[I])
|
|
nonzero(seq)
| 返回序列中的非零元素 | 处理二维矩阵，需要找出非零元素时使用 |
A = eye(3); I2 = nonzero(A); A[I2] = 3; print(A)
|
|
where(cond, x, y)
| 根据条件
cond
从
x
和
y
中选择元素 | 当需要根据条件在两个数组中选择元素时使用 |
up = arange(10); down = arange(10, 0, -1); highest = where(up > down, up, down); print(highest)
|
|
select(cond, vals, default=0)
| 允许设置多个条件，根据条件选择相应的元素 | 有多个条件需要判断，根据不同条件选择不同值时使用 |
up = arange(10); ramp = select([up < 4, up > 7], [4, 7], up); print(ramp)
|
信号处理流程总结
下面是一个 mermaid 格式的流程图，展示了从信号生成到信号检测和峰值检测的完整流程：
graph LR
A[生成信号参数设置] --> B[生成时间向量 t 和初始值向量 y]
B --> C[随机放置三角脉冲]
C --> D[添加噪声]
D --> E[绘制信号与噪声图]
E --> F[设置阈值进行信号检测]
F --> G[绘制检测结果图]
G --> H[进行峰值检测]
H --> I[绘制峰值检测结果图]
不同插值方法的对比
在前面我们介绍了多项式插值和样条插值，这里我们通过一个表格来对比它们的特点：
| 插值方法 | 优点 | 缺点 | 适用场景 |
| — | — | — | — |
| 多项式插值 | 实现简单，能通过已知插值点 | 当插值点较多时，可能会出现龙格现象，导致插值函数在区间边缘波动较大 | 插值点较少，对插值函数的光滑性要求不高的情况 |
| 样条插值 | 分段多项式插值，结果平滑 | 计算相对复杂 | 需要得到光滑插值函数的情况 |
非线性方程求解方法的选择
对于
scipy.optimize
模块中的
fsolve(f, x0)
、
bisect(f, a, b)
和
newton(f, x0)
这三个求解非线性方程的函数，我们可以根据不同的情况来选择使用：
-
fsolve(f, x0)
：
-
适用场景
：当我们有一个初始猜测值
x0
，并且希望找到函数
f(x)
的根时可以使用。它可以处理较为复杂的非线性函数。
-
操作步骤
：
1. 定义函数
f(x)
。
2. 选择一个接近根的初始猜测值
x0
。
3. 调用
optimize.fsolve(f, x0)
函数求解。
-
bisect(f, a, b)
：
-
适用场景
：当我们知道根所在的区间
[a, b]
时使用。该方法要求函数在区间
[a, b]
上连续，并且
f(a)
和
f(b)
的符号不同。
-
操作步骤
：
1. 定义函数
f(x)
。
2. 确定根所在的区间
[a, b]
。
3. 调用
optimize.bisect(f, a, b)
函数求解。
-
newton(f, x0)
：
-
适用场景
：当我们有一个初始猜测值
x0
，并且函数
f(x)
可导时使用。牛顿法收敛速度较快，但可能会出现不收敛的情况。
-
操作步骤
：
1. 定义函数
f(x)
。
2. 选择一个接近根的初始猜测值
x0
。
3. 调用
optimize.newton(f, x0)
函数求解。
特殊函数的应用拓展
scipy.special
模块中的特殊函数在许多领域都有广泛的应用，例如：
-
贝塞尔函数
：在波动问题、热传导问题、电磁学等领域有重要应用。例如，在圆柱形波导中，贝塞尔函数可以用来描述电磁波的传播模式。
-
伽马函数
：在概率论、统计学、物理学等领域经常出现。例如，在计算某些概率分布的积分时，会用到伽马函数。
-
误差函数
：在统计学、物理学、工程学等领域有应用。例如，在描述正态分布的累积分布函数时，会用到误差函数。
波形函数的实际应用
前面介绍的波形函数，如
sawtooth()
、
square()
、
gausspulse()
和
chirp()
，在信号处理和通信领域有很多实际应用：
-
锯齿波（
sawtooth()
）
：常用于音频合成、示波器测试等。在音频合成中，锯齿波可以产生尖锐的声音效果。
-
方波（
square()
）
：在数字电路、通信系统中广泛应用。例如，在数字信号传输中，方波可以表示二进制的 0 和 1。
-
高斯脉冲（
gausspulse()
）
：在雷达系统、无线通信中用于产生短脉冲信号。高斯脉冲具有良好的时域和频域特性。
-
线性调频信号（
chirp()
）
：在雷达、声纳等系统中用于目标检测和测距。线性调频信号的频率随时间线性变化，可以提高系统的分辨率。
通过合理运用 Python 中的这些工具和函数，我们可以在科学计算、信号处理等领域解决各种实际问题，并且可以根据具体需求选择合适的方法和函数，以达到最佳的效果。
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
