傅里叶里傅

前言

本来只打算简单介绍下傅里叶变换的, 但是一开始写文章, 就发现自己很多点都没有搞清楚. 陆陆续续花了两周多的时间才感觉自己是真的明白了. 果然以教为学是最好的学习方法.

这篇文章没有非常严谨的推导证明, 只是用我觉得 intuitive 的理解来介绍相关知识. 因为傅里叶分析主要用于信号处理, 所以在网上查找资料时, 涉及到很多信号处理的相关知识, 什么采样定理, 梳状函数等等, 我会尽量不去涉及这些相关知识 (因为我也不懂啊).

看这篇文章前最好对傅里叶变换有个大概的了解, 推荐看这篇文章: 傅里叶分析之掐死教程

概述

前言中推荐的那篇文章对于傅里叶分析的介绍已经很好了, 文章中涉及到了几个概念, 1) 傅里叶级数 (FS) 2) 傅里叶变换 (FT) 3) 离散傅里叶变换 (DFT), 如果查找相关资料你还会看到 4) 离散时间傅里叶变换 (DTFT), 5) 离散傅里叶级数 (DFS) 6) 快速傅里叶变换 (FFT), 7) 连续时间傅里叶变换 (CTFT) 等等, 我当时的疑惑点就是这些概念彼此之间究竟是什么关系? 离散的是什么? 连续的是什么? 离散时间又是什么意思?

上面提到的几个概念中有些其实描述的是同一个东西, 只是名字不一样. 一般说 傅里叶变换 就是指 连续时间傅里叶变换, 通常也认为 离散傅里叶级数 和 离散傅里叶变换 是同一个概念, 最后, 快速傅里叶变换 只是 离散傅里叶变换 的一种快速算法实现. 最终上面的概念就剩下四个 1) 傅里叶级数 2) 傅里叶变换 3) 离散傅里叶变换 4) 离散时间傅里叶变换. 傅里叶分析的核心思想就是对一个函数进行分析, 对于一个函数我们可以从多个维度去描述它, 傅里叶分析关心的是函数它连续还是离散? 周期还是非周期? 两两组合就分别对应我们那四个概念.

  连续 离散
周期 傅里叶级数 离散傅里叶变换
非周期 傅里叶变换 离散时间傅里叶变换

离散时间的意思是指信号在时域上离散, 换句话说就是待分析的函数是离散的. 离散傅里叶变换 没有 时间 这个限定词, 就意味着它既在时间上离散, 还在其它地方离散, 和时间对应的就是频率. 后面我们会看到, 在频域上离散就表示待分析的函数是周期的.

傅里叶级数 (Fourier Series)

先给出定义: f(x) 是周期为 T 的函数

这个公式的含义就是, 任何周期函数可以转换成一系列的正余弦函数的累加. 对于这个公式, MIT 6.024 线性代数公开课里给出的解释我觉得是最漂亮的. 任何一个 n 维的空间中, 一定可以找到一组正交基, 空间中所有向量都可以由这组正交基进行线性运算来得到. 所谓正交基, 就是一组两两正交的基向量. 在二维和三维空间中, 两向量点乘为 0 就说明这两个向量垂直. 将这一性质拓展到 n 维空间就是, 两个向量的点乘为 0 则两向量正交.

现在我们把函数当成一个无限维空间 (希尔伯特空间) 中的向量, 类比有限维空间中点乘的公式, 定义希尔伯特空间的点乘的公式为:

我们可以发现傅里叶级数中的各项, 两两之间点乘的结果为 0, 所以傅里叶级数就是希尔伯特空间中的一组正交基.

利用傅里叶级数正交这一性质, 求各项的系数就非常简单了, 公式两边各乘上要求系数的项, 然后积分就好了 (因为等式两边均为周期为 T 的函数, 所以只需要在周期 T 内积分就可以了):

p.s. 项除以 2 的原因是, 如果不除以 2, 其它各项乘以自身之后积分的结果都是 , 而 1 在周期内的积分是 T, 这里为了统一, 所以 需要除以 2.

傅里叶级数还有另外一个形式. 根据欧拉公式, 我们有 , 所以:

将上式代入原本的傅里叶级数, 我们有:

我们可以把之前求的 带入, 也可以利用正交性质来求解. 需要注意的是, 涉及到复数时, 向量点乘的定义为 , 其中的 表示共轭转置 (Hermitian), 其实就是把向量中的 i 变成 -i.

