在学习3Blue1Brown关于傅里叶变换的视频后,想要尝试重复他所进行的实验。本文中的代码仅使用了numpy库和matplotlib库
信号的波形图
先从最简单的信号波形图开始。时间区间以及信号向上平移的默认值采用了视频中相应的近似值。fl 是频率的列表,表示这些频率的余弦信号的叠加;b 是每个余弦信号向上平移的距离;l,r 为时间范围。
def drawSignal(fl, b=1.1, l=0, r=4.5):
x = np.linspace(l, r, 1000)
y = np.zeros(len(x)) + b*len(fl)
for f in fl:
y += np.cos(2*np.pi*f*x)
plt.plot(x, y)
plt.show()
如图是绘制 2Hz+3Hz 的波形图,drawSignal([2, 3])
缠绕图像
然后是缠绕图像的绘制。f 是缠绕频率;fl 是原始信号的频率列表;n 是缠绕圈数,比如 n = 4,将缠绕 4 圈,线圈转过 8 π \pi π;bound 是绘制出的图像的坐标轴范围;b 是每个余弦信号向上平移的距离;sz 是采样的总数。
def drawWire(fw, fl, n=4, bound=2.5, b=1.1, sz=10000):
plt.axis('equal')
plt.xlim((-bound, bound))
plt.ylim((-bound, bound))
t = n*1/fw # 总共经过了 nT 的时间
x = np.linspace(0, t, sz)
y = np.zeros(len(x)) + b*len(fl)
for fs in fl:
y += np.cos(2*np.pi*fs*x)
theta = -np.linspace(0, 2*n*np.pi, sz) # 总共转过了 2*n*pi 弧度
wx = y * np.cos(theta) # 通过参数方程得到x,y
wy = y * np.sin(theta)
plt.plot(wx, wy)
plt.show()
如图是绘制 f = 3 的缠绕图像,drawWire(3, [3], 1)
缠绕图像的周期问题
在绘制这个图像的时候遇到一个问题,缠绕图像经过多少个周期,它的图像就完整了,并开始重复了呢?
缠绕周期
T
w
=
1
f
w
T_w=\frac1{f_w}
Tw=fw1 原始信号周期
T
s
=
1
f
s
T_s=\frac1{f_s}
Ts=fs1 如果满足周期性,则有
M
T
w
=
N
T
s
MT_w = NT_s
MTw=NTs 其中
M
,
N
M, N
M,N 为正整数,即求最小的
M
M
M 使得
M
T
w
T
s
=
M
f
s
f
w
\frac{MT_w}{T_s} = \frac{Mf_s}{f_w}
TsMTw=fwMfs 为正整数。
考虑到本实验中原始信号的频率可以设为正整数,而缠绕频率希望能够连续变化。可以令
f
w
=
p
q
f_w = \frac pq
fw=qp 其中
p
,
q
p,q
p,q 为互素的正整数,则等价于求最小的
M
M
M 使得
M
f
s
q
p
\frac{Mf_sq}{p}
pMfsq 为正整数。
注意到
f
s
q
f_sq
fsq 在该假设下是正整数,不妨设
g
g
g 为
f
s
q
f_sq
fsq 和
p
p
p 的最大公因数,且
p
=
g
∗
h
p = g*h
p=g∗h,则
M
f
s
q
p
=
M
f
s
q
/
g
∗
g
g
∗
h
=
M
f
s
q
/
g
h
\frac{Mf_sq}{p} = \frac{Mf_sq/g*g}{g*h} = \frac{Mf_sq/g}{h}
pMfsq=g∗hMfsq/g∗g=hMfsq/g 其中
f
s
q
/
g
f_sq/g
fsq/g 是正整数,所以
M
M
M 可以取最小值
h
h
h
例如刚才的例子 f = 3, fl = [3],M 取 1即可。另外,原始信号可能是多个频率的余弦信号的组合,只需要对每一个信号,都计算一个M,并取其中最大的那一个,则可以得到最终的 M.
代码实现如下:
# 计算最大公因数
def gcd(p, q):
r = p % q
while r != 0:
p = q
q = r
r = p % q
return q
# 计算最小的重复周期数
def minN(p, q, fl):
n = 0
for f in fl:
n = max(n, p // gcd(f*q, p))
return n
这样就可以方便地画出完整图像:
p = 2
q = 5
fl = [3]
n = minN(p, q, fl)
drawWire(p/q, fl, n)
质心变化图像
在绘制质心时,一开始出了差错;后来发现错误原因是我计算的是完整图像的质心,相当于金属丝在无限地缠绕;而视频中给出的是一定时间内的缠绕图像,即金属丝缠绕出的可能是完整图像,也可能是不完整的。
比如
f
w
=
1.03
f_w = 1.03
fw=1.03, 取
n
=
t
0
∗
f
w
n = t_0 * f_w
n=t0∗fw 那么实际缠绕的时间为
n
T
=
t
0
∗
f
w
/
f
w
=
t
0
nT = t_0*f_w/f_w=t_0
nT=t0∗fw/fw=t0 得到的图像可能是这样的:
但这个缠绕频率需要很多个周期才能显示出周期性,缠绕12圈的图像如下(仍然不完整):
而我估计出视频中采用的时间长度大概是4.4,实现代码如下:
# 计算质心 返回结果
def calcWire(f, fl, n=4, b=1.1, sz=10000):
t = n*1/f
x = np.linspace(0, t, sz)
y = np.zeros(len(x)) + b*len(fl)
for fs in fl:
y += np.cos(2*np.pi*fs*x)
theta = -np.linspace(0, 2*n*np.pi, sz)
wx = y * np.cos(theta)
wy = y * np.sin(theta)
return np.mean(wx)
# 在一个频率范围采样若干个频率 得到一组(频率,质心横坐标)数据
def drawCenter(fl, b=1.1, draw=True):
q = 100 # 每单位Hz采样的数量
x = np.linspace(1/q, 5, 5*q)
xl = []
for p in range(1, 5*q+1):
xl.append(calcWire(p/q, fl, p/q*4.4, b=b, sz=100000))
xl = np.array(xl)
if draw:
plt.plot(x, xl)
plt.show()
return x, xl
fl=[3] 时的质心变换图像如下,drawCenter([3])
在f=3处确实有一个尖峰,而且其纵坐标值确实为0.5左右
fl=[2,3] 时的质心变换图像如下,drawCenter([2,3])
确实有f=2, f=3两个尖峰
质心变化图像的叠加
对于fl = [2,3]的质心变换图像,其图像是fl=[2],fl=[3]两个质心变换图像的和,实验如下:
def drawAdding():
x2, xl2 = drawCenter([2], b=1.1, draw=False)
x3, xl3 = drawCenter([3], b=1.1, draw=False)
x23, xl23 = drawCenter([2, 3], b=1.1, draw=False)
plt.subplot(211)
plt.plot(x2, xl2+xl3)
plt.subplot(212)
plt.plot(x2, xl23)
plt.show()
结果如下,上图为分解信号的质心变换的和,下图是复合信号的质心变换。
感觉对傅里叶变换又更理解了一些呢。