傅里叶变换

视频地址：www.bilibili.com/video/BV1T3…

前言

不知道大家小时候有没有玩过万花尺。

其实我们用傅里叶变换就可以模拟出万花尺的效果。

最近为了研究傅里叶变换，我看了许多相关的书和视频。

书都是教科书，学校里的那种；

视频是B站里的“国家精品课”— 数字图像处理。

这些看着都很好，很有意义，但就是知其然，不知其所以然。

其中有很多知识点一开始看得明明白白，逐步深入后，就不知不觉的迷失自我了。

仔细想来，导致这个结果的，大致有两点：

首先，是我自身基础有些薄弱，而且也不够才思敏捷，不是一个学霸；

其次，应该是老师和他的教科书没有把授人以渔的事情做好。

后来，我发现了3Blue1Brown 的课程，瞬间惊为天人，他用几何图形把复杂的数学知识给可视化了，让人真正的理解了各种公式和定律的前因后果。

接下来，我会参考其《形象展示傅里叶变换》，用js为大家完成从傅里叶变换到代码实现的最后一公里配送。

1-傅里叶变换的概念

傅里叶变换(Fourier transform）是一种线性积分变换，用于信号在时域和频域之间的变换，这在物理学和工程学中用得很多。

大家可以暂且将时域理解为我们当前所感知的真实世界，而频域就是真实世界的另一种表现形式。

傅里叶变换可以对自然界中的各种信号进行分析，确定其基本成分，比如分解声音信号：

接下来咱们就解释一下声音的分解原理。

2-声音的分解原理

我们通过一个单频信号来解释声音的分解原理。

1.绘制时域图，用余弦曲线建立一个单频信号，信号时长为4.5秒，信号频率为3拍/秒。

2.绘制信号缠绕图，将信号缠绕在一个圆上。

信号在缠绕的时候，是有一个缠绕频率的。

下图是按照0.3圈/s 的缠绕频率绘制的。

我们可以对照着时域图看一下：

在时域图中，越高的点，离信号缠绕图的圆心就越远。

当前这张信号缠绕图是受两种不同的频率影响的：

• 信号频率：3拍/s
• 信号缠绕频率：0.3圈/s

信号缠绕频率是可以自由变换的。

当信号缠绕频率和信号频率一样时，即信号缠绕频率是3圈/s，就会有神奇的事情发生：

由上图可见，时域图中所有高处的点都在信号缠绕图的x轴右侧，低处的点在x轴左侧。

这个规律非常重要，我们可以利用它来做一些事情。

我们可以把信号缠绕图里的信号线视之为一个有质量的东西，比如一朵花。

在这朵花里面必然有一个质量的中心点，也就是质心。

当我们改变缠绕频率的时候，图像的缠绕方式会变化，质心也会变化。

对于大部分缠绕频率来说，信号图像的峰、谷会均匀的分布在信号缠绕图的原点周围，从使得质心一直在原点附近。

然而，当信号的缠绕频率等于信号频率时，也就是3圈/s，质心会偏向右侧。

为了更形象的展示这个现象，我们可以画一个频域图，用来就记录信号缠绕图中，随着缠绕频率的改变，质心到原点的距离变化。

从上图中，我们可以明显看到，在缠绕频率为3的地方有一个局部极大值。

在此还要注意一下缠绕频率为0的地方，这个地方之所以有一个距离为1的极大值，是因为时域图的波动区间是[0,2]。

缠绕频率为0时可以理解为把时域里的信号曲线沿x轴挤压成一捆长度为2线，其质心就在2/2=1的位置。

如果我们将时域图的波动区间设置为[-1,1]，那么当绕频率为0时，质心就在(-1+1)/2=0 的位置。

由此可见，我们想要分解频率，就需要重点关注缠绕频率在3处的点。

生成上面这张频域图的变换方式，就叫做时域的近傅里叶变换。

近傅里叶变换和真正的傅里叶变换还是有一点差距的，这个我们稍后说。

3.以同样的原理，再换一个2拍/s 的单频率的信号试试。

通过上面的实验，我们可以发现一个重要规律：

在频域图中，当单频信的缠绕频率与信号频率一致时，会有一个局部极大值，而在其它地方都接近于0。

4.以此原理，我们还可以绘制双频信号的频域图。

我们把缠绕频率2和3记下来，换成信号频率，便可以还原出原本的两个单频信号了。

到目前为止，我们就已经知道了如何分解声音信号，其关键就是如何在信号缠绕图里获取质心，以及质心到原点的距离。

接下来，我们再说一下如何把时域图缠绕到圆上。

3-把时域图缠绕到圆上

把时域图缠绕到圆上的关键在于将欧拉公式与时域图相乘。

欧拉公式是复分析领域的公式，它将三角函数与复指数函数关联到一起。

• e ：自然对数的底数
• i ：虚数单位，即-1的2次方根
• cos φ ：φ的余弦值
• sin φ ：φ的正弦值
• φ ：弧度

对于欧拉公式的来龙去脉，咱们先不细说。

大家可以将上面的欧拉公式理解为一个单位向量(cos φ,sin φ)。

这个单位向量是在基于弧度φ做旋转的：

• φ递增时，单位向量逆时针旋转；
• φ递减时，单位向量顺时针旋转。

为了让向量以时间为单位做旋转，我们可以把欧拉公式写成这样：

• 2π：一个整圆的弧度
• t：时间

如果t是以秒为单位的，那么上面的公式代表的就是每t秒旋转一圈。

若是觉得上面的公式会转得太快或太慢，我们可以再加一个频率f，用来控制求旋转速度：

当f 代表的是时域图的缠绕频率时，上式就可以理解为按照f 圈/s 的频率旋转。

现在我们已经在渐渐接近傅里叶变换了。

在傅里叶变换中有一个习惯，通常认为旋转是顺时针的。

所以我们需要在e的指数前面在加个负号：

这样，时间递增时，欧拉公式里的单位向量会顺时针旋转。

我们把时域图乘以上面的欧拉公式，便可以画出信号曲线绕圆形旋转的效果。

算法如下：

• g(t) ：时域函数

现在我们已经知道了如何把信号曲线缠绕到圆上，接下来我们就可以去寻找这一团曲线的质心了。

4-寻找质心

1.在信号缠绕图上采样，样采样本点越多，计算结果越精确。

2.对采样点做一个向量加法的运算，然后将运算结果除以样本数，就可以得到质心位置了。

公式如下：

• N ：采样数量
• tk ：当前采样点的采样时间

我们也可以将其写成积分的形式：

• t1 ：采样开始时间
• t2 ：采样结束时间
• t2-t1：时间长度
• dt ：采样的时间间隔，取极小值

在上式中，对复函数做定积分运算是有点奇怪的，所以我们可以先重点理解其描述的意思即可。

现在上面的积分公式和真正的傅里叶变换就只有一处不同了，将除以时间长度的步骤删掉即可：

用傅里叶变换计算出来的质心位置和之前的公式并没有本质的区别，它们只是一个时间长度的缩放差异。

注：当时域的函数为g(t) 时，其对应的频域函数就会在g(t)上面加一个^。

有了质心之后，我们接下来就可以根据质心到原点的距离绘制频域图了。

总结一下我们到目前为止所说的几个要点：

• 时域图是以时间为自变量进行绘图的。

  

• 频域图是以信号曲线的缠绕频率为自变量进行绘图的。

现在傅里叶变换的基本原理已经搞定，接下来我们就用代码实现一下。

5-代码实现

1.封装傅里叶对象Fourier，这个对象就是用于傅里叶变换的。

const PI2 = Math.PI * 2

const defAttr = () => ({
    // 时间长度
    duration: 4.5,
    // 缠绕频率的极大值
    maxTwining: 4.5,
    // 时间采样距离
    dTime: 0.01,
    // 缠绕频率采样距离
    dTwining: 0.01,
    // 信号波动幅度
    amplitude: 1,
    // 信号偏移
    offsetY: 0,
    // 信号频率集合
    signalFrequencies: [3],
    // 信号缠绕频率，圈/s
    turns: 3,
    // 信号生成器
    signal: () => {},
})
export default class Fourier {
    constructor(attr = {}) {
        Object.assign(this, defAttr(), attr)
        this.update()
    }
    // 更新
    update() {
        //信号生成器
        this.signal = this.SignalGenerator()
    }
    // 时域图的顶点集合
    timePoints(fn = (p) => p) {
        const { duration, dTime } = this
        const points = []
        for (let t = 0; t < duration; t += dTime) {
            const p = [t, this.signal(t)]
            points.push(fn(p))
        }
        return points
    }
    // 信号缠绕图的顶点集合
    eulerPoints(fn = (p) => p) {
        const { duration, dTime, turns } = this
        const points = []
        for (let t = 0; t < duration; t += dTime) {
            const p = this.twining(t, turns)
            points.push(fn(p))
        }
        return points
    }
    // 频域图的顶点集合
    frequencyPoints(fn = (p) => p) {
        const { duration, maxTwining, dTime, dTwining } = this
        // 在时域图上的采样数量
        const sampleNum = this.duration / this.dTime
        // 顶点集合
        const points = []
        // 遍历缠绕频率
        for (let f = 0; f < maxTwining; f += dTwining) {
            //质心位置
            let [cx, cy] = [0, 0]
            // 遍历时间
            for (let t = 0; t < duration; t += dTime) {
                const [x, y] = this.twining(t, f)
                cx += x
                cy += y
            }
            // 质心距离
            const length = Math.sqrt(cx * cx + cy * cy) / sampleNum
            points.push(fn([f, length]))
        }
        return points
    }
    // 把时域图缠绕在圆上
    twining(time, twiningFrequency) {
        const ang = -PI2 * time * twiningFrequency
        const cos = Math.cos(ang)
        const sin = Math.sin(ang)
        const y = this.signal(time)
        return [cos * y, sin * y]
    }
    // 信号生成器的生成函数
    SignalGenerator() {
        const { signalFrequencies, amplitude, offsetY } = this
        // 信号曲线在1s内的总弧度
        const rads = signalFrequencies.map((frequency) => PI2 * frequency)
        // 信号的频道数量
        const len = signalFrequencies.length
        // 返回信号生成器
        return function (t) {
            let y = 0
            for (let i = 0; i < len; i++) {
                y += amplitude * Math.cos(t * rads[i]) + offsetY
            }
            return y / len
        }
    }
}

在我上面的代码里都有详细的注释，所以我就不再做过多赘述了。若大家哪里不懂，可以微信沟通1051904257。

2.使用echarts，把信号从时域转频域。

<!DOCTYPE html>
<html lang="en">
    <head>
        <meta charset="UTF-8" />
        <title>傅里叶变换</title>
        <style>
            html {
                height: 100%;
            }
            body {
                margin: 0;
                overflow: hidden;
                height: 100%;
            }
            #main {
                height: 100%;
            }
        </style>
    </head>
    <body>
        <div id="main"></div>
        <script src="https://lib.baomitu.com/echarts/5.3.3/echarts.min.js"></script>
        <script type="module">
            import Fourier from './Fourier.js'

            /* 实例化echarts */
            const main = document.getElementById('main')
            const chart = echarts.init(main)

            /* 傅里叶对象 */
            // 时间长度
            const duration = 4.5
            // 缠绕频率的极大值
            const maxTwining = 4.5
            // 时间采样距离
            const dTime = 0.01
            // 缠绕频率采样距离
            const dTwining = 0.01
            // 信号波动幅度
            const amplitude = 1
            // 信号偏移
            const offsetY = 0
            // 信号频率集合
            const signalFrequencies = [3]
            // 当前信号缠绕频率，圈/s
            let turns = 0.3
            // 实例化傅里叶对象
            const fourier = new Fourier({
                duration,
                maxTwining,
                dTime,
                dTwining,
                amplitude,
                offsetY,
                signalFrequencies,
                turns,
            })

            /* 多绘图区 */
            let option = {
                xAxis: [
                    {
                        min: 0,
                        max: duration,
                        name: '时间(s)',
                        gridIndex: 0,
                    },
                    {
                        min: -amplitude - offsetY,
                        max: amplitude + offsetY,
                        name: 'x',
                        gridIndex: 1,
                    },
                    {
                        min: 0,
                        max: maxTwining,
                        name: '缠绕频率',
                        gridIndex: 2,
                    },
                ],
                yAxis: [
                    {
                        min: -amplitude + offsetY,
                        max: amplitude + offsetY,
                        name: '信号值',
                        gridIndex: 0,
                    },
                    {
                        min: -amplitude - offsetY,
                        max: amplitude + offsetY,
                        name: 'y',
                        gridIndex: 1,
                    },
                    {
                        min: 0,
                        max: fourier.amplitude + offsetY,
                        name: '质心距离',
                        gridIndex: 2,
                    },
                ],
                grid: [
                    { bottom: '73%' },
                    {
                        top: '36%',
                        bottom: '36%',
                        width: main.clientHeight * 0.28,
                    },
                    { top: '73%' },
                ],
            }
            chart.setOption(option)

            /* 时域图 */
            const time = {
                type: 'polyline',
                shape: {
                    points: fourier.timePoints((p) => toPixel(0, p)),
                },
            }

            /* 信号缠绕图 */
            const euler = {
                id: 2,
                type: 'polyline',
                shape: {
                    points: fourier.eulerPoints((p) => toPixel(1, p)),
                },
            }

            /* 频域图 */
            const frequency = {
                type: 'polyline',
                shape: {
                    points: fourier.frequencyPoints((p) => toPixel(2, p)),
                },
            }

            /* 频域图标记线 */
            // 鼠标划入样式
            const markLineOver = {
                id: 4,
                style: {
                    lineWidth: 2,
                    shadowColor: '#68e1ff',
                    shadowOffsetY: 3,
                    shadowBlur: 3,
                },
            }
            // 鼠标划出样式
            const markLineOut = {
                id: 4,
                style: {
                    lineWidth: 1,
                },
            }
            // 标记线
            const markLine = {
                id: 4,
                type: 'polyline',
                shape: {
                    points: markLinePoints(),
                },
                style: {
                    stroke: '#00acec',
                },
                draggable: 'horizontal',
                onmouseover: () => {
                    chart.setOption({
                        graphic: [markLineOver],
                    })
                },
                ondrag: ({ offsetX }) => {
                    // 像素坐标位转图表坐标位
                    const turns = chart
                        .convertFromPixel(
                            {
                                gridIndex: 2,
                            },
                            [offsetX, 0]
                        )[0]
                        .toFixed(3)
                    // 吸附整数
                    const roundTurns = Math.round(turns)
                    fourier.turns =
                        Math.abs(roundTurns - turns) < 0.02 ? roundTurns : turns

                    // 更新信号缠绕图
                    chart.setOption({
                        graphic: [
                            markLineOver,
                            {
                                id: 5,
                                ...textPos(),
                                style: {
                                    ...textDefaultStyle,
                                    text: fourier.turns,
                                },
                            },
                            {
                                id: 2,
                                shape: {
                                    points: fourier.eulerPoints((p) => toPixel(1, p)),
                                },
                            },
                        ],
                    })
                },
                onmouseout: () => {
                    chart.setOption({
                        graphic: [markLineOut],
                    })
                },
            }
            /* 缠绕频率提示文字 */
            const textDefaultStyle = {
                text: fourier.turns,
                font: '16px cursive',
                textAlign: 'center',
                textVerticalAlign: 'bottom',
                fill: '#000',
                stroke: '#fff',
                lineWidth: 5,
            }
            const text = {
                id: 5,
                type: 'text',
                ...textPos(),
                style: {
                    ...textDefaultStyle,
                    text: fourier.turns,
                },
                z: 1000,
            }

            // 绘图
            chart.setOption({
                graphic: [time, euler, frequency, markLine, text],
            })

            // 图表坐标位转像素坐标位
            function toPixel(gridIndex, p) {
                return chart.convertToPixel({ gridIndex }, p)
            }
            // 标记线的顶点集合
            function markLinePoints() {
                const { turns, amplitude, offsetY } = fourier
                return [
                    toPixel(2, [turns, 0]),
                    toPixel(2, [turns, amplitude + offsetY]),
                ]
            }
            // 缠绕频率提示文字位置
            function textPos() {
                const { turns, amplitude, offsetY } = fourier
                const [x, y] = toPixel(2, [turns, amplitude + offsetY])
                return { x, y: y - 5 }
            }
        </script>
    </body>
</html>

其最终效果就是我们之前截图所看见的样子。

我们可以修改一下上面的参数，比如在信号频率集合里再加一个频率，修改一下缠绕圈数。

// 信号频率集合
const signalFrequencies = [3, 2]
// 信号缠绕频率，圈/s
let turns = 1.75

效果如下：

我们可以把傅里叶当成万花尺玩：

参数如下：

// 时间长度
const duration = 20
// 采样样距离
const dx = 0.01
// 信号波动幅度
const amplitude = 1
// 信号偏移
const offsetY = 1.5
// 信号频率集合
const signalFrequencies = [6, 4]
// 信号缠绕频率，圈/s
let turns = 1.0526

好啦，关于傅里叶变换，咱们就说到这，后面我们会再对傅里叶变换进行逐步深入，做出更多好玩的效果。