傅里叶变换学习笔记

在学习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=\frac{1}{f_w}$ 原始信号周期 $T_s=\frac{1}{f_s}$ 如果满足周期性，则有 $M{T}_w=N{T}_s$ 其中 $M,N$ 为正整数，即求最小的 $M$ 使得 $\frac{M{T}_w}{T_s}=\frac{M{f}_s}{f_w}$ 为正整数。
考虑到本实验中原始信号的频率可以设为正整数，而缠绕频率希望能够连续变化。可以令 $f_w=\frac{p}{q}$ 其中 $p,q$ 为互素的正整数，则等价于求最小的 $M$ 使得 $\frac{M{f}_{s}q}{p}$ 为正整数。
注意到 $f_{s}q$ 在该假设下是正整数，不妨设 $g$ 为 $f_{s}q$ 和 $p$ 的最大公因数，且$p=g\ast h$，则 $\frac{M{f}_{s}q}{p}=\frac{M{f}_{s}q/g\ast g}{g\ast h}=\frac{M{f}_{s}q/g}{h}$ 其中 $f_{s}q/g$ 是正整数，所以 $M$ 可以取最小值 $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$, 取 $n=t_0\ast f_w$ 那么实际缠绕的时间为 $nT=t_0\ast f_w/f_w=t_0$ 得到的图像可能是这样的：

但这个缠绕频率需要很多个周期才能显示出周期性，缠绕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()

结果如下，上图为分解信号的质心变换的和，下图是复合信号的质心变换。

感觉对傅里叶变换又更理解了一些呢。