matplotlib

fft

import numpy as np
import matplotlib.pyplot as plt
# matplotlib 绘图窗口中文显示问题
# from pylab import mpl
plt.rcParams['font.sans-serif'] = ['SimHei']  # 步骤一（替换sans-serif字体）
plt.rcParams['axes.unicode_minus'] = False   # 步骤二（解决坐标轴负数的负号显示问题）

# -------------------------
# 采样频率
Fs = 50
N = 1024
# -------------------------
# 数据--->采集的用来计算频率的数据
# Note:此处对N/Fs做一下解释：
# 一共采样了N个点，每个点的时间为1/Fs
# 因此此处采用N*1/Fs的方式
t = np.linspace(0, N/Fs, N)
# 直流
A0 = 5
# 波形1
A1 = 4
f1 = 1  # 1Hz
# 波形2
A2 = 2
f2 = 3  # 3Hz
# x = A1*np.sin(2*np.pi*f1*t) 		#用标准sin测试一下
# x = A0 + A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t)
x = A0 + A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi *
                                             f2*t) + (-5)*np.random.random(N)

# --------------------------
# 进行fft变换
y = np.fft.fft(x)
freq = np.fft.fftfreq(N, d=(t[0]-t[1]))

# -----------------------------------
# 能不能将高频率的波形滤掉
# 直接删除可好？----这个...太简洁了吧？
# ------------------------------------
# 定义一个滤波函数

def filter(freq, ft, flag=0, threshold=5):
    for i in range(N):
        if flag == 1:
            # 滤掉直流
            if np.abs(freq[i]) == 0:
                ft[i] = 0
        elif flag == 2:
            # 滤掉低频部分
            if np.abs(freq[i]) < threshold and np.abs(freq[i]) != 0:
                ft[i] = 0
        elif flag == 3:
            # 滤掉高频部分
            if np.abs(freq[i]) > threshold:
                ft[i] = 0
        elif flag == 4:
            # 尝试滤掉低振幅的杂波
            # 效果不是很好？是不是方法太low？
            if np.abs(ft[i])/N < 0.2:
                ft[i] = 0
        elif flag == 5:
            # 若已经知道当前最大频率
            # 则再高的频率直接滤掉
            # 相对于方法4而言，效果要好一些
            if np.abs(freq[i]) > 3.5:
                ft[i] = 0
        else:
            print("no filter")

figure, ax = plt.subplots()
figure.suptitle('频谱&滤波')
# --------------------------
# 1.0 显示原始图形
plt.subplot(221)
plt.axis([0, 5, -10, 10])
plt.plot(t, x, 'b')

# 1.1 频谱
plt.subplot(222)
plt.axis([-5, 5, 0, 3])
y1 = np.abs(y)/N  # 显示绝对值
plt.plot(freq, y1, 'r')

# 滤波
filter(freq, y, 5, threshold=2)
# 进行ifft反变换
x2 = np.fft.ifft(y)

# 1.2 频谱
plt.subplot(224)
plt.axis([-5, 5, 0, 3])
y1 = np.abs(y)/N  # 显示绝对值
plt.plot(freq, y1, 'r')

# 2.1fft逆变换之后
plt.subplot(223)
plt.axis([0, 5, -10, 10])
plt.plot(t, x2.real, 'b')  # 直接x2这里有warning，是不是应该x2.real
# 因为x2为Complex

plt.show()

fft2

# 安装：pip3.6 install opencv-pytho
import cv2
import numpy as np
import matplotlib.pyplot as plt

# ------------------------------
# 1.0读入图片
plt.figure('a beautiful girl')
beauty = cv2.imread('beauty.jpg', 0)
plt.subplot(221)
plt.imshow(beauty, 'gray')

# ------------------------------
# 1.1 fft2变换
y = np.fft.fft2(beauty)
y = np.fft.fftshift(y)
y1 = np.log(np.abs(y))
plt.subplot(222)
plt.imshow(y1, 'gray')