最终得到:

傅里叶变换 (Fourier Transform)

傅里叶级数要求 f(x) 是周期函数, 那如果不是周期函数是不是就没法处理了呢? 换个角度看, 非周期函数其实就是周期 T 无穷大的的周期函数.

所以对于傅里叶变化, 我们只需要把 带入到上面的傅里叶级数中就可以了.

为了方便, 我们定义 , 其中 表示频率.

时:

将上面这些带入之前的傅里叶级数有:

因为 , 原先 中的系数 变成 “挪” 到了 f(x) 中, 同时将 f(x) 中的求和符号换成了积分符号.

离散时间傅里叶变换 (Discrete-Time Fourier Transform)

之前我们看的都是连续的函数, 但是现实中我们获取到的都是离散的数值, 或者说是数字化了的信号. 离散时间傅里叶变换就是离散化的傅里叶变换, f(x) 变成了一组离散的数 f[k] 其中 k 为采样点序号.

怎么处理呢? 想像一下, 所谓的连续函数不过就是采样间隔无穷小的函数, 令每秒采样数为 N:

因为 f[k] 是离散的, 所以把 中的积分符号换成求和符号, 注意, dx 不能去掉, 为什么? 因为 dx 就是个常量, 一个无穷小的常量, 和积分是没有关系的.

分析下 你会发现它是一个周期函数, 也就是说 . 为什么? 首先 是个周期为 的函数, 要让 成立, 就需要对于每个 k, 我们有 , 显然只要 为非零整数即可. 我一开始这里没有明白, 就很疑惑, 为什么同一个函数, 离散的在频域是周期的而连续的在频域不是周期的, 现在看来就很简单了, 因为 f(x) 中的 x 是实数, 你没办法找到一个数与任意实数相乘都为非零整数, 而对于离散的函数, 你只要 是非零整数即可, 因为 k 本身就是整数.

dx “挪” 到 f[k] 中, 上面式子就变成:

现在还有个问题, 前面分析过了, 是周期函数, 所以显然 也是个周期函数, 对一个周期函数, 在 上积分, 显然值只可能是: 无穷大, 无穷小或者零. 这种情况下, 再乘以一个无穷小的 dx, 最终的值究竟是什么呢? 我们来看下这个 dx 究竟是什么. 我们反过来看, 是个连续周期函数, 那么我们是不是也可以对它进行傅里叶级数分解呢? 所以:

可以看到 , 最后我们利用 dx 的周期为 这一性质, 可以将在 的积分变成:

所以离散时间傅里叶变换的最终结果就是:

离散傅里叶变换 (Discrete Fourier Transform)

有了离散时间傅里叶变换的基础, 离散傅里叶变换就很简单了, 只需要将 转换回 即可:

同样的道理, 利用 f[k] 周期为 NT 的性质, 同时将 “挪” 到 中:

最终离散傅里叶变换的结果为:

我们可以看到, 离散傅里叶变换的形式才是傅里叶分析的一般形式, 要连续傅里叶变换就让 , 要非周期就让 .

所谓连续不过就是采样间隔无穷小的离散.

快速傅里叶变换 (Fast Fourier Transform)

快速傅里叶变换就不多说了, 就是基于离散傅里叶变换的线性代数形式的一些特征, 将计算复杂度由 降为 的一种算法.

其他

可以发现时域的函数和频域的函数在傅里叶分析中是非常对称的, 比如时域连续, 在频域就是非周期, 在时域非周期, 在频域就是连续. 写这篇文章的时候, 时域和频域的表达形式, 我也刻意用了相对称的表述方式.

最后, 我们可以将最开始的表格扩充下:

时域 连续 离散  
周期 傅里叶级数 离散傅里叶变换 离散
非周期 傅里叶变换 离散时间傅里叶变换 连续
  非周期 周期 频域

p.s. 对于带复数的函数求导, 只需要将 i 当成一个常数即可

epoll 比 select 高效的原因

epoll 比 select 高效的原因还是要从他们的实现上来说.

首先, linux 下文件描述符都可以在其状态就绪 (有新的数据, 可以写入数据等) 时, 执行一个指定的回调函数. select 和 epoll 的实现都是基于这一机制, 区别就在于如何利用这一回调机制.

select

