Radix-4 FFT IFFT inverse FFT 基4 的FFT python 版本

最新推荐文章于 2024-05-10 09:30:00 发布
最新推荐文章于 2024-05-10 09:30:00 发布
阅读量1.9k 收藏 10
点赞数 1
文章标签： 算法 信号处理 python
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Maxwave13/article/details/105620518
版权

Radix-2 FFT 网上的资源很多，但是Radix-4 FFT的资源很少，我只找到一个C++版本的，而且网上几乎没有 IFFT 的代码。
我实现了一个python版本的Radix-4，包括正变换，和反变换，方便理解和学习。
其中第一个版本的复杂度更低，第二个版本更方便理解。
建议先从第二个版本看起。
最好是先理解Radix-2 FFT再学习Radix-4的版本，Radix-2 FFT是普通的碟式算法FFT，网上很多资料。

参考资料：

首先，定义位翻转函数：

def bit_reverse(X, Y, N):
    j = 0
    N2 = N//4
    for i in range(N-1):
        if(i<j):
            X[i],X[j] = X[j],X[i]
            Y[i],Y[j] = Y[j],Y[i]

        N1 = N2
        while j>= 3 * N1:
            j = j - 3 * N1
            N1 = N1 // 4
        j = j + N1
    return X, Y

第一个版本的 Radix-4 FFT

def FFT_Radix4(X, Y, N, M):
    
    X_ = np.copy(X)
    Y_ = np.copy(Y)
    
    N2 = N
    for k in range(M):
        E = 2 * np.pi / N2
        N1 = N2     
        N2 = N2//4
        A = 0

        for j in range(N2):
            A1 = j * E
            A2 = 2 * A1
            A3 = 3 * A1

            CO1, CO2, CO3 = np.cos(A1), np.cos(A2), np.cos(A3)
            SI1, SI2, SI3 = np.sin(A1), np.sin(A2), np.sin(A3)
  
            for i in range(j,N,N1):

                i0 = i + 0 * N2
                i1 = i + 1 * N2
                i2 = i + 2 * N2
                i3 = i + 3 * N2            

                R1 = X_[i0] + X_[i2]
                R2 = X_[i1] + X_[i3]
                R3 = X_[i0] - X_[i2]
                R4 = X_[i1] - X_[i3]
                S1 = Y_[i0] + Y_[i2]
                S2 = Y_[i1] + Y_[i3]
                S3 = Y_[i0] - Y_[i2]
                S4 = Y_[i1] - Y_[i3]


                X_[i0] = R1 + R2;
                R2 = R1 - R2;
                R1 = R3 - S4;
                R3 = R3 + S4;
                Y_[i0] = S1 + S2;
                S2 = S1 - S2;
                S1 = S3 + R4;
                S3 = S3 - R4;

                X_[i1] = CO1*R3 + SI1*S3;
                Y_[i1] = CO1*S3 - SI1*R3;
                X_[i2] = CO2*R2 + SI2*S2;
                Y_[i2] = CO2*S2 - SI2*R2;
                X_[i3] = CO3*R1 + SI3*S1;
                Y_[i3] = CO3*S1 - SI3*R1;
    X_, Y_ = bit_reverse(X_, Y_, N)
    return X_, Y_

写反变换函数，这里主要修改了两个地方，一个是 E = -2*pi /N2加了个负号，第二就是将R3，S3与R1，S1调换。

def IFFT_Radix4(X, Y, N, M):
    
    X_ = np.copy(X)
    Y_ = np.copy(Y)
    
    N2 = N
    for k in range(M):
        E = -2 * np.pi / N2
        N1 = N2     
        N2 = N2//4
        A = 0

        for j in range(N2):
            A1 = j * E
            A2 = 2 * A1
            A3 = 3 * A1

            CO1, CO2, CO3 = np.cos(A1), np.cos(A2), np.cos(A3)
            SI1, SI2, SI3 = np.sin(A1), np.sin(A2), np.sin(A3)

            for i in range(j,N,N1):

                i0 = i + 0 * N2
                i1 = i + 1 * N2
                i2 = i + 2 * N2
                i3 = i + 3 * N2            

                R1 = X_[i0] + X_[i2]
                R2 = X_[i1] + X_[i3]
                R3 = X_[i0] - X_[i2]
                R4 = X_[i1] - X_[i3]
                S1 = Y_[i0] + Y_[i2]
                S2 = Y_[i1] + Y_[i3]
                S3 = Y_[i0] - Y_[i2]
                S4 = Y_[i1] - Y_[i3]


                X_[i0] = R1 + R2;
                R2 = R1 - R2;
                R1 = R3 - S4;
                R3 = R3 + S4;
                Y_[i0] = S1 + S2;
                S2 = S1 - S2;
                S1 = S3 + R4;
                S3 = S3 - R4;

                X_[i1] = CO1*R1 + SI1*S1;
                Y_[i1] = CO1*S1 - SI1*R1;
                X_[i2] = CO2*R2 + SI2*S2;
                Y_[i2] = CO2*S2 - SI2*R2;
                X_[i3] = CO3*R3 + SI3*S3;
                Y_[i3] = CO3*S3 - SI3*R3;
    X_, Y_ = bit_reverse(X_, Y_, N)
    return X_, Y_

