numpy 使用Python绘制正弦曲线

q3aa0525  于 2023-10-19  发布在  Python
关注(0)|答案(2)|浏览(140)

我在编写一个通用程序时遇到了一个问题,如果给定频率和曲线的总持续时间,该程序可以自动调整绘制正弦曲线所需的点数**。
根据Mariagboffi在另一个SO问题上提供的建议,我编写了以下通用代码。

import numpy as np
import matplotlib.pyplot as plt
import math

# Take formula delta-t <= 1/(5*omega) from SO answer "https://stackoverflow.com/a/76895268/2386113"
frequency = 0.02 #Two HZ
delta_t = math.floor( 1/(5*2*np.pi*frequency))

total_duration = 300
num_points = math.ceil(total_duration/delta_t)
t = np.linspace(0, total_duration, num_points)
sine_t = np.sin(t)

#Only for test: to plot a line of constant amplitude
straight_line = t/t

# red for numpy.sin()
plt.plot(t, sine_t, color = 'red')
plt.plot(t, straight_line, color = 'blue')
plt.title("numpy.sin()")
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

**问题:**上面的代码对于frequency = 2 Hz工作正常,但是如果我设置了frequency = 0.02 Hz,那么我又得到了一个错误的图(可以通过增加硬编码的点数来纠正)。我如何根据频率和持续时间计算点的数量?
频率= 2 Hz(正确)

频率= 0.02 Hz(错误)

更新:

根据注解中的chrslg建议(这是有道理的),我更正了代码并使用了np.sin(t*2*np.pi*frequency)。代码的更新部分是(注意频率现在设置为2):

frequency = 2 #Two HZ
delta_t = max(math.floor( 1/(5*2*np.pi*frequency)), 1)

total_duration = 300
num_points = math.ceil(total_duration/delta_t)
t = np.linspace(0, total_duration, num_points)
sine_t = np.sin(2*np.pi*frequency*t) #sin(omega * t)

由于上面的代码,我得到了下面的正弦曲线,这又是错误的,因为frequency = 2意味着每秒应该有2个正弦波,因此在300秒的持续时间内应该有600个正弦波,而下面的截图中不是这样的:

mpgws1up

mpgws1up1#

我认为你需要取两倍于频率的样本。这是奈奎斯特采样定理的结果(更多信息here和这里)。简单地说,如果周期函数的频率为f/Hz,而你想完全重建周期波,你将需要至少2f每秒才能完成。正弦函数是一个周期函数。
如果你有Hz的频率,你可以使用下面的代码片段(对你的原始代码做了轻微的修改):

import numpy as np
import matplotlib.pyplot as plt
import math

frequency_hz = 1

total_duration = 4
multiplicator = 40  # nyquist recommends to samples at least 2*frequency in Hz per second
num_points = math.floor(total_duration * frequency_hz * multiplicator) + 1
print(num_points)
t = np.linspace(0, total_duration, num_points)
sine_t = np.sin(t * 2 * np.pi * frequency_hz)

# Only for test: to plot a line of constant amplitude
straight_line = t / t

# red for numpy.sin()
plt.plot(t, sine_t, color="red")
plt.plot(t, straight_line, color="blue")
plt.title("numpy.sin()")
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

生成此图:

我建议将multiplicator更改为一个较大的数字,如20,以便始终有一个清晰的图。
我使用了floor函数,使得num_points是一个整数。+1num_points并不是真的需要。
我不确定0.02是什么单位,但我认为这段代码也可以适应那个单位!

hsvhsicv

hsvhsicv2#

我不知道为什么建议5ωT的点数(T是总持续时间)。这在一定程度上是有道理的,因为它与要绘制的周期数成正比。因此,如果5ωT对于ω的一个值是可以的,那么对于任何其他值都应该是可以的。这意味着每一个完整的周期,无论在任何时间间隔中有多少个周期,都是由31.416点绘制的。
这确实是足够的,看到正弦。但是.416仍然意味着每个周期不是在相同的相对横坐标上绘制的。这意味着,对你来说更重要的是,如果我正确地得到了你想要的,每个周期的最高点并不总是相同的(有时,它几乎完全在sin(π/2)=1结束。但是,下一个周期峰值由sin(π/2+2π31/31.416)sin(π/2+2π32/31.416)表示。
delta_T计算中使用floor部分地解决了这个问题,确保在每个周期内以整数个点结束(因此,每个周期具有相同的采样,无论好坏,因此整个正弦看起来具有恒定的幅度)
但是,floor也是你现在的问题所在。因为,对于任何足够大的频率(超过0.031,听起来并不那么大),delta_t将四舍五入为0。我想这就是为什么你添加了max(...,1),它可以防止除以0的错误,但不是问题:你的delta_t是一个非常粗略的近似值。
如果你想每个周期有一个恒定的采样数,那么,就这样做吧。

K=8 # K=4 is the strict minimum to see 0, 1, 0, -1, repeating
# K=a factor of 4 is needed if you want two of the samples to be 1 and -1, that is to see the whole amplitude

num_points = math.ceil(K*frequency*total_duration)
# either frequency*total_duration is an integer (that is there is an 
# exact number of periods in the interval, which was the case in both of
# your examples. 300×0.02 = 6 complete periods. 300×2 = 600 of them)
# and then, floor is useless (except for typing reasons)
# or it isn't, and then, just using `math.ceil` doesn't solve the problem: 
# it isn't the correct sampling: we need to also adjust total_duration, so 
# that an exact number of kth of period is represented by `num_points`
total_duration = num_points/K/frequency

# Another important point: you need `num_points` from 0 included to
# total_duraction excluded. Or else you don't distributed evenly your points
# (total_duration is the 1st point of the next period, not the last of this one)
# Or, if you really need to include `total_duration` abscissa, then,
# you need `num_points+1` not `num_points` points
t = np.linspace(0, total_duration, num_points, False) # False meaning "excluded"
# or
t = np.linspace(0, total_duration, num_points+1, True) # True, meaning included, is the default, I explicit it here to show the difference

sine = np.sin(2*np.pi*frequency*t)
plt.plot(t, sine)
plt.show()

极简K=4, frequency=0.02结果x1c 0d1x
K=12,频率=0.02。由于K是4的因子,我们仍然看到完整的-1,1振幅,有更多的点

K=5,频率=0.02。由于K不是4的因子,因此在最高点处的采样不精确地处于π/2相位,也就是说,我们不具有完整的-1,1幅度。但至少,它总是相同的,因为每个周期正好有5个点。

K=32,频率=0.02。这是非常接近2×5π你似乎已经从另一个答案。但至少它是准确的。和4的倍数。

K=32,频率=0.017。由于在这个频率下,这不是(0,300)间隔中的确切周期数(甚至不是周期的1/32),因此我们调整该间隔,以便每个周期仍然具有确切的样本数(因此我们在这里实际上从0到301.47绘制,以便仍然能够具有每个周期K=32个样本的重复模式)

相关问题