from math import log, ceil, cos, sin, pi
import matplotlib.pyplot as plt
import numpy as np
 
 
# 這兩行代碼解決 plt 中文顯示的問(wèn)題
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
 
 
def fft(x, N=None):
    # DIT-FFT 函數(shù)說(shuō)明
    # x: 時(shí)域序列
    # N: N點(diǎn)DFT, 理論上N=2**M
    # 返回值為序列x的DFT
    if N is None:
        N = len(x)
    elif N < len(x):
        N = len(x)
    
    if N == 2:
        return [x[0]+x[1], x[0]-x[1]]
    
    # 補(bǔ)0使得N=2**M
    M = ceil(log(N, 2))
    N = 2**M
    x = x + [0] * (N-len(x))
    
    # 遞歸地計(jì)算偶數(shù)項(xiàng)和奇數(shù)項(xiàng)的DFT
    X1 = fft(x[0::2])
    X2 = fft(x[1::2])
    X = [0] * N
    for i in range(N//2):
        # 蝶形計(jì)算
        tmp = (cos(2*pi/N*i)-1j*sin(2*pi/N*i))*X2[i]
        X[i] = X1[i] + tmp
        X[i+N//2] = X1[i] - tmp
    
    return X
 
 
if __name__ == '__main__':
    x = [1]*10
    y = fft(x, 1024)
    # print(y)
    z = [abs(i) for i in y]
    # print(z)
    plt.plot(np.arange(len(z))*2/len(z), z, label='10點(diǎn)矩形窗函數(shù)的FFT')
    plt.title("幅度譜")
    plt.xlabel(r'單位：$\pi$')
    plt.ylabel(r'$|H(j\omega)|$')
    plt.grid(linestyle="-.")
    plt.legend()
    plt.show()

from math import log, ceil, cos, sin, pi
import matplotlib.pyplot as plt
import numpy as np
 
 
# 這兩行代碼解決 plt 中文顯示的問(wèn)題
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
 
 
def fft(x, N=None):
    # DIT-FFT 函數(shù)說(shuō)明
    # x: 時(shí)域序列
    # N: N點(diǎn)DFT, 理論上N=2**M
    # 返回值為序列x的DFT
    """
    采用流式計(jì)算方法，只用了一個(gè)N(N=2**M)點(diǎn)的數(shù)組內(nèi)存
    """
    if N is None:
        N = len(x)
    elif N < len(x):
        N = len(x)
    
    # 補(bǔ)0使得：N=2**M
    M = ceil(log(N, 2))
    N = 2**M
    x = x + [0] * (N-len(x))
    
    fm = "{:0"+f"{M}"+"b}"
    X = [0] * N
    for i in range(N//2):
        index1 = eval('0b'+fm.format(i*2)[::-1])
        index2 = eval('0b'+fm.format(i*2+1)[::-1])
        X[2*i] = x[index1] + x[index2]
        X[2*i+1] = x[index1] - x[index2]
    
    for i in range(1, M):
        # 第i步表示將2**i點(diǎn)DFT合成2**(i+1)點(diǎn)的DFT
        # 蝶形寬度width
        width = 2**i
        
        """
        將X(k)序列進(jìn)行分組，每組2**(i+1)個(gè)點(diǎn)，
        便于將每組中兩組2**i點(diǎn)DFT合成一組2**(i+1)點(diǎn)的DFT
        """
        # num=2*width=2**(i+1), 表示每組點(diǎn)數(shù)
        num = 2*width
        # 組數(shù)groups
        groups = N//num
        
        for j in range(groups):
            # 對(duì)每組將2**i點(diǎn)DFT合成2**(i+1)=num點(diǎn)的DFT
            for k in range(num//2):
                # 旋轉(zhuǎn)因子
                W = cos(2*pi/num*k) - 1j * sin(2*pi/num*k)
                # 第j組第k個(gè)
                index = j*num + k
                tmp = W * X[index+width]    # 每個(gè)蝶形一次復(fù)數(shù)乘法
                X[index], X[index+width] = X[index]+tmp, X[index]-tmp
                
    return X
    
 
if __name__ == '__main__':
    x = [1]*10
    y = fft(x, 1024)
    # print(y)
    z = [abs(i) for i in y]
    # print(z)
    plt.plot(np.arange(len(z))*2/len(z), z, label='10點(diǎn)矩形窗函數(shù)的FFT')
    plt.title("幅度譜")
    plt.xlabel(r'單位：$\pi$')
    plt.ylabel(r'$|H(j\omega)|$')
    plt.grid(linestyle="-.")
    plt.legend()
    plt.show()

import numpy as np
from numpy.fft import fft, ifft
# from scipy.fftpack import fft, ifft
import matplotlib.pyplot as plt
 
 
# 這兩行代碼解決 plt 中文顯示的問(wèn)題
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
 
 
if __name__ == '__main__':
    x = 2*np.sin(np.pi/2*np.arange(100))+np.sin(np.pi/5*np.arange(100))
    z = [abs(i) for i in fft(x, 2048)]
    # print(z)
    L = len(z)
    plt.plot((np.arange(L)*2/L)[:L//2], z[:L//2], label='兩個(gè)不同頻率正弦信號(hào)相加的DFT')
    plt.title("幅度譜")
    plt.xlabel('$\pi$')
    plt.ylabel('$|H(j\omega)|$')
    plt.grid(linestyle="-.")
    plt.legend()
    plt.show()
    
    print('max(abs(ifft(fft(x))-x)) = ', end='')
    print(max(abs(ifft(fft(x))-x)))

max(abs(ifft(fft(x))-x)) = 9.01467522361575e-16