当我们调用 select, 并传入一批文件描述符, 假设当前所有文件描述符都不就绪, select 会在所有这些文件描述符上加上一个回调函数, 然后将当前进程睡眠. 回调函数的功能很简单, 就是唤醒睡眠的进程. 进程唤醒后继续执行, 会轮询检查监控的文件描述符是否就绪, 同时将文件描述符上对应的回调函数移除.

epoll

和 select 不同的地方是, epoll 不需要频繁的在监控的文件描述符上添加移除回调函数. 只有当你调用 epoll_ctl 添加监控的文件描述符时 epoll 才会在对应的文件描述符上添加一个回调函数. 相应的, 调用 epoll_ctl 删除监控的文件描述符时会将回调函数移除. 另外需要说明的是, 调用 epoll_create 的时候, epoll 会创建一个文件描述符, 当你调用 epoll_wait 的时候, 会在 epoll 创建的文件描述符上插入一个回调函数, 然后将当前进程睡眠.

epoll 插入到监控的文件描述符上的回调函数所做的事情是, 将当前就绪的文件描述符放到 epoll 自己维护的一个红黑树里, 同时执行 epoll 自己创建的文件描述符上的回调函数, 这个回调函数会将睡眠的进程唤醒. 唤醒的进程会继续执行 epoll_wait 函数, 将之前插入的回调函数移除.

总结

可以发现, epoll 比 select 高效的原因是, a) 不需要频繁的在监控的文件描述符上添加移除回调函数 b) 可以直接获取就绪的文件描述符, 不需要遍历检查每个文件描述符

等待队列

上面的 select 和 epoll 的实现都基于文件描述符在就绪时执行回调函数, 那么这个又是怎么实现的呢?

其实也很简单, 每个文件描述符维护了一个等待队列 (wait_queue_head). 所有的文件描述符都对外暴露一个 file_operations 的结构体 (具体实现是各个 IO 设备的驱动实现的), 其中有一个 poll 方法 (和系统调用 poll 不是同一个东西). 调用这个 poll 方法, 会在等待队列中插入一个等待队列项 (wait_queue), 等待队列项中就有我们要执行的回调函数. 这样当对应的文件描述符状态改变 (写入, 读取等操作) 时就会遍历等待队列并执行相应的回调函数.

其他

用惯了高层语言, 看 select 和 epoll 的源码的时候, 总是会带入一些已经抽象过的概念, 比如进程. 对于内核代码来说, 进程不过就是内存中的一块数据而已, 所有的进程调度, 都是需要内核代码自己实现.

很久不用 c 了, 对 c 的很多语法都忘了. 比如一开始好奇 epoll 的回调函数里怎么拿到对应的文件描述符的, 后来发现原来是因为 wait_queue 是属于一个结构体的, 回调函数中直接基于偏移获取 wait_queue 所属的结构体, 再获取结构体中的文件描述符, 这个不就是闭包的底层实现嘛!

看源码的时候总是看了后面忘前面, 后来发现还是要把代码图形化, 画 UML 图是个好方法. 所以记忆还是要基于图形化!

源码解读Linux等待队列

源码解读poll/select内核机制

源码解读epoll内核机制

ETF 知识点整理

ETF (Exchange-traded Fund) 是跟踪指数的基金. 如何跟踪的呢? 基金严格按照指数的股票组成比例进行投资, 一旦基金的价格与跟踪的股票出现价差, 套利者就会将价差抹平.

净值是基金单位份额对应的资产的价值. 要注意净值和市场价不是同一个概念, 两者不一定相等. 买入卖出的价格叫做市场价, 认购赎回的价格叫做净值.

分级基金 (Structured Func) 和普通的基金一样也是把基金投资于各种资产获得收益. 只是人为的将基金 (也叫做母基金) 分成 A, B 两个子基金, 其中 A 拿的是固定收益, B 的收益取决于母基金整体的收益. 当母基金获利时, 收益扣除 A 的固定收益后就是 B 的收益. 当母基金亏损时, B 的亏损就是 A 的固定收益加上整体的亏损. 可以理解为 B 是向 A 约定利息借钱来投资母基金. 为了防止 B 亏损太多 (通常为亏损 75%, 也就是母基金亏损 37.5%) 无法支付 A 的收益, 会对 B 进行折算, 也就是强制平仓. 因为 A 的收益是固定收益, 所以 A 的净值为 . A 的收益以母基金份额的形式发放. 在折算时, 和股票分红一样, A 的市场价是要做除权的.

