FFT代码编写（python）

小谢啥都不会

已于 2023-10-16 15:08:08 修改

阅读量898

点赞数 4

文章标签： python

于 2023-10-15 22:56:05 首次发布

本文链接：https://blog.csdn.net/qq_67725674/article/details/133844362

版权

文章目录

前言

数字信号处理刚学完FFT，老师让我们编程来加深对FFT的理解

一、FFT代码设计步骤

（1）判断

判断输入是否为2的整次幂，若不为则填充0。

# 检查一个整数是否是2的整次幂
def is_power_of_two(n):
    m=1
    while(True):
        m=m*2
        if m==n:
            return True
        elif m>n:
            return False

def pad_to_power_of_two(x):
    N = len(x)
    if is_power_of_two(N):
        # 如果输入长度已经是2的整次幂，不需要填充，直接返回
        return x
    else:
        # 找到最接近N且大于N的2的整次幂
        next= 1
        while next< N:
            next*= 2
        # 填充零值到输入列表的末尾
        padded_x = x + [0] *(next - N)
        return padded_x

这里判断2的整次幂还有一种简单的方法：n & (n - 1)

（2）比特翻转

def bit_reverse(x):
    N = len(x)
    num_bits =int(math.log2(N))
    # 通过位反转重新排列输入数据
    reversed_x = [0] * N
    for i in range(N):
        reversed_i = int(f"{i:0{num_bits}b}"[::-1], 2)
        reversed_x[i] = x[reversed_i]
    return reversed_x

reversed_i = int(f"{i:0{num_bits}b}"[::-1], 2) 是一个将整数索引 i 转换为比特反转后的索引 reversed_i 的方法，它的具体步骤如下：

f"{i:0{num_bits}b}" 这部分是将整数 i 格式化为一个二进制字符串。{i:0{num_bits}b} 的意思是将整数 i 转换为一个二进制字符串，其中 num_bits 表示二进制字符串的最小宽度，前面会用零填充。这确保了二进制字符串的长度为 num_bits。

[:: -1] 是字符串切片，它将生成一个颠倒顺序的字符串。即，它将二进制字符串中的字符反转。

最后，int(..., 2) 用于将反转后的二进制字符串转换回整数。2 表示我们将使用二进制表示法进行转换

（3）FFT核心计算

def fft(x):
    N = len(x)
    # 位反转重新排列输入数据
    x = bit_reverse(x)
    for s in range(1, N.bit_length()):
        m = 2 ** s
        for k in range(0, N, m):#k当前组的起始索引
            for j in range(0, m // 2):#每组的一半
                if (k + j + m // 2==N):
                    break;
                W= cmath.exp(-2j * cmath.pi * j*(2 ** (N.bit_length()-s-1)) / N)
                t = W* x[k + j + m // 2]
                u = x[k + j]
                x[k + j] = u + t
                x[k + j + m // 2] = u - t
    return x

二、其它

（4）绘图

def visual(x,X):
    # 计算频谱（幅度谱）
    X_magnitude = np.abs(X)
    # 生成频率坐标轴
    N = len(x)
    #freq = np.fft.fftfreq(N)#范围从负半频率（负采样率的一半）到正半频率（正采样率的一半）
    # 创建频谱图
    plt.figure(figsize=(10, 6))
    plt.subplot(2, 1, 1)#两行子图，一列子图，改子图放第一
    plt.stem(x, use_line_collection=True)
    plt.title('Input Signal')
    plt.subplot(2, 1, 2)
    plt.stem(X_magnitude, use_line_collection=True)
    plt.title('FFT Spectrum')
    plt.tight_layout()
    plt.show()