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

末蓝、 2024-04-17 19:47 80阅读 0赞

# 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')