折算的本质就是将基金的净值恢复到 1, 所以折算都是针对净值来算的. 分级基金会在固定日期, 基金净值上涨和下跌的时候进行折算, 分别叫做 定折, 上折 和 下折. 以下折为例, 甲持有 B 一万份, 现在 B 净值 0.25. 折算后, B 净值 1, 甲持有 B 两千五百份.

分级基金的运作原理,定折、下折、上折操作,以及分级A、分级B的套利原理

死锁

定义

一组进程处于死锁状态是指: 如果在一个进程集合中的每个进程都在等待只能由该集合中的其他一个进程才能引发的事件, 则称一组进程或系统此时发生了死锁.

必要条件

1971 年 Coffman 总结出了系统产生死锁的四个 必要条件:

  1. 互斥条件 (mutual exclusion)

    任一时刻, 一个资源仅能被一个进程独占. 另一个进程请求已被占用的资源时, 它将被设置成等待状态, 直到占用者释放资源.

  2. 占有和等待条件 (hold and wait)

    一个进程因请求资源得不到满足而等待时, 不释放已占有的资源.

  3. 不剥夺条件 (no preemption)

    任一进程不能从另一进程那里抢夺资源.

  4. 循环等待条件 (circular wait)

    存在一个循环等待链.

注意了, 这里是死锁产生的必要条件, 不是充分条件. 也就是说当这四个条件满足时, 并不是说系统此时一定发生死锁了, 但是当系统发生死锁时, 这四个条件一定是满足的.

解决

可以从三个方面来解决死锁问题,它们是 死锁防止 (deadlock prevention); 死锁避免 (deadlock avoidance); 死锁检测和恢复 (deadlock detection and recovery).

比较难区分的是死锁防止和死锁避免, 我的理解是, 死锁防止是 客观 上防止四个必要条件满足, 死锁避免是 主观 上避免四个必要条件满足.

如果我们把死锁比喻成在高速公路上开车撞到障碍物导致死亡, 那么死锁防止等同于政府派环卫工人去清理障碍物; 死锁避免是驾驶员通过驾驶技巧来规避障碍物; 死锁检测与恢复是撞倒障碍物之后通过安全气囊防止死亡.

死锁防止

层次化分配是一种有效的死锁防止策略. 简单来说就是将资源排序, 请求资源时只能按照顺序请求. 这样就破坏了四个必要条件中的循环等待条件.

以银行转账为例, 我们可以让程序必须按照用户 ID 大小顺序来锁定账户. (呃, 这个好像是主观避免了)

死锁避免

银行家算法是最有名的死锁避免算法, 但是因为算法要求进程提前知道自己所需的最大资源数量, 这就导致该算法并不具有多大的实际价值. 算法并不复杂, 直接上代码:

import numpy as np

n_processes = int(input('Number of processes? '))
n_resources = int(input('Number of resources? '))

available_resources = [int(x) for x in input('Claim vector? ').split(' ')]

currently_allocated = np.array([[int(x) for x in input('Currently allocated for process ' + str(i + 1) + '? ').split(' ')] for i in range(n_processes)])
max_demand = np.array([[int(x) for x in input('Maximum demand from process ' + str(i + 1) + '? ').split(' ')] for i in range(n_processes)])

total_available = available_resources - np.sum(currently_allocated, axis=0)

running = np.ones(n_processes) # An array with n_processes 1's to indicate if process is yet to run

while np.count_nonzero(running) > 0:
	at_least_one_allocated = False
	for p in range(n_processes):
		if running[p]:
			if all(i >= 0 for i in total_available - (max_demand[p] - currently_allocated[p])):
				at_least_one_allocated = True
				print(str(p) + ' is running')
				running[p] = 0
				total_available += currently_allocated[p]
	if not at_least_one_allocated:
		print('Unsafe')
		exit()
				
print('Safe')

死锁检测和恢复

这是一种比较有实际意义的解决方案, 当检测到系统发生死锁了, 可以强制破坏上述四个必要条件, 比如, 让优先级低的进程释放已占用的资源.

总结

总而言之, 解决死锁就从四个必要条件出发.

操作系统教程

Banker’s algorithm

Makefile 生成文件追踪

很多大型开源项目使用 make 来构建, 有时候我们想知道一个文件是怎么被构建出来的, 或者说是从哪些文件构建出来的, 知道这些对于找到程序的入口是很有帮助的. 那么该怎么做呢?

