关于fft，svd，卷积，图像卷积

1.内容

fft，svd

2.代码

2.1. Ensumble.py

#!/usr/bin/python
# -*- coding:utf-8 -*-

import operator


def c(n, k):
    return reduce(operator.mul, range(n-k+1, n+1)) / reduce(operator.mul, range(1, k+1))


def bagging(n, p):
    s = 0
    for i in range(n / 2 + 1, n + 1):
        s += c(n, i) * p ** i * (1 - p) ** (n - i)
    return s


if __name__ == "__main__":
    for t in range(10, 101, 10):
        print t, '次采样正确率：', bagging(t, 0.6)

2.2 fft

# !/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

# 三角波函数
def triangle_wave(size, T):
    t = np.linspace(-1, 1, size, endpoint=False) #最后一点设置为开，因为周期性
    # where
    # y = np.where(t < 0, -t, 0)
    # y = np.where(t >= 0, t, y)
    y = np.abs(t)
    y = np.tile(y, T) - 0.5
    x = np.linspace(0, 2*np.pi*T, size*T, endpoint=False)
    return x, y


def sawtooth_wave(size, T):
    t = np.linspace(-1, 1, size)
    y = np.tile(t, T)
    x = np.linspace(0, 2*np.pi*T, size*T, endpoint=False)
    return x, y


def triangle_wave2(size, T):
    x, y = sawtooth_wave(size, T)
    return x, np.abs(y)


def non_zero(f):
    f1 = np.real(f)
    f2 = np.imag(f)
    eps = 1e-4
    return f1[(f1 > eps) | (f1 < -eps)], f2[(f2 > eps) | (f2 < -eps)]


if __name__ == "__main__":
    mpl.rcParams['font.sans-serif'] = [u'simHei'] # 使用中文
    mpl.rcParams['axes.unicode_minus'] = False
    np.set_printoptions(suppress=True) #显示小数

    x = np.linspace(0, 2*np.pi, 16, endpoint=False)  #不要终点
    print '时域采样值：', x
    y = np.sin(2*x) + np.sin(3*x + np.pi/4)
    # y = np.sin(x)

    N = len(x)
    print '采样点个数：', N
    print '\n原始信号：', y
    f = np.fft.fft(y)
    print '\n频域信号：', f/N
    a = np.abs(f/N)
    print '\n频率强度：', a

    iy = np.fft.ifft(f) #做逆fft
    print '\n逆傅里叶变换恢复信号：', iy
    print '\n虚部：', np.imag(iy)
    print '\n实部：', np.real(iy)
    print '\n恢复信号与原始信号是否相同：', np.allclose(np.real(iy), y) #判等函数

    plt.subplot(211)
    plt.plot(x, y, 'go-', lw=2) 
    plt.title(u'时域信号', fontsize=15)
    plt.grid(True)
    plt.subplot(212)
    w = np.arange(N) * 2*np.pi / N
    print u'频率采样值：', w
    plt.stem(w, a, linefmt='r-', markerfmt='ro')
    plt.title(u'频域信号', fontsize=15)
    plt.grid(True)
    plt.show()

    # 三角/锯齿波
    x, y = triangle_wave(20, 5)
    # x, y = sawtooth_wave(20, 5)
    N = len(y)
    f = np.fft.fft(y)
    # print '原始频域信号：', np.real(f), np.imag(f)
    print '原始频域信号：', non_zero(f)
    a = np.abs(f / N)

    # np.real_if_close
    f_real = np.real(f)
    eps = 0.1 * f_real.max()
    print eps
    f_real[(f_real < eps) & (f_real > -eps)] = 0
    f_imag = np.imag(f)
    eps = 0.1 * f_imag.max()
    print eps
    f_imag[(f_imag < eps) & (f_imag > -eps)] = 0
    f1 = f_real + f_imag * 1j
    y1 = np.fft.ifft(f1)
    y1 = np.real(y1)
    # print '恢复频域信号：', np.real(f1), np.imag(f1)
    print '恢复频域信号：', non_zero(f1)

    plt.figure(figsize=(8, 8), facecolor='w')
    plt.subplot(311)
    plt.plot(x, y, 'g-', lw=2)
    plt.title(u'三角波', fontsize=15)
    plt.grid(True)
    plt.subplot(312)
    w = np.arange(N) * 2*np.pi / N
    plt.stem(w, a, linefmt='r-', markerfmt='ro')
    plt.title(u'频域信号', fontsize=15)
    plt.grid(True)
    plt.subplot(313)
    plt.plot(x, y1, 'b-', lw=2, markersize=4)
    plt.title(u'三角波恢复信号', fontsize=15)
    plt.grid(True)
    plt.tight_layout(1.5, rect=[0, 0.04, 1, 0.96])
    plt.suptitle(u'快速傅里叶变换FFT与频域滤波', fontsize=17)
    plt.show()

