1 波形合成

假定给一系列振幅和一系列频率，要求构建一个信号，此信号是这些频率元素的和。这样的操作就是合成

def synthesize(amps, fs, ts):
	"""
	amps 振幅数组
	fs 频率数组
	ts 采样时间点
	"""
    # ts 和 fs 的外积， m*n 矩阵
    # 每行表示 ts 的一个元素，每列表示 fs 的一个元素 
    # 每个元素表示时间和频率的乘积 ft
    args = np.outer(ts, fs)
    print(args.shape) #(11025, 4)
    # m*n 矩阵,每个元素表示 cos(2πft)
    M = np.cos(PI2*args)
    print(M.shape) #(11025, 4)
    # 点乘 m*n × n*1 = m*1
    ys = np.dot(M, amps)
    print(ys.shape) #(11025,)
    return ys
# 生成一秒时长的波形
duration=1
# 采样率
framerate=11025
fs = [100, 200, 300, 400]
amps = np.array([0.6, 0.25, 0.1, 0.005])
# 样本数
n = round(duration * framerate)
# 采样时间点
ts = np.arange(n) / framerate

synthesize(amps, fs, ts)

代码中 M 矩阵的每一行对应的时间从 0.0s 到 1.0 秒，共11025 行，表示 11025 个时间点。列表示相同时间点对应的不同频率，$np.dot(M,amps)$ 会将 M 的一行也即频率cos值乘以振幅 amps，随后加和。最后返回的 ys 即是每一个时间点不同频率cos值对应的振幅加和。

计算公式表示如下：

关于外积

2 分析

现在求上面的逆问题，也就是知道一个波形由一系列给定频率的余弦信号加和而成的，如何计算出各个频率元素对应的振幅？换种方式说就是，给定 ys、ts 和 fs, 如何算出 amps?

在合成中我们用到了如下公式：

2.1 通过线性方程组求解

反过来在逆运算中：

同样首先计算矩阵 M = cos(2πt⊗f )

随后我们需要找到 a 使得 y = Ma

M 有 11025 行， 4列(4 个频率元素)。它的每一行都是 4 个频率乘以对应的振幅加和得到的，所以我们可以取 M 的前 4 行进行计算，4 个未知数 4个方程，相当于求线性方程组。

NumPy 提供 linalg.solve 求 a

def analyze(ys, fs, ts):
    args = np.outer(ts, fs)
    M = np.cos(PI2*args)
    amps = np.linalg.solve(M, ys)
    return amps

ys = synthesize(amps, fs, ts)

# 频率个数
n = len(fs)
analyze(ys[:n], fs, ts[:n])    

out:array([0.6  , 0.25 , 0.1  , 0.005])

2.2 通过正交矩阵求解

上面解线性方程组是一种比较耗时的操作，需要的时间和 n³ 成正比，所以我们需要更好的方法。

如果 M 是方阵且可逆的话，在线性代数我们可以对 M 求逆得到 a： $a=M^{-1}y$

这就意味着如果能够高效地计算 M^-1,就能使用简单的矩阵操作得到 a。这样的耗时正比于 n²

求逆也是一个耗时的操作，但如果 M 是正交矩阵(I = MM^T)，那么 M 的转置就是 M 的逆矩阵。求转置是一个常数操作时间。

通过选择合适的 ts 和 fs, 可以让 M 成为正交矩阵。方法有多种，这也是离散余弦变换有多个版本的原因。

一种简单的方法是平移 ts 和 fs 半个单位。这个版本称为 DCT-IV，意思是 DCT 的 8 个版本 中的第 4 个。

DCT-IV

def dct_iv(ys):
    N = len(ys)
    ts = (0.5 + np.arange(N))/N
    fs = (0.5 + np.arange(N))/2
    args = np.outer(ts, fs)
    M = np.cos(PI2*args)
    amps = np.dot(M, ys)/2
    return amps

amps = np.array([0.6, 0.25, 0.1, 0.05])
N = 4.0
ts = (0.5 + np.arange(N))/N
fs = (0.5 + np.arange(N))/2
ys = synthesize(amps, fs, ts)
dct_iv(ys)

scipy.fft.dct 也提供类似计算

参考

漫画傅里叶解析
Python数字信号处理应用