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 个时间点。列表示相同时间点对应的不同频率, n p . d o t ( M , a m p s ) np.dot(M, amps) 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 通过正交矩阵求解
上面解线性方程组是一种比较耗时的操作,需要的时间和 n3 成正比,所以我们需要更好的方法。
如果 M 是方阵且可逆的话,在线性代数我们可以对 M 求逆得到 a: a = M − 1 y a = M^{-1}y a=M−1y
这就意味着如果能够高效地计算 M-1,就能使用简单的矩阵操作得到 a。这样的耗时正比于 n2
求逆也是一个耗时的操作,但如果 M 是正交矩阵(I = MMT),那么 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 也提供类似计算