2.3 svd

#!/usr/bin/python
#  -*- coding:utf-8 -*-

import numpy as np
import os
from PIL import Image
import matplotlib.pyplot as plt
import matplotlib as mpl
from pprint import pprint


def restore1(sigma, u, v, K):  # 奇异值、左特征向量、右特征向量
    m = len(u)
    n = len(v[0])
    a = np.zeros((m, n))
    for k in range(K):
        uk = u[:, k].reshape(m, 1)
        vk = v[k].reshape(1, n)
        a += sigma[k] * np.dot(uk, vk)
    a[a < 0] = 0
    a[a > 255] = 255
    # a = a.clip(0, 255)
    return np.rint(a).astype('uint8')


def restore2(sigma, u, v, K):  # 奇异值、左特征向量、右特征向量
    m = len(u)
    n = len(v[0])
    a = np.zeros((m, n))
    for k in range(K+1):
        for i in range(m):
            a[i] += sigma[k] * u[i][k] * v[k]
    a[a < 0] = 0
    a[a > 255] = 255
    return np.rint(a).astype('uint8')


if __name__ == "__main__":
    A = Image.open("6.son.png", 'r')
    output_path = r'.\Pic'
    if not os.path.exists(output_path):
        os.mkdir(output_path)
    a = np.array(A)
    K = 50
    u_r, sigma_r, v_r = np.linalg.svd(a[:, :, 0])
    u_g, sigma_g, v_g = np.linalg.svd(a[:, :, 1])
    u_b, sigma_b, v_b = np.linalg.svd(a[:, :, 2])
    plt.figure(figsize=(10,10), facecolor='w')
    mpl.rcParams['font.sans-serif'] = [u'simHei']
    mpl.rcParams['axes.unicode_minus'] = False
    for k in range(1, K+1):
        print k
        R = restore1(sigma_r, u_r, v_r, k)
        G = restore1(sigma_g, u_g, v_g, k)
        B = restore1(sigma_b, u_b, v_b, k)
        I = np.stack((R, G, B), 2)
        Image.fromarray(I).save('%s\\svd_%d.png' % (output_path, k))
        if k <= 12:
            plt.subplot(3, 4, k)
            plt.imshow(I)
            plt.axis('off')
            plt.title(u'奇异值个数：%d' % k)
    plt.suptitle(u'SVD与图像分解', fontsize=18)
    plt.tight_layout(2)
    plt.subplots_adjust(top=0.9)
    plt.show()

股票数据
2.4 convolve.py

# !/usr/bin/python
# -*- coding:utf-8 -*-

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt


