一个 周期信号 分解为 若干个 正弦信号
一个 周期信号 分解为 若干个 正弦信号, 就是 傅里叶级数, 不过 我对 傅里叶级数 了解不多, 一方面 是 懒得去 细看, 一方面 也是 为了 保持 神秘感 。
我们把 一个 周期信号, 甚至 非周期信号, 记为 y = Src ( t ) , 也称为 源信号 。
傅里叶级数 的 一般形式 可以 写为 : Sin [ 1 ] + Sin [ 2 ] + Sin [ 3 ] + …… + Sin [ n ] , n -> 无穷
Sin [ n ] = An * sin ( ωn * t + ψn ) + bn , A 为 振幅, ω 为 角速度, t 为 时间, ψ 为 初始相位, b 为 增量, n 为 项 的 序号(下标), An 是 第 n 项 的 振幅, ωn 是 第 n 项 的 角速度, ψn 是 第 n 项 的 初始相位, bn 是 第 n 项 的 增量 。
那么, 将 y = Src ( t ) 展开 为 傅里叶级数 可以 这样表示 :
Src ( t ) = Sin [ 1 ] + Sin [ 2 ] + Sin [ 3 ] + …… + Sin [ n ] , n -> 无穷
只要 确定了 每一项 的 A 、ω 、ψ 、b , 就 得到 源信号 对应 的 傅里叶级数 了 。
那么, 每一项 的 A 、ω 、ψ 、b 怎么 确定 ?
记 Sins ( t ) = Sin [ 1 ] + Sin [ 2 ] + Sin [ 3 ] + …… + Sin [ n ] , n -> 无穷
可以写一个 定积分, ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] ,
[ t1, t2 ] 是 定积分 的 区间, 也是 源信号 的 区间,
| Src (t) - Sins (t) | 表示 Src (t) - Sins (t) 的 绝对值 。
我们 只要 找出 让 ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] 这个 定积分 的 值 等于 0 的 条件 就可以了 。
就是说, 我们要 为 每一项 找到 合适 的 A 、ω 、ψ 、b , 使得 ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] = 0 。
这是一个 泛函 问题 。
ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] 这个 定积分 的 意义是 Src (t) 和 Sins (t) 的 波形 的 差异 有多大, 若 这个 定积分 为 0, 则 Src (t) 和 Sins (t) 的 波形 完全相同 。
因为 Sins (t) 是 无穷级数, 由 n 个 项 组成, n -> 无穷, 所以,
ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] = 0 (1) 式
也可以写为 ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] -> 0 (2) 式
(1) 式 是 等于 0, (2) 式 是 趋于 0 。
(1) 式 和 (2) 式 是 一个 泛函方程, 也是 一个 积分方程 。
(1) 式 (2) 式 的 解 是 每一项 的 A 、ω 、ψ 、b, 可以 写成 矩阵 :
A1 , ω1 , ψ1 , b1
A2 , ω2 , ψ2 , b2
A3 , ω3 , ψ3 , b3
……
An , ωn , ψn , bn
Oh …… 终于 知道 矩阵 有什么用了 ……
一般的, 泛函 的 常见问题 是 求 最小积分条件, 所以, (1) 式 (2) 式 还可以 写成 :
ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] -> min (3) 式
表示 ʃ | Src (t) - Sins (t) | dt , [ t1, t2 ] 趋于 最小,
(3) 式 和 (1) 式 (2) 式 的 意思 差不多 。 下面 我们 就以 (1) 式 为 代表 了 。
(1) 式 这个 泛函方程 要怎么解? 不知道 。 基本上, 这种 方程 已经 达到 无从下手 的 境界 了 。
是 什么 造成了 这种 妖怪方程 ? 数学 的 抽象 。
可以 对 (1) 式 作一些 简化, 先把 Sins (t) 简化为 Sins (t) = A1 sin ( ω1 t ) + A2 sin ( ω2 t ) + A3 sin ( ω3 t ) + …… + An sin ( ωn t ) ,
这个简化 忽略了 ψ 和 b 。
于是, (1) 式 变为 :
ʃ | Src (t) - [ A1 sin ( ω1 t ) + A2 sin ( ω2 t ) + A3 sin ( ω3 t ) + …… + An sin ( ωn t ) ] | dt , [ t1, t2 ] = 0
再简化一点, 把 绝对值号 去掉,
ʃ Src (t) - [ A1 sin ( ω1 t ) + A2 sin ( ω2 t ) + A3 sin ( ω3 t ) + …… + An sin ( ωn t ) ] dt , [ t1, t2 ] = 0
积分 也 去掉 算了 ,
Src (t) - [ A1 sin ( ω1 t ) + A2 sin ( ω2 t ) + A3 sin ( ω3 t ) + …… + An sin ( ωn t ) ] = 0
A1 sin ( ω1 t ) + A2 sin ( ω2 t ) + A3 sin ( ω3 t ) + …… + An sin ( ωn t ) = Src (t)
现在 这 方程 能 解 了 吗 ? 好像还不能, 那 把 sin ( ω1 t ) , sin ( ω2 t ) , …… sin ( ωn t ) 都 换成 常量 k1, k2, …… kn, 把 Src (t) 也换成常量 S ,
k1 A1 + k2 A2 + k3 A3 + …… kn An = S , A1 , A2, A3 , …… An 为 未知数, k1, k2, k3, …… kn 为 常数, S 为 常数 (4) 式
这个 方程 能 解 了 吗? 可以 进一步简化, 让 n 为 有限大 的 自然数 ,
k1 A1 + k2 A2 + k3 A3 + …… kn An = S , A1 , A2, A3 , …… An 为 未知数, k1, k2, k3, …… kn 为 常数, S 为 常数, n 为有限大的自然数 (5) 式
(5) 式 是 一个 n 元方程, 也是一个 不定方程, 也是一个 丢番图 方程 。
(5) 式 方程 怎么解 ? 当然 (5) 式 有 无穷多个 解, 那 比如, 求 最小解, 比如 A1 + A2 + A3 + …… + An 的 和 最小 的 解 。
(4) 式 (5) 式 和 (1) 式 有 什么关系 ? 没有关系, 但可以 通过 (4) 式 (5) 式 类比 看看 (1) 式 要 怎么解 。
那么, 要怎么来 求 得 傅里叶级数 呢 ? 可以 想点 办法 。
比如, 先把 问题 简化 为 求 第一项 Sin [1] , 可以 先 求出 一个 正弦信号, 波形 和 源信号 相差 最小, 这就是 Sin [1] , 然后, 让 源信号 减去 Sin [1] , 得到 的 差 记为 z :
z = Src (t) - Sin [1]
接下来 再 求 一个 正弦信号, 波形 和 z 相差 最小, 这就是 第二项 Sin [2] , 然后, 让 z = z - Sin [2] , 以此类推, 求 Sin [3] , Sin [4] , …… , Sin [n] 。
先来看 怎么求 Sin [1] , 和 上文 (1) 式 的 原理 一样, 可以 写一个 定积分 来 表示 源信号 和 Sin [1] 的 波形 差异 大小 :
ʃ | Src (t) - Sin [1] | dt , [ t1, t2 ]
让 这个 定积分 最小,
ʃ | Src (t) - Sin [1] | dt , [ t1, t2 ] -> min (6) 式
(6) 式 就是 求 Sin [1] 的 泛函方程 。
因为 Sin [1] = A1 sin ( ω1 t + ψ1 ) + b1 , 代入 (6) 式 ,
ʃ | Src (t) - A1 sin ( ω1 t + ψ1 ) - b1 | dt , [ t1, t2 ] -> min (7) 式
对于 一些 简单 、规则 的 周期信号, 比如 方波, 三角波, 这些 信号 作为 源信号 的 话, 可以 不考虑 ψ1 、b1 ,
同时, 显然, 它们 的 基波 的 周期 和 它们 一样, 即 ω1 就是 源信号 的 ω ,
于是 (7) 式 可以 简化 为 :
ʃ | Src (t) - A1 sin ( ω t ) | dt , [ t1, t2 ] -> min (8) 式
ω 已知, 就是 源信号 的 ω , 所以, (8) 式 中 只有一个 待定系数 A1 , 只要 求出 A1 就可以, 这样的话, 可以试试 变分法, 看用 欧拉方程 能不能 解出来 。
但 欧拉方程 的 解 是 一个 函数, 这里 的 A1 是 一个 系数, 是 常量, 这就 尴尬 了, 呵呵呵呵 。 这两者(函数 和 常量) 的 矛盾 能不能 调和 ? 期待 数学天才们 的 表演, 拭目以待 。
我们 可以 画 一个 图 把 上面 的 过程 形象 的 表示出来 :
蓝色 三角波 是 源信号, 红色 正弦波 是 基波, 也就是 傅里叶级数 的 第一项 Sin [1] , 红线 和 蓝线 之间 用 绿线 标出 的 区域 就是 两者 波形 的 相差, 绿线 标出 的 区域 面积 大, 则 两者波形 相差 大, 面积小, 则 波形相差小 。
绿线 标出 的 区域 面积 就是 定积分 ʃ | Src (t) - A1 sin ( ω t ) | dt , [ t1, t2 ] , (8) 式 就是 让 这个 面积 最小 。
如图, 在这个 场景 里, 我们 只要 确定 A1, 使得 绿线 标出 的 区域 面积 最小 就可以 。
上文 提到 可以 考虑 用 变分法, 大家 可以 自己 试试, 这里 先不讨论 。 我们 看看 用 离散 线性 样本 的 方法 。
可以 大概 给 A1 划一个 范围, 把 这个 范围 切割 成 9 等分, 这样 可以 有 10 个 值, 把 这 10 个 值 代入 定积分 ʃ | Src (t) - A1 sin ( ω t ) | dt , [ t1, t2 ] , 看 哪个值 得到 的 定积分 最小, 就 取 这个 值 作为 A1, 切割 的 等分 越多, 则 结果 越 精确 。
[t1, t2] 在 这里 是 一个周期 。
当然, 对于 横坐标 t, 在 [ t1, t2 ] 区间 里 也要 切割一些 等分 来 近似计算积分 。 假设 [ t1, t2 ] 切割了 100 个 等分, 乘上 A1 的 10 个 值, 计算 的 时间复杂度 就是 10 * 1000 = 1 万 。
这个方法 对于 一般 的 源信号, 也许 有 不错 的 近似 效果 。 但是 对于 一些 源信号, 比如 突变性 的, 就 不适用 了 。
比如, 上图 蓝色 的 源信号, 突变 很明显 , 合成 它 的 基波 和 谐波, 大概 要 像 图上 这样 比较合适 。
红色 是 基波, 也就是 傅里叶级数 第一项, 橙色 是 一个 谐波, 它的 频率 是 基波 的 3 倍, 可以 认为 是 傅里叶级数 的 第三项 。
显然, 基波 和 源信号 的 波形 差异 并不是 最小, 甚至 还 有点 大 。 从这里看出, 构造 傅里叶级数 需要一些 “总揽全局” 的 规划 。
又或者说, 上面 一开始 提出 的 一项 一项 求 的 方法, 割裂了 傅里叶级数 的 “整体性” 。
如果 把 傅里叶级数 的 所有项 放到一起来 考虑, 可以 简化一点, 让 n 等于一个 有限的 不太大的 自然数, 比如 n = 5, 这样就是 由 5 个 正弦信号 来 组成 源信号 。
5 个 正弦信号 写成 矩阵 :
A1 , ω1 , ψ1 , b1
A2 , ω2 , ψ2 , b2
A3 , ω3 , ψ3 , b3
A4 , ω4 , ψ4 , b4
A5 , ω5 , ψ5 , b5
每个 正弦信号 的 A , ω , ψ , b 均 取 10 个值, 这样来匹配, 则 每个 正弦信号 会 产生出 10^4 = 1 万 个 样本,
有 5 个 正弦信号, 5 个 正弦信号 的 样本 在一起 匹配, 会 产生出 ( 1 万 ) ^ 5 = 10^20 个 样本,
10^20 个 样本 中 波形 和 源信号 差异 最小 的 那个 样本 就是 最优解 。
10 ^ 20 , 这个 计算量 太大 了 , 可以简化一点, 让 ω 固定, 5 个 正弦信号 的 ω 依次为 ω , ω * 2 , ω * 3 , ω * 4 , ω * 5 , 忽略 b 。
这样 写成 矩阵 :
A1 , ω , ψ1
A2 , 2 * ω , ψ2
A3 , 3 * ω , ψ3
A4 , 4 * ω , ψ4
A5 , 5 * ω , ψ5
因为 ω 是 固定 的, 每个 正弦信号 参与 取值 匹配 的 只有 A 、ψ, 所以, 每个 正弦信号 的 样本 是 10^2 = 100 个 ,
5 个 正弦信号 组合匹配 产生 的 样本 是 100 ^ 5 = 10^10 = 100 亿 个 。
这个 计算量 还是 太大, 呵呵 。 所以 这个 算法 似乎 有点 不科学, 也有点 白痴, 哈哈 。
到目前为止, 这个 算法 处理的, 是 简单信号 (Simple and Pure) , 简单信号 是指 波形 单一 、特点明显 的 周期信号, 常见 的 人造信号 很多是 简单信号, 比如 方波 、三角波 、梯形波 , 等等 。
简单信号 比较 容易 进行 傅里叶级数 分解 。
自然界 中 的 信号 大多 是 复杂的, 具有 复杂 和 多样 的 波形, 即使 有 周期性, 但 也不是 一成不变 的 。
自然界 中 的 信号 比如 自然界 中 的 声音, 人声 、动物的声音 、乐器的声音, 以及 自然界 中 的 各种声音 。
以 声音信号 为例, 波形 表示 音色 和 音质 , 还传达着信息, 人类语言 、动植物的声音 、大自然的声音 都 传达 着 各种各样 的 信息 。
我们来看一个 声音信号, 这个 声音信号 是 我 假想 的 。
可以 把 这个 信号 分为 3 个 信号, 也是 3 个 分量 :
分量 1 :
分量 2 :
分量 3 :
分量 1 的 频率 和 源信号 一样, 音调(音高) 听起来 大概 和 源信号 一样, 但 脉冲宽度窄, 这表示 脉冲本身的 频率高, 所以, 音色 中 带有 尖锐 的 特质 。 另外, 脉冲宽度窄 也 意味着 振动 能量 小, 所以, 听起来 会 比较 弱小 、细小, 总的来说, 听起来 是一个 尖细 的 声音 。
分量 1 合成到 源信号 中, 会 让 源信号 的 音色 听起来 更丰富 。
我感觉, 分量 1 听起来 是 “嘀 -” , 分量 2 听起来 是 “嘟 -” , 分量 3 听起来 是 “咚 -” 。 哎 ? 怎么 都是 声母 D 开头 的 ?
上面 3 个 分量 都是 简单信号, 这样 的 分量 称为 频域 的 自然分量, 简称 自然分量, 又名 特征分量 。 将 信号 分解 为 自然分量 的 方法 称为 自然分解法, 又名 特征分解法 。
自然分量 和 傅里叶级数 不一样, 傅里叶级数 是 一组 正弦函数, 且 每一个 正弦函数 的 频率 是 指定 的 。
自然分量 不是 傅里叶级数 。
要把 自然界 的 信号 分解 为 傅里叶级数 并不容易, 首先一个 问题 就是 基波 的 频率 是多少? 怎么确定 ?
自然界 的 信号 通常 是 分解为 自然分量 来 进行 分析, 以 声音信号 为例, 比如 测量 音调(音高), 分离出 基因 泛音, 声纹识别 、语音识别 、各种广义的声音识别 和 声音特征分析 。
自然界 的 信号 分解 为 傅里叶级数 反而 会 丢失 原始 的 特征信息 。
也许, 通常, 傅里叶级数 只适合于 分解 简单信号 。
将 自然界 的 信号 分解为 自然分量, 自然分量 是 简单信号, 可以 进一步 分解为 傅里叶级数 。
当然, 自然界 的 信号 不是 绝对规则的, 分解 得到 的 自然分量 的 每个周期 的 波形 和 周期 也不一定 完全一样, 可以 近似等价 为 每个周期 都一样 的 规则的 、理想的 简单信号, 再 分解 为 傅里叶级数 。
也可以 单独 取 一个 周期 的 波形 来 分解 为 傅里叶级数, 当然 也可以取 任意 一段 波形 来 分解 为 傅里叶级数 。
也可以 将 任意 一段波形 当作 一个 周期 来 分解 为 傅里叶级数, 就像 对于 非周期信号, 将 整个信号 的 定义域 作为 一个 周期 。
自然界 的 声音信号 , 可以 通过 特征分解法 分离 出 一组 自然分量, 这组 分量 中, 振幅 明显 大于 其它 分量 的 那个 分量 称为 基音, 其它 的 分量 称为 泛音 。
一个 声音 的 音调 由 基音 决定, 基音 的 音调 代表 声音 的 音调 。 泛音 的 音调 让 声音 听起来 层次丰富 。
如果 一个 声音信号 的 分量 中, 有 2 个 分量 的 振幅 不相上下, 彼此 没有 明显优势, 那么, 这可能是 2 个 声源 的 声音, 也有可能是 一个 声源 发出了 频率 不同 的 2 个 声音 混在一起 , 这 2 个 分量 可以认为 是 2 个 声音 的 基音 。
但 事情 比这 复杂, 那 2 个 声音 各自 的 泛音 要 怎么 区分, 怎么 找出来 呢 ?
这就 涉及 到 声音识别 的 问题 了, 这 涉及 到 机器识别 , 或者说 人工智能 。 比如 人类 , 或者 动物 可以 在 一堆 声音 里 识别 出 不同 声源 发出 的 声音, 可以 同时 区分 出 不同的 人 说话, 乐器, 以及 各种声音 。
这涉及到 特征识别 、特征提取, 这涉及到 人工智能 和 仿生学 。
为什么说 仿生学 呢 ? 因为 生物 有一套 “程序” 来 区分 和 提取 出 声音信号 中 的 各种声音 。
又比如, 一群 蝙蝠 在一起, 并不会 把 其它 蝙蝠 的 超声波 和 自己 发出 的 超声波 混淆起来 。
要 怎样 用 特征分解法 来 将 信号 分解为 自然分量 ? 可以看看 《卷积 毫无意义》 https://www.cnblogs.com/KSongKing/p/12839957.html 。
第一次 扫描 分离出 的, 是 频率 最高 的 分量, 甚至是一些 噪点,
第二次 扫描 分离出 的, 是 频率 第二高 的 分量,
第三次 扫描 分离出 的, 是 频率 第三高 的 分量,
……
最后一次 扫描 分离出 的, 是 频率 最低 的 分量 。
说到这里, 会想到, 有时候, 特征点 和 噪点 仅 一纸之隔 。
我写了一个 演示程序, 可以把 几个 正弦信号 合成为一个 合成信号 。 可以作为 学习测试 工具 。
程序 是 用 Html 5 + javascript 写的, 项目地址 : https://github.com/kelin-xycs/FourierStudy 。
进入 项目页面 后, 点击 右边 的 “Clone or download” 绿色按钮, 就可以 下载 项目 了 。
项目 里 的 程序文件 是 一个 Html 文件 SinesCompose.html , 用 浏览器 打开 就可以运行 。
2020-05-29 补充 :
上文 有 这样一段话 :
“
ω 已知, 就是 源信号 的 ω , 所以, (8) 式 中 只有一个 待定系数 A1 , 只要 求出 A1 就可以, 这样的话, 可以试试 变分法, 看用 欧拉方程 能不能 解出来 。
但 欧拉方程 的 解 是 一个 函数, 这里 的 A1 是 一个 系数, 是 常量, 这就 尴尬 了, 呵呵呵呵 。 这两者(函数 和 常量) 的 矛盾 能不能 调和 ? 期待 数学天才们 的 表演, 拭目以待 。
”
其实不是这么回事 。 求 A1 是 一个 函数极值问题, 不需要用到 欧拉-拉格朗日 方程 。
可以 先 把 (8) 式 里 的 定积分, 也就是 ʃ | Src (t) - A1 sin ( ω t ) | dt , [ t1, t2 ] 这个 定积分 的 表达式 求出来(如果能求出来的话), 怎么求呢 ? 先求 ʃ | Src (t) - A1 sin ( ω t ) | dt 这个 不定积分 的 表达式(如果能求出来的话), 代入 t1, t2 求 不定积分 的 差 就是 定积分, 把 定积分 记为 y = D ( A1 ) , 即 把 定积分 看作 A1 的 函数, 求 y = D ( A1 ) 的 极值, 或者说 求 y = D ( A1 ) 的 极值条件, 即 A1 = ? 时, D ( A1 ) 取 极值 。
D ( A1 ) 取 最小值 时 的 A1 , 就是 要求 的 A1 。
极值 怎么求 ? 函数极值 出现在 极值点 和 折点 。
极值点 是 导数 为 0 的 点, 且 点 的 两边 的 导数 异号 。
折点 是 导数 为 无穷, 但 函数值 不是 无穷, 且 点 的 两边 的 导数 异号 的 点 。 折点 也可以 称为 不光滑极值点 。
还有一种情况 是 单边折点, 单边折点 是 导数 为 无穷, 函数值 不是 无穷, 且 只在 点 的一边 有 函数, 另一边 没有 函数 的 点 。
比如, y = 根号 ( x ) , 当 x = 0 时, y = 0 , y ′ = 无穷 , 当 x >= 0 时, y 存在, 当 x < 0 时, y 不存在 。
所以, x = 0 是 y = 根号 ( x ) 的 单边折点 。