另外，这里根据 反变换的共轭求法，给出一个共轭版本。

def IFFT_Radix4_conjugate(X, Y, N, M):
    Y = - Y
    res_X, res_Y = FFT_Radix4(X, Y, N, M)
    return res_X, -res_Y

接下来是实验，可以看到，经过正反变换之后，得到的结果和 numpy自带函数一致。

import numpy as np

M = 2
N = 4 ** M 
X = np.arange(N)
Y = np.arange(N, 2*N)
X = X.astype(np.float)
Y = Y.astype(np.float)


X_fft, Y_fft = FFT_Radix4(X, Y, N, M)
X_ifft, Y_ifft = IFFT_Radix4(X_fft, Y_fft, N, M)

tmp_fft = np.fft.fft(X+1j*Y)
tmp_ifft = np.fft.ifft(tmp_fft)

print('input:',X)
# print('fft2 :',np.real(tmp_fft))
# print('fft4 :',X_fft)

print('ifft2:',np.real(tmp_ifft))
print('ifft4:',X_ifft/16)

运行结果：
在这里插入图片描述
第二个版本更容易理解：
首先是正变换函数：

def FFT_Radix4_v2(X, Y, N, M):
    X_ = np.copy(X)
    Y_ = np.copy(Y)

    N2 = N
    for k in range(M):
        E = 2 * np.pi / N2
        N1 = N2     
        N2 = N2//4
        A = 0

        for j in range(N2):
            A1 = j * E
            A2 = 2 * A1
            A3 = 3 * A1
            co1, co2, co3 = np.cos(A1), np.cos(A2), np.cos(A3)
            si1, si2, si3 = np.sin(A1), np.sin(A2), np.sin(A3)
            
            CO1, CO2, CO3 = np.cos(A1), np.cos(A2), np.cos(A3)
            SI1, SI2, SI3 = np.sin(A1), np.sin(A2), np.sin(A3)
            
            for i in range(j,N,N1):
                i0 = i + 0 * N2
                i1 = i + 1 * N2
                i2 = i + 2 * N2
                i3 = i + 3 * N2     
                
                xa, xb, xc, xd = X_[i0], X_[i1], X_[i2], X_[i3]
                ya, yb, yc, yd = Y_[i0], Y_[i1], Y_[i2], Y_[i3]

                xa_ = xa + xb + xc + xd    
                ya_ = ya + yb + yc + yd   
                xb_ = (xa+yb-xc-yd)* co1 + (ya-xb-yc+xd)* (si1)    
                yb_ = (ya-xb-yc+xd)* co1 - (xa+yb-xc-yd)* (si1)    
                xc_ = (xa-xb+xc-xd)* co2 + (ya-yb+yc-yd)* (si2)    
                yc_ = (ya-yb+yc-yd)* co2 - (xa-xb+xc-xd)* (si2)    
                xd_ = (xa-yb-xc+yd)* co3 + (ya+xb-yc-xd)* (si3)    
                yd_ = (ya+xb-yc-xd)* co3 - (xa-yb-xc+yd)* (si3)   
                X_[i0], X_[i1], X_[i2], X_[i3] = xa_, xb_, xc_, xd_
                Y_[i0], Y_[i1], Y_[i2], Y_[i3] = ya_, yb_, yc_, yd_ 

    X_, Y_ = bit_reverse(X_, Y_, N)
    return X_, Y_

然后是反变换函数，修改和上面的一样。