if __name__ == "__main__":
    stock_max, stock_min, stock_close, stock_amount = np.loadtxt('6.SH600000.txt', delimiter='\t', skiprows=2, usecols=(2, 3, 4, 5), unpack=True) #unpack代表返回一个还是多个
    N = 100
    stock_close = stock_close[:N]
    print stock_close

	# 滑动平均，每n个值取平均得到新的值
    n = 5
    weight = np.ones(n)
    weight /= weight.sum()
    print weight
    stock_sma = np.convolve(stock_close, weight, mode='valid')  # simple moving average

    weight = np.linspace(1, 0, n)
    weight = np.exp(weight)
    weight /= weight.sum()
    print weight
    stock_ema = np.convolve(stock_close, weight, mode='valid')  # exponential moving average

    t = np.arange(n-1, N)
    poly = np.polyfit(t, stock_ema, 10)
    print poly
    stock_ema_hat = np.polyval(poly, t)

    mpl.rcParams['font.sans-serif'] = [u'SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.plot(np.arange(N), stock_close, 'ro-', linewidth=2, label=u'原始收盘价')
    t = np.arange(n-1, N)
    plt.plot(t, stock_sma, 'b-', linewidth=2, label=u'简单移动平均线')
    plt.plot(t, stock_ema, 'g-', linewidth=2, label=u'指数移动平均线')
    plt.legend(loc='upper right')
    plt.grid(True)
    plt.show()

    plt.figure(figsize=(9, 6))
    plt.plot(np.arange(N), stock_close, 'r-', linewidth=1, label=u'原始收盘价')
    plt.plot(t, stock_ema, 'g-', linewidth=2, label=u'指数移动平均线')
    plt.plot(t, stock_ema_hat, 'm-', linewidth=3, label=u'指数移动平均线估计')
    plt.legend(loc='upper right')
    plt.grid(True)
    plt.show()

图像卷积
2.5 image_convolve.py

#!/usr/bin/python
#  -*- coding:utf-8 -*-

import numpy as np
import os
from PIL import Image


def convolve(image, weight):
    height, width = image.shape
    h, w = weight.shape
    height_new = height - h + 1
    width_new = width - w + 1
    image_new = np.zeros((height_new, width_new), dtype=np.float)
    for i in range(height_new): #提取出子图，然后跟权值进行按位的乘
        for j in range(width_new):
            image_new[i,j] = np.sum(image[i:i+h, j:j+w] * weight)
    image_new = image_new.clip(0, 255)
    image_new = np.rint(image_new).astype('uint8')
    return image_new

# image_new = 255 * (image_new - image_new.min()) / (image_new.max() - image_new.min())

if __name__ == "__main__":
    A = Image.open("6.son.png", 'r')
    output_path = '.\\Pic\\'
    if not os.path.exists(output_path):
        os.mkdir(output_path)
    a = np.array(A)
    soble_x = np.array(([-1, 0, 1], [-2, 0, 2], [-1, 0, 1])) #给出各种算子
    soble_y = np.array(([-1, -2, -1], [0, 0, 0], [1, 2, 1]))
    soble = np.array(([-1, -1, 0], [-1, 0, 1], [0, 1, 1]))
    prewitt_x = np.array(([-1, 0, 1], [-1, 0, 1], [-1, 0, 1]))
    prewitt_y = np.array(([-1, -1,-1], [0, 0, 0], [1, 1, 1]))
    prewitt = np.array(([-2, -1, 0], [-1, 0, 1], [0, 1, 2]))
    laplacian = np.array(([0, -1, 0], [-1, 4, -1], [0, -1, 0]))
    laplacian2 = np.array(([-1, -1, -1], [-1, 8, -1], [-1, -1, -1]))
    weight_list = ('soble_x', 'soble_y', 'soble', 'prewitt_x', 'prewitt_y', 'prewitt', 'laplacian', 'laplacian2')
    print '梯度检测：'
    for weight in weight_list:
        print weight, 'R',
        R = convolve(a[:, :, 0], eval(weight)) # 对不同的通道分别做处理，权值进行遍历
        print 'G',
        G = convolve(a[:, :, 1], eval(weight))
        print 'B'
        B = convolve(a[:, :, 2], eval(weight))
        I = 255 - np.stack((R, G, B), 2) #将每三个通道合成一个
        Image.fromarray(I).save(output_path + weight + '.png') #储存png文件

    # # X & Y
    # print '梯度检测XY：'
    # for w in (0, 2):
    #     weight = weight_list[w]
    #     print weight, 'R',
    #     R = convolve(a[:, :, 0], eval(weight))
    #     print 'G',
    #     G = convolve(a[:, :, 1], eval(weight))
    #     print 'B'
    #     B = convolve(a[:, :, 2], eval(weight))
    #     I1 = np.stack((R, G, B), 2)
    #
    #     weight = weight_list[w+1]
    #     print weight, 'R',
    #     R = convolve(a[:, :, 0], eval(weight))
    #     print 'G',
    #     G = convolve(a[:, :, 1], eval(weight))
    #     print 'B'
    #     B = convolve(a[:, :, 2], eval(weight))
    #     I2 = np.stack((R, G, B), 2)
    #
    #     I = 255 - np.maximum(I1, I2)
    #     Image.fromarray(I).save(output_path + weight[:-2] + '.png')