当前位置:网站首页>数字信号处理-9-离散余弦变换

数字信号处理-9-离散余弦变换

2022-11-24 21:47:14小何才露尖尖角

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=M1y

这就意味着如果能够高效地计算 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 也提供类似计算

参考

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

原网站

版权声明
本文为[小何才露尖尖角]所创,转载请带上原文链接,感谢
https://blog.csdn.net/weixin_40994552/article/details/128005076

随机推荐