# -------------------------------
# 2.0 插入滤波器
# 边上是1，是高通
# 中心是1，是低通
row, col = beauty.shape
mask = np.zeros(beauty.shape, np.uint8)
radius = 20
# 这里使用//是因为整数除法，若直接/是小数除法
mask[row//2-radius:row//2+radius, col//2-radius:col//2+radius] = 1
y2 = y*mask
y2 = np.fft.ifftshift(y2)  # ...?

# -------------------------------
# 2.1 ifft2逆变换
x1 = np.fft.ifft2(y2)
plt.subplot(223)
x1 = np.abs(x1)
plt.imshow(x1, 'gray')

# 2.2 mask图像
plt.subplot(224)
mask2 = np.abs(mask)
plt.imshow(mask2, 'gray')

plt.show()

mplot3d

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 假设圆锥底面积的半径为R
R = 10
theta = np.arange(0, 2 * np.pi, 0.01)

X = R * np.cos(theta)
Y = R * np.sin(theta)
X, Y = np.meshgrid(X, Y)
Z = np.sqrt(X**2 + Y**2)

fig = plt.figure()
ax = Axes3D(fig)
# 用取样点(x,y,z)去构建曲面
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.cm.coolwarm)

# 给三个坐标轴注明
ax.set_xlabel('x label', color='r')
ax.set_ylabel('y label', color='g')
ax.set_zlabel('z label', color='b')

plt.show()

动态

# -*- coding: utf_8 -*-
# 注意，使用utf-8的文件，添加中文。必须包含上述
# 文件编码类型声明
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

# animation
fig = plt.figure()
ax = plt.axes(xlim=(0, 2), ylim=(-2, 2))
# create simply plot an empty line
# 注意加逗号的作用，表示其为tuple的对象
line, = ax.plot([], [], lw=2)
line2, = ax.plot([], [], lw=2)

def init():
    line.set_data([], [])
    line2.set_data([], [])
    return line, line2,

def animate(i):
    x = np.linspace(0, 2, 1000)
    y = np.sin(2 * np.pi * (x-0.01*i))
    y2 = np.cos(2 * np.pi * (x-0.01*i))
    line.set_data(x, y)
    line2.set_data(x, y2)
    return line, line2,

anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=200, interval=20, blit=True)

# 若是生成im.mp4似乎需要ffmpeg支持。
# 此处就先保存为*.htm/*.html看下效果
anim.save('im.htm', metadata={'artist': 'Guido'})
plt.show()

三维

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

figure = plt.figure()
# axes = Axes3D(figure)

R = 10
t = np.linspace(-np.pi, np.pi, 100)
X = R * np.cos(t)
Y = R * np.sin(t)
# X, Y = np.meshgrid(X, Y)
# Z = R**2 - X**2 - Y**2
X = R * np.outer(np.cos(t), np.sin(t))
Y = R * np.outer(np.sin(t), np.sin(t))
Z = R * np.outer(np.ones(np.size(t)), np.cos(t))
"""
    不明白为什么这样绘制出来的是球形
"""
axes = figure.add_subplot(121, projection='3d')
axes.plot_surface(X, Y, Z, rstride=4, cstride=4, color='b')
axes = figure.add_subplot(122, projection='3d')
axes.plot_wireframe(X, Y, Z)
plt.show()

螺旋

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 解决中文乱码问题，方便显示
plt.rcParams['font.sans-serif'] = ['Fangsong']
# 解决保存图像是负号'-'显示为方块的问题,或者转换负号为字符串
plt.rcParams['axes.unicode_minus'] = False

R = 10
b = 1
theta = np.arange(0, 4 * np.pi, 0.01)

X = R * np.cos(theta)
Y = R * np.sin(theta)
Z = b * theta


def plane_view():
    # 主要是为了查看一下
    # 三视图
    plt.axis('equal')
    plt.plot(X, Y, 'r', label='俯视图')
    plt.plot(X, Z, 'g', label='正视图')
    plt.plot(Y, Z, 'b', label='侧视图')
    plt.legend()    # 显示图例
    # plt.show()


def view_3d():
    # 显示实际的三维曲线
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    # 用取样点(x,y,z)去构建曲面
    ax.plot(X, Y, Z, c='r', label='螺旋线')

    # 给三个坐标轴注明
    ax.set_xlabel('x label', color='r')
    ax.set_ylabel('y label', color='g')
    ax.set_zlabel('z label', color='b')
    plt.legend()    # 显示图例
    # plt.show()


if __name__ == '__main__':
    plane_view()
    view_3d()
    plt.show()