第一反应肯定是直接看 Makefile, 但是大型项目的 Makefile 都很复杂, 各种 include 和变量. 虽然 make 提供 -n-p 选项, 但是因为构建命令里还会调用 make 命令, 当使用 -n 选项时, 这些命令并不会执行, 所以看不到最终执行的所有指令.

那么该怎么办呢? 其实很简单, 直接构建然后从构建输出的日志里来寻找答案, 因为 make 默认会将构建时执行的命令打印出来 (如果使用了 -s 选项或者命令加了 @ 前缀就不会输出).

下面我们以查看 OpenJDK 中 bin/java 是如何构建出来的为例来看下具体的操作:

  1. 首先将项目完整的构建一遍

  2. 将你想要追踪的文件删除, 这里我们将 bin/java 删掉

  3. 重新构建一遍项目同时将输出重定向到一个文件里

    步骤 1 和 2 主要是为了减少构建日志输出的量从而方便搜索

  4. 最后就是从输出日志里搜索相关的信息

bin/java 为例, 我们在构建日志里搜索 bin/java, 找到如下日志:

/Applications/Xcode.app/Contents/Developer/usr/bin/llvm-gcc -o /Users/yoncise/Downloads/openjdk/build/macosx-x86_64/bin/java  -m64 -Xlinker -rpath -Xlinker @loader_path/.  -Xlinker -install_name -Xlinker @rpath/java -L/Users/yoncise/Downloads/openjdk/build/macosx-x86_64/lib  -Wl,-all_load /Users/yoncise/Downloads/openjdk/build/macosx-x86_64/tmp/java/jli/obj64/static/libjli.a -framework Cocoa -framework Security -framework ApplicationServices -sectcreate __TEXT __info_plist ../../../../src/macosx/lib/Info-cmdline.plist \
	/Users/yoncise/Downloads/openjdk/build/macosx-x86_64/tmp/java/java/obj64/main.o -pthread -lz  -pthread 

我们可以发现 bin/java 的构建依赖 /Users/yoncise/Downloads/openjdk/build/macosx-x86_64/tmp/java/java/obj64/main.o.

根据目录我们可以知道, 这个文件是第一次完整构建的时候生成的并不是程序的源码, 删掉它, 再重复上面的步骤构建一遍.

再一次构建完成后, 我们在构建日志里搜索 obj64/main.o 可以找到如下日志:

/Applications/Xcode.app/Contents/Developer/usr/bin/llvm-gcc  -Os   -fno-strict-aliasing -fPIC -W -Wall  -Wno-unused -Wno-parentheses -pipe -m64 -fno-omit-frame-pointer -D_LITTLE_ENDIAN -F/System/Library/Frameworks/JavaVM.framework/Frameworks -F/System/Library/Frameworks/ApplicationServices.framework/Frameworks -I/usr/X11R6/include -D_DARWIN_UNLIMITED_SELECT  -DNDEBUG -DARCH='"x86_64"' -Dx86_64 -D_ALLBSD_SOURCE -DRELEASE='"1.7.0-internal"' -DFULL_VERSION='"1.7.0-internal-yoncise_2017_06_13_23_02-b00"' -DJDK_MAJOR_VERSION='"1"' -DJDK_MINOR_VERSION='"7"' -D_LARGEFILE64_SOURCE -D_GNU_SOURCE -D_REENTRANT -DMACOSX -D_LP64=1 -I/usr/X11R6/include -D_DARWIN_UNLIMITED_SELECT -I. -I/Users/yoncise/Downloads/openjdk/build/macosx-x86_64/tmp/java/java/CClassHeaders -I../../../../src/macosx/javavm/export -I../../../../src/share/javavm/export -I../../../../src/share/bin -I../../../../src/macosx/bin -I../../../../src/solaris/bin -DPACKAGE_PATH='"/opt/local"' -DPROGNAME='"java"' -DEXPAND_CLASSPATH_WILDCARDS -DLAUNCHER_NAME='"openjdk"'    -c -o /Users/yoncise/Downloads/openjdk/build/macosx-x86_64/tmp/java/java/obj64/main.o \
		-DRELEASE='"1.7.0-internal"' -DFULL_VERSION='"1.7.0-internal-yoncise_2017_06_13_23_02-b00"' -DJDK_MAJOR_VERSION='"1"' -DJDK_MINOR_VERSION='"7"' ../../../../src/share/bin/main.c

至此, 我们发现 bin/java 的入口文件是 src/share/bin/main.c !