def IFFT_Radix4_v2(X, Y, N, M):
    X_ = np.copy(X)
    Y_ = np.copy(Y)

    N2 = N
    for k in range(M):
        E = 2 * np.pi / N2
        N1 = N2     
        N2 = N2//4
        A = 0

        for j in range(N2):
            A1 = j * E
            A2 = 2 * A1
            A3 = 3 * A1
            co1, co2, co3 = np.cos(A1), np.cos(A2), np.cos(A3)
            si1, si2, si3 = np.sin(A1), np.sin(A2), np.sin(A3)
            
            CO1, CO2, CO3 = np.cos(A1), np.cos(A2), np.cos(A3)
            SI1, SI2, SI3 = np.sin(A1), np.sin(A2), np.sin(A3)
            
            for i in range(j,N,N1):
                i0 = i + 0 * N2
                i1 = i + 1 * N2
                i2 = i + 2 * N2
                i3 = i + 3 * N2     
                
                xa, xb, xc, xd = X_[i0], X_[i1], X_[i2], X_[i3]
                ya, yb, yc, yd = Y_[i0], Y_[i1], Y_[i2], Y_[i3]

                xa_ = xa + xb + xc + xd    
                ya_ = ya + yb + yc + yd    
                xb_ = (xa-yb-xc+yd)* co1 - (ya+xb-yc-xd)* (si1)    
                yb_ = (ya+xb-yc-xd)* co1 + (xa-yb-xc+yd)* (si1)    
                xc_ = (xa-xb+xc-xd)* co2 - (ya-yb+yc-yd)* (si2)    
                yc_ = (ya-yb+yc-yd)* co2 + (xa-xb+xc-xd)* (si2)    
                xd_ = (xa+yb-xc-yd)* co3 - (ya-xb-yc+xd)* (si3)    
                yd_ = (ya-xb-yc+xd)* co3 + (xa+yb-xc-yd)* (si3)    


                X_[i0], X_[i1], X_[i2], X_[i3] = xa_, xb_, xc_, xd_
                Y_[i0], Y_[i1], Y_[i2], Y_[i3] = ya_, yb_, yc_, yd_ 

    X_, Y_ = bit_reverse(X_, Y_, N)
    return X_, Y_

进行同样的实验：

import numpy as np

M = 2
N = 4 ** M 
X = np.arange(N)
Y = np.arange(N, 2*N)
X = X.astype(np.float)
Y = Y.astype(np.float)


X_fft, Y_fft = FFT_Radix4_v2(X, Y, N, M)
X_ifft, Y_ifft = IFFT_Radix4_v2(X_fft, Y_fft, N, M)

tmp_fft = np.fft.fft(X+1j*Y)
tmp_ifft = np.fft.ifft(tmp_fft)

print('input   :',X)
# print('fft2    :',np.real(tmp_fft))
# print('fft4_v2 :',X_fft)

print('ifft2   :',np.real(tmp_ifft))
print('ifft4_v2:',np.round(X_ifft/16))

下面是运行结果：
在这里插入图片描述
可以看到，两次结果均正确。
Good Job！

