代码演示傅里叶合成演示

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

很久前写过相关傅里叶学习笔记傅里叶入门–动手演示波形叠加，演示正余弦信号叠加合成矩形信号。当时对代码及资料理解浅薄，至于很多地方描述不正确，其中信号合成后和原先信号不一致，此篇笔记就是重新探讨这个问题。

---

一、涉及公式

资料来源：

《实用数字信号处理：从原理到应用》作者：StevenW.Smith 第8章

式一：计算时域信号 xi

式二：计算频域信号 Re和Im

式三：计算 ReX 上划线 ImX 上划线

特殊情况：

二、实践

1.环境

我的环境 python3.7。
需要 matplotlib , numpy 两个组件

如果没有组件先安装，比如

pip install matplotlib
1

2.完整代码

import math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.lines import Line2D

n = 100
N = 2*n  # 200点时域信号
kl = n + 1

pi = math.pi  #  数学pi 约= 3.14159

x = np.zeros((N, ))
# x[50:150] = 200 # 构造一个矩形信号 x   (这两行代码使用其中一行即可)
x[50:100] = 100; x[100:150] = 200   # 阶梯信号   (这两行代码使用其中一行即可)

Re = np.zeros((kl, ))
Im = np.zeros((kl, ))
Rel = np.zeros((kl, ))
Iml = np.zeros((kl, ))

# 计算频域 Re 和 Im
# 对应 式二
def calculateReIm():
    for k in range(kl):
        for i in range(N):
            Re[k] = Re[k] + x[i] * math.cos(2 * pi * k * i / N)
            Im[k] = Im[k] - x[i] * math.sin(2 * pi * k * i / N)

# 计算 Re上划线 和 Im上划线
# 对应 式三
def calculateReImSHX():
    global Rel
    global Iml

    Rel = Re / (N / 2)
    Iml = -Im / (N / 2)
    Rel[0] = Re[0] / N
    Iml[n] = Im[n] / N

# 计算时域 x
# 对应 式一
xi = np.zeros((N, ))
def calculateX():

    for k in range(kl):
        for i in range(N):
            xi[i] = xi[i] + Rel[k] * math.cos(2 * pi * k * i / N)
            xi[i] = xi[i] + Iml[k] * math.sin(2 * pi * k * i / N)


if __name__=='__main__':
    calculateReIm()
    calculateReImSHX()
    # calculateX()  # 此画图暂时用不到， 可放开作为验证 xi

'''
以上部分属于基本数值计算
==============================================
以下布局都属于画图
'''
xiline = np.zeros((N, ))
cos = []
sin = []
sumQx = np.zeros((N, ))
linex = []

temp = 0

linecos = Line2D([], [], color='b')
linesin = Line2D([], [], color='g')
linesumQx = Line2D([], [], color='r')
linesumX = np.arange(N)
def up(v):
    global cos, sin, sumQx, temp, linex
    k = v[0]
    i = v[1]
    if temp != k:
        cos = []
        sin = []
        linex = []
        temp = k
    c = Rel[k] * math.cos(2 * pi * k * i / N)
    s = Iml[k] * math.sin(2 * pi * k * i / N)
    cos.append(c)
    sin.append(s)
    sumQx[i] += c + s
    linex.append(i)

    linecos.set_data(linex, cos)
    linesin.set_data(linex, sin)
    linesumQx.set_data(linesumX, sumQx)
    time_text.set_text(time_template % (i, k))
    return linecos, linesin, linesumQx, time_text,


def emitter():
    for i in range(kl):
        for j in range(N):
            yield i, j


def init():
    linecos.set_data([], [])
    linesin.set_data([], [])
    linesumQx.set_data(linesumX, sumQx)
    time_text.set_text('i = 0 k = 0')
    return linecos, linesin, linesumQx, time_text,


fig, ax = plt.subplots()
ax.add_line(linecos)
ax.add_line(linesin)
ax.add_line(linesumQx)
ax.set_ylim(-250, 250)
ax.set_xlim(0, N)
ax.set_title('red: sum, blue: cos, green: sin')
time_template = 'i = %d k = %d'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
ani = animation.FuncAnimation(fig, up, emitter, init_func=init, interval=10, blit=True, repeat=False, save_count=10000)
plt.show()
# ani.save("fuliye-jx.mp4")  # 保存成 mp4
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

---

总结

保存MP4

如果要保存MP4, 先安装ffmpeg，然后放开 save 那行代码。 ffmpeg安装教程参考
结果演示：
矩形：

阶梯：