1. 程式人生 > >python模擬離散傅立葉變換

python模擬離散傅立葉變換

code:

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

# 產生三角波的函式
def triangle(size, T):
	# 生成-1到1之間的size個時間點
	t = np.linspace(-1, 1, size, endpoint=False)
	# 這裡使用y=|x|函式生成倒三角波一樣的影象
	y = np.abs(t)
	# 上面已經生成一個三角波影象,然後進行復制T個
	y = np.tile(y, T) - 0.5
	# 接著吧上面總共取樣的T個週期的所有采樣點集合到x變數
	x = np.linspace(0, 2 * np.pi * T, size * T, endpoint=False)
	return x, y

# 把取樣點的0點去掉
def delete_zero(f):
	f1 = np.real(f)
	f2 = np.imag(f)
	# 設定一個極小接近0的值
	e_min = 1e-5
	# 同上面的極小值比較去0
	return f1[(f1 > e_min) | (f1 < -e_min)], f2[(f2 > e_min) | (f2 < -e_min)]

if __name__ == "__main__":

    x = np.linspace(0, 2*np.pi, 32, endpoint=False)
    print('時域上訊號取樣值:\n', x)
    # y = np.sin(2*x) + np.sin(3*x + np.pi/4) + np.sin(5*x)
    y = np.sin(x)

    N = len(x)
    print('總共取樣點個數:\n', N)
    print('\n原始訊號值:\n', y)
    f = np.fft.fft(y)
    print('\n頻域上訊號取樣值:\n', f/N)
    a = np.abs(f/N)
    print('\n頻率取樣點強度:\n', a)

    iy = np.fft.ifft(f)
    print('\n逆傅立葉變換恢復訊號:\n', iy)
    print('\n虛部:\n', np.imag(iy))
    print('\n實部:\n', np.real(iy))
    print('\n恢復訊號與原始訊號是否相同:\n', np.allclose(np.real(iy), y))

    plt.subplot(211)
    plt.plot(x, y, 'go-', lw=2)
    plt.title('time domain signal', fontsize=15)
    plt.grid(True)
    plt.subplot(212)
    w = np.arange(N) * 2*np.pi / N
    print('頻率取樣值:\n', w)
    plt.stem(w, a, linefmt='r-', markerfmt='ro')
    plt.title('frequency domain signal', fontsize=15)
    plt.tight_layout(2)  
    plt.subplots_adjust(top=0.9)
    plt.grid(True)
    plt.show()

    # 模擬三角鋸齒波
    x, y = triangle(30, 7)
    N = len(y)
    f = np.fft.fft(y)
    print("原始的三角頻域訊號:\n", np.real(f), np.imag(f))
    print("原始去0後的訊號:\n", delete_zero(f))
    a = np.abs(f/N)
    
    f_real = np.real(f)
    e_min = 0.3 * f_real.max()
    print("e_min:\n", e_min)
    # 去0處理
    f_real[(f_real < e_min) & (f_real > -e_min)] = 0
    
    f_imag = np.imag(f)
    ee_min = 0.3 * f_imag.max()
    f_imag[(f_imag < ee_min) & (f_imag > -ee_min)] = 0

    print("ee_min:\n", ee_min)
    new_f = f_real + f_imag
    new_y = np.fft.ifft(new_f)
    new_y = np.real(new_y)

    print("恢復的頻率訊號:\n", new_f)
    print("恢復的去0後的頻率訊號:\n", delete_zero(new_f))

    plt.figure(figsize=(8, 8), facecolor='w')
    plt.subplot(311)
    plt.plot(x, y, 'g-', linewidth=2)
    plt.title('triangle signal', fontsize=15)
    plt.grid(True)
    plt.subplot(312)
    w = np.arange(N) * 2*np.pi / N
    plt.stem(w, a, lineformat='r-', markerformat='ro')
    plt.title('frequency domain signal', fontsize=15)
    plt.grid(True)
    plt.subplot(313)
    plt.plot(x, new_y, 'b-', lw=2, markersize=4)
    plt.title('triangle restore signal', fontsize=15)
    plt.grid(True)
    plt.tight_layout(1.5, rect=[0, 0.04, 1, 0.96])
    plt.suptitle('quickly fft translation and frequency fliter signal', fontsize=17)
    plt.show()



  傅立葉變換是訊號處理領域的一種極其重要的手段,很多的濾波演算法也是基於傅立葉變換而來,通過這個實驗我們可以簡單模擬傅立葉變換極其逆變換,同時也可以看出低通和高通濾波無非也是講個別訊號抑制,通過設定閾值在濾波中濾掉。