Maxwave13
关注
  • 1
    点赞
  • 踩
  • 10
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
FFT算法介绍及OAI中FFT的实现
qq_44113393的博客
04-23 2333
FFT算法介绍及OAI中FFT的实现OAI中FFT的实现FFT介绍DIT radix-2 FFT（时域抽取基-2 FFT）DIF radix-2 FFT (频域抽取基-2 FFT)DIF radix-4 FFT (频域抽取基-4 FFT)OAI中FFT的实现 OAI中FFT的实现 FFT介绍 FFT是DFT（离散傅立叶变换）的一种快速算法。 DFT是一种常用的信号处理方式，可以将时域信号转换到频域...
Cooley-Tukey算法 （蝶形算法）
热门推荐
xxinliu的专栏
04-08 2万+
Cooley－Tukey算法区别于其他FFT算法的一个重要事实就是N的因子可以任意选取。这样也就可以使用N＝rS的Radix－r算法了。最流行的算法都是以r＝2或r＝4为基的，最简单的DFT不需要任何乘法就可以实现。例如：在S级且r＝2的情形下，下列索引映射的结果是： 　　S＞2时的－个一般惯例是，在信号流程图中2点DFT是以蝶形图的形式绘出的，图1给出了8点变换的图示。信号流程图己经简化成
基4快速傅里叶变换，基4快速傅里叶反变换 Radix4 FFT
qq_37198521的博客
08-17 300
就将一个N 点得DFT转化为四个N/4点DFT来计算。依此类推，直至分解到最后一级。以上就是按频率抽取的基4快速傅里叶变换算法。与基2FFT相比，基4FFT的乘法量约减少25%，加法量也略有减少。
基于Radix-4的Booth乘法器
最新发布
qq_43433724的博客
05-10 1046
无符号数的最高位1被当作符号位进行运算了，这与Booth乘法器的原理不符合，接下来最高位的1将作为Bi-1参与Cofe的计算，那么就需要额外的两bit数，因此在高位需要扩展两个0，这是需要的最多的额外bit的情况。如果乘数是7，那么扩展低位后变为1110，这时候变为如下的情况，只需要1bit即可，但是扩展两位并不影响运算。因此每次计算时，只需要从最低位置开始，按照基系数Coef计算出相邻三个bit的和，然后根据和为-2、-1、0、1、2来得到乘法结果并移位，之后将所有得到的乘法结果相加即可。
radix 4 FFT verilog代码 c代码及相应uvm验证平台.zip
08-25
FFT有以下特性： l支持2^N复数点FFT/IFFT运算，其中4<= N <= 10 l支持数据input和output并行 l采用Raidx-4 butterfly设计 l支持添加循环前缀 l支持自动休眠（低功耗） 验证平台基于windos（questasim），包含与c model的自动比对
fft.js:最快的JS Radix-4Radix-2 FFT实现
05-25
FFT.js Radix-4 FFT的实现。 用法 const FFT = require ( 'fft.js' ) ; const f = new FFT ( 4096 ) ; const input = new Array ( 4096 ) ; input . fill ( 0 ) ; const out = f . createComplexArray ( ) ; 如果data只有实数（如使用toComplexArray时的情况），则可以运行实数FFT将其计算速度提高25％： const realInput = new Array ( f . size ) ; f . realTransform ( out , realInput ) ; realTransform仅填充out的左半部分，因此，如果需要整个频谱（对称的），请执行以下操作： f . completeSp
Radix4Radix2FFT实现
08-12
**Radix-4与Radix-2快速傅里叶变换(FFT)实现** 在数字信号处理和计算领域，快速傅里叶变换（Fast Fourier Transform，FFT）是一种高效计算离散傅里叶变换（DFT）及其逆变换的方法。FFT算法大大减少了计算量，尤其...
fft.rar_RADiX 2 fft_radix-2 dit
09-21
标题中的"fft.rar_RADiX 2 fft_radix-2 dit"表明这是一个关于快速傅里叶变换（FFT）的资源，特别提到了基2的直接型（Direct Inverse Transform，DIT）算法。FFT是一种高效的计算离散傅里叶变换（DFT）的方法，它将...
fft.zip_fft radix4 dif_sittingfdn
09-24
标题中的"fft.zip_fft radix4 dif_sittingfdn"表明这是一个关于快速傅里叶变换（FFT）的程序，特别地，它使用了基4的分解（radix-4）算法，并且可能涉及到一种称为“sittingfdn”的特定实现或功能。在描述中，我们...
快速傅立叶变换：使用 radix2 的 fft 算法。-matlab开发
06-01
编写此函数是为了理解和/或解释 fft 算法和 radix2 过程的用法。 此函数取您要查找 FFT 的信号，为 m 位 您也可以在不输入信号或点数的情况下调用该函数，因为 y 和 m 有默认值； y = sin(2*pi*50*t)+sin(2*pi*120*t)； 米=3； 您可以看到用作内置函数的顺序反转步骤，但没有其他步骤。 您还将看到采取蝴蝶树的步骤。 有关更多信息，请阅读函数的帮助文件。 将来这可以作为选择性基数或您建议的任何其他功能进行改进，请注意我的任何建议。
FFT原理及实现(Radix-2)
nmzhang0909的专栏
06-17 1623
基2 FFT总的思想是将输入信号对半分割, 再对半分割, 再再对半分割(以下省略10000个再再...J) 直至分割到2点. 两点DFT简化假设输入为x[0],x[1]; 输出为X[0],X[1].  伪代码如下 : // ------------------------------------------------------------------#define    N
java inverse_Java版本的FFT和Inverse FFT | 学步园
weixin_42468413的博客
02-24 217
最近有些朋友在一些项目中需要用到Java版本的FFT和Inverse FFT，玄机逸士在网上找到了一个版本，供大家参考，现抄录如下：/**************************************************************************Compilation:javac FFT.java*Execution:java FFT N*Dependencie...
Java版本的FFT和Inverse FFT
玄机逸士的专栏
05-04 7656
最近有些朋友在一些项目中需要用到Java版本的FFT和Inverse FFT，玄机逸士在网上找到了一个版本，供大家参考，现抄录如下：（原文地址：http://www.cs.princeton.edu/introcs/97data/FFT.java.html）/******************************************************************
基4的DIF的fft算法全部弄懂
GoodGoodStudy
08-25 1万+
刚刚想清楚旋转因子的问题，现在算法已经全部清楚了，下面看汇编程序。
4.8 IFFT/FFT
fighting2019的博客
12-14 1万+
4.8.1 IFFT/FFT原理 1. IFFT、FFT在OFDM系统中的应用 2、IFFT/FFT原理  3、DIT的基2FFT算法    4、DIF的基2FFT运算  5、基于的频率抽选FFT算法   4.8.2 基于 DIF FFT的硬件 结构 1、各级的蝶形运算硬件实现    2、transfer12、 transfer34、trans...
从头到尾彻底理解傅里叶变换算法
05-03 6452
http://blog.csdn.net/v_JULY_v/article/details/6196862 作者：July、dznlong   二零一一年二月二十日 前言 第一部分、  DFT 第一章、傅立叶变换的由来 第二章、实数形式离散傅立叶变换（Real DFT） 从头到尾彻底理解傅里叶变换算法、下 第三章、复数 第四章、复数形式离散傅立叶变换   前言： “关于傅立叶
数学专题小结：FFT算法
Skyline
09-02 1万+
快速傅里叶变换（FFT，Fast Fourier Transform）是信号处理的常用手段，可以把时域信号变成频域信号，时域的卷积运算对应于频域就成了简单的乘法运算。由于两个多项式的乘积，其系数的运算实际上也是一种卷积运算，因此可以用FFT来计算多项式的乘法。网上关于FFT算法的讲解大多都是工程领域的，这里从算法竞赛的角度来剖析一下怎么简洁而高效的实现FFT。 一，为什么要用FFT来计算
短时逆傅里叶变换inverse short-time fft
weixin_40304245的博客
04-03 2974
Matlab短时逆傅里叶变换函数 本节内容部分来自于Matlab参考文档 Matlab中短时逆傅里叶函数是signal processing toolbox工具箱中的istft函数，调用格式x = istft(s,fs,Name,Value)； 若想使istft后的函数完美复现原信号，必须满足以下三点： 输入信号：要满足k = (length(x) - noverlap)/(length(win...
radix-8 fft 代码实现
08-08
下面是一个简单的示例，展示如何实现基于 radix-8 的快速傅里叶变换（FFT）的代码： ```python import numpy as np def radix8fft(x): N = len(x) if N <= 1: return x else: even = radix8fft(x[0::2]) odd = radix8fft(x[1::2]) twiddle = np.exp(-2j * np.pi * np.arange(N) / N) x = np.concatenate([even + twiddle[:N//8] * odd, even + twiddle[N//8:N//4] * odd, even + twiddle[N//4:3*N//8] * odd, even + twiddle[3*N//8:N//2] * odd, even + twiddle[N//2:5*N//8] * odd, even + twiddle[5*N//8:3*N//4] * odd, even + twiddle[3*N//4:7*N//8] * odd, even + twiddle[7*N//8:] * odd]) return x # 示例输入信号 x = np.array([0, 1, 2, 3, 4, 5, 6, 7]) # 应用 radix-8 FFT X = radix8fft(x) # 打印结果 print(X) ``` 这个代码示例使用了递归的方式来实现 radix-8 FFT。在每次递归中，输入序列被分成偶数索引和奇数索引的两个子序列，然后通过应用旋转因子（twiddle factors）将它们重新组合起来。最终得到的结果就是输入序列的 radix-8 FFT。 请注意，这只是一个简单的示例，实际的 radix-8 FFT 实现可能会更复杂，并且通常会使用更高效的算法和优化技术来提高性能。
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值