精华内容
下载资源
问答
  • 基2基3定点混合基fft

    2010-12-14 15:26:35
    定点基2,基3 的混合基fft,没有加入基5的部分,不过留的有空,有兴趣的谈朋友可以自己加一下
  • 三、FFT的C语言实现(时域) 四、FFT的C语言实现(频域) ---------------------------------------------------------------------------------------------------------------- 一、理解离散傅立叶变换...
    目录:
    一、理解离散傅立叶变换(FFT结果的物理意义(编程参考用))
    二、基二频域算法原理
    三、基二FFT的C语言实现(时域)
    四、基二FFT的C语言实现(频域)
     
    ----------------------------------------------------------------------------------------------------------------
    一、理解离散傅立叶变换(FFT结果的物理意义(编程参考用))
    ----------------------------------------------------------------------------------------------------------------

    二、基二频域算法原理 

    原理请查看“按时间抽取基2的FFT算法的实现--杜义君” 

    ----------------------------------------------------------------------------------------------------------------

    三、基二FFT的C语言实现(时域) 

    #include "math.h"
    #include "stdio.h"

    struct compx
    { double real;
      double imag;
    } compx ;


    struct compx EE(struct compx b1,struct compx b2)//复数相乘
    {
    struct compx b3;
    b3.real=b1.real*b2.real-b1.imag*b2.imag;
    b3.imag=b1.real*b2.imag+b1.imag*b2.real;
    return(b3);
    }

    void FFT(struct compx *xin,int N)
    {
    int f,m,LH,nm,i,k,j,L;
    double p , ps ;
    int le,B,ip;
    float pi;
    struct compx v,w,t;
    LH=N/2;
    f=N;
    for(m=1;(f=f/2)!=1;m++){;}//求出m为log2 N
    nm=N-2;   
    j=N/2;
    //变址运算,对时间进行奇偶分解

    for(i=1;i<=nm;i++)//即xin第一位和最后一位不用操作,不用变址,其余各位根据码位倒置
    {
    if(i
    k=LH;
    while(j>=k){j=j-k;k=k/2;}
    j=j+k;
    }

    {
    for(L=1;L<=m;L++)//运行m级蝶形运算

    le=pow(2,L);
    B=le/2;
    pi=3.14159;
     for(j=0;j<=B-1;j++)
      {
       p=pow(2,m-L)*j;//k的取值(0,1,2,...,(pow(2,L)/2)-1),当在第一级时,p只为0,当在第二

    //级时,p为0和2,当为第三级时,p为0,1,2,3.在递升的级别中不断延续,见数字信号处理(华中

    //科技大学p83图的系数WN),且相邻不同种基本蝶形的蝶形节系数的增量为(2*pi/N)*pow(2,m-L)
       ps=(2*pi/N)*p;
       w.real=cos(ps);
       w.imag=-sin(ps);//求出WNk
       for(i=j;i<=N-1;i=i+le)//确定与蝶形系数相乘的Xm(q)的下标m,此方法为蝶形图的特点来得到,

    //相邻同种基本蝶形的间距为2的L次方。
         {
          ip=i+B;//即对频谱进行前后分解
          t=EE(xin[ip],w);//复数相乘
          xin[ip].real=xin[i].real-t.real;//基本蝶式运算
          xin[ip].imag=xin[i].imag-t.imag;
          xin[i].real=xin[i].real+t.real;
          xin[i].imag=xin[i].imag+t.imag;
         }
      }
    }
    }  
    return ;
    }


    //输入时域数据点为num,则输出频域数据点同为num,num是数据长度,必须为2的整数次幕,

    //其大小由数据采样定理来决定

    #include
    #include
    #include

    float  result[257];//振幅,其平方为功率谱
    struct  compx s[257];
    int   Num=16;//数据长度,必须为2的整数次幕
    const float pp=3.14159;

    main()
    {

    int i;
    for(i=0;i<16;i++)
    {
    s[i].real=sin(pp*i/32);
    s[i].imag=0;
    }

    FFT(s,Num);

    for(i=0;i<16;i++)
    {
    printf("%.4f",s[i].real);
    printf("+%.4fj\n",s[i].imag);
    result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));//pow功 能: 指数函数(x的y次方) 用 法: double pow(double x, double y);
    }


    }

     

    时域公式

     

     

     

    频域公式

     

    ----------------------------------------------------------------------------------------------------------------

    四、基二FFT的C语言实现(频域)

        基于时域和基于频域的算法很相似的,只不过时域里是先乘后加减,对输入序列进行倒序,而频域的算法先加减后乘,输入序列不用倒序,对输出序列进行倒序。调试成功之后,发现两者结果相差很小了。

    #include "math.h"
    #include "stdio.h"

    struct compx
    { double real;
      double imag;
    } compx ;


    struct compx EE(struct compx b1,struct compx b2)
    {
    struct compx b3;
    b3.real=b1.real*b2.real-b1.imag*b2.imag;
    b3.imag=b1.real*b2.imag+b1.imag*b2.real;
    return(b3);
    }

    void FFT(struct compx *xin,int N)
    {
    int f,m,LH,nm,i,k,j,L;
    double p , ps ;
    int le,B,ip;
    float pi;
    struct compx v,w,t;
    LH=N/2; 
    f=N;
    for(m=1;(f=f/2)!=1;m++){;}  //2^m=N{
    for(L=m;L>=1;L--)    //这里和时域的也有差别


    le=pow(2,L);
    B=le/2; //每一级碟形运算间隔的点数

    pi=3.14159;
     for(j=0;j<=B-1;j++)
      {
       p=pow(2,m-L)*j;
       ps=2*pi/N*p;
       w.real=cos(ps);
       w.imag=-sin(ps);
       for(i=j;i<=N-1;i=i+le)
         {
          ip=i+B;  
          t=xin[i];
          xin[i].real=xin[i].real+xin[ip].real;
          xin[i].imag=xin[i].imag+xin[ip].imag;  
          xin[ip].real=xin[ip].real-t.real;
          xin[ip].imag=xin[ip].imag-t.imag;     
          xin[ip]=EE(xin[ip],w);
         }
      }
    }
    }
    //变址运算

    nm=N-2;   
    j=N/2;
    for(i=1;i<=nm;i++)
    {
    if(i
    k=LH;
    while(j>=k){j=j-k;k=k/2;}
    j=j+k;
    }

    }


    //main programe

    #include
    #include
    #include

    float  result[257];
    struct  compx s[257];
    int   Num=16;
    const float pp=3.14159;

    main()
    {

    int i;
    for(i=0;i<16;i++)
    {
    s[i].real=sin(pp*i/32);
    s[i].imag=0;
    }

    FFT(s,Num);

    for(i=0;i<16;i++)
    {
    printf("%.4f",s[i].real);
    printf("+%.4fj\n",s[i].imag);
    result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
    }

    }

    ----------------------------------------------------------------------------------------------------------------

    展开全文
  • 本文代码中FFT使用递归版本实现 FFT加速多项式乘法原理不多说了,直接贴代码如下: 在vs2017上测试成功 #include "pch.h" #define _CRT_SECURE_NO_WARNINGS #include "stdlib.h" #include "math.h" #include ...

    本文代码中FFT使用递归版本实现

    FFT加速多项式乘法原理不多说了,直接贴代码如下:

    在vs2017上测试成功

    #include "pch.h"
    #define _CRT_SECURE_NO_WARNINGS
    #include "stdlib.h"
    #include "math.h"
    #include "stdio.h"
    
    
    #define N 8
    #define MAXN 100
    
    #define Pi  3.1415927   //定义圆周率Pi
    #define LEN sizeof(struct Compx)  //定义复数结构体大小
    
    //-----定义复数结构体-----------------------
    typedef struct Compx
    {
    	double real;
    	double imag;
    }Compx;
    
    //-----复数乘法运算函数---------------------
    struct Compx mult(struct Compx b1, struct Compx b2)
    {
    	struct Compx b3;
    	b3.real = b1.real*b2.real - b1.imag*b2.imag;
    	b3.imag = b1.real*b2.imag + b1.imag*b2.real;
    	return(b3);
    }
    
    //-----复数减法运算函数---------------------
    struct Compx sub(struct Compx a, struct Compx b)
    {
    	struct Compx c;
    	c.real = a.real - b.real;
    	c.imag = a.imag - b.imag;
    	return(c);
    }
    
    //-----复数加法运算函数---------------------
    struct Compx add(struct Compx a, struct Compx b)
    {
    	struct Compx c;
    	c.real = a.real + b.real;
    	c.imag = a.imag + b.imag;
    	return(c);
    }
    
    void fft(Compx *a, int n, int inv);
    
    int main()
    {
    	int i;
    
    	int x[N] = { 0 }, y[N] = { 0 };
    
    	printf("\nN=%d\n" , N);
    	printf("\n输入两个多项式的系数,输入系数为N多项式长度的一半\n");
    	printf("\n输入第一个多项式的系数\n");
    	for (i = 0; i < N/2; i++)
    	{
    		scanf("%d", &x[i]);
    	}
    
    	printf("\n输入第二个多项式的系数\n");
    	for (i = 0; i < N/2; i++)
    	{
    		scanf("%d", &y[i]);
    	}
    
    	struct  Compx * a = (struct Compx *)malloc(N*LEN);			//为结构体分配存储空间
    	struct  Compx * b = (struct Compx *)malloc(N*LEN);
    	struct  Compx * c = (struct Compx *)malloc(N*LEN);
    
    	//初始化======================================= 
    	printf("\na多项式初始化:\n");
    	for (i = 0; i < N; i++)
    	{
    		a[i].real = x[i];
    		a[i].imag = 0;
    		printf("%.4f ", a[i].real);
    		printf("+%.4fj  ", a[i].imag);
    		printf("\n");
    	}
    	printf("\nb多项式初始化:\n");
    	for (i = 0; i < N; i++)
    	{
    		b[i].real = y[i];
    		b[i].imag = 0;
    		printf("%.4f ", b[i].real);
    		printf("+%.4fj  ", b[i].imag);
    		printf("\n");
    	}
    	printf("\n结果c多项式初始化:\n");
    	for (i = 0; i < N; i++)
    	{
    		c[i].real = 0;
    		c[i].imag = 0;
    		printf("%.4f ", c[i].real);
    		printf("+%.4fj  ", c[i].imag);
    		printf("\n");
    	}
    
    	fft(a, N, 1);
    
    	printf("\n第一个多项式FFT计算结果:\n");
    	for (i = 0; i < N; i++)
    	{
    		printf("%.4f ", a[i].real);
    		printf("+%.4fj  ", a[i].imag);
    		printf("\n");
    	}
    
    	fft(b, N, 1);
    
    	printf("\n第二个多项式FFT计算结果:\n");
    	for (i = 0; i < N; i++)
    	{
    		printf("%.4f ", b[i].real);
    		printf("+%.4fj  ", b[i].imag);
    		printf("\n");
    	}
    
    	for (i = 0; i < N; i++)
    		a[i] = mult(a[i] , b[i]);
    
    	fft(a, N, -1);
    	for (i = 0; i < N; i++) {
    		c[i].real = a[i].real / N;
    		c[i].imag = a[i].imag / N;
    	}
    		
    	printf("\n乘积多项式结果:\n");
    	for (i = 0; i < N; i++)
    	{
    		printf("%.4f ", c[i].real);
    		printf("+%.4fj  ", c[i].imag);
    		printf("\n");
    	}
    
    	return 0;
    }
    
    void fft(Compx *a, int n, int inv) {
    	if (n == 1)return;
    	int mid = n / 2;
    	static Compx b[MAXN];
    	int i;
    	for (i = 0; i < mid; i++) {
    		b[i] = a[i * 2];
    		b[i + mid] = a[i * 2 + 1];
    	}
    		
    	for (i = 0; i < n; i++)
    		a[i] = b[i];
    	fft(a, mid, inv);
    	fft(a + mid, mid, inv);//分治
    	for (i = 0; i < mid; i++)
    	{
    		Compx x;
    		x.real = cos(2 * Pi*i / n);
    		x.imag = inv * sin(2 * Pi*i / n);
    
    		b[i] = add(a[i], mult(x, a[i + mid]));
    		b[i + mid] = sub(a[i] , mult(x , a[i + mid]));
    	}
    	for (i = 0; i < n; i++)
    		a[i] = b[i];
    }

     

    展开全文
  • 浮点混合基fft程序

    2010-12-14 15:24:02
    基2,基3的混合基fft,基5的部分没有加入,不过留的有地方,有兴趣的可以自己加一下,跟基2基3部分差不多
  • 快速傅立叶变换(FFT)是信号处理和数据分析中最重要的算法之一,很多人只是调用现成的库如FFTW,但为了知其所以然,加深对算法的理解,我们有必要搞懂FFT算法是怎么计算的,这里不讨论傅里叶变换的理论和推导,只...

    FFT算法

    快速傅立叶变换(FFT)是信号处理和数据分析中最重要的算法之一,很多人只是调用现成的库如FFTW,但为了知其所以然,加深对算法的理解,我们有必要搞懂FFT算法是怎么计算的,这里不讨论傅里叶变换的理论和推导,只讨论实际工程中怎样计算,由于python代码的可读性以及计算的方便性,使用python代码展示FFT计算过程

    傅里叶变换

    傅里叶变换FT(fourier transform)用于将时域信号x(t)和频域信号X(f)之间变换,公式如下所示:

    X(f)=x(t)ej2πftdt=X(f)ej2πftdfX(f) = \int^{\infty}_{-\infty}x(t)e^{-j2\pi ft}dt = \int^{\infty}_{-\infty}X(f)e^{j2\pi ft}df

    离散傅里叶变换

    由于计算机只能处理有限长的离散信号,因此
    因此必须建立对应的离散傅里叶变换DFT(Discrete Fourier Transform):

    Xk=n=0N1xnei2πknNX_k = \sum\limits_{n=0}^{N-1} x_n \cdot e^{-i2 \pi k \frac{n}{N}}

    如果我们定义一个矩阵MM
    Mkn=ei2πknNM_{kn} = e^{-i2 \pi k \frac{n}{N}}

    则很明显DFT的公式只是一个简单的线性变换:
    X=xMX = x \cdot M
    因此简单的使用矩阵乘法就能计算出DFT的结果,我们可以很容易的写出DFT的python代码

    import numpy as np
    
    def DFT(x):
        N = len(x)
        k,n = np.meshgrid(np.arange(N),np.arange(N))
        W = np.exp(-1j*2*np.pi*k*n/N)   
        return np.dot(x,W)
    

    我们以2048点DFT为例,与numpy中内置的FFT做对比,看看速度相差多少

    x=np.random.random(2048)
    X=np.fft.fft(x)
    X=DFT(x)
    
    compute 2048 points dft using np.fft cost: 0.001677 ms
    compute 2048 points dft using DFT cost: 0.321912 ms
    

    可以看到,速度相差了差不多2000倍,对于每个值X(k)X(k)的计算需要N个复数乘法(4N个乘法和2N个加法)和N-1个复数加法(2N-2个加法),因此DFT的总计算量需要N2N^2个复数乘法和N2NN^2-N个复数加法,复杂度是O[N2]\mathcal{O}[N^2],是不利于计算机进行实时信号处理的,因此为了优化DFT的计算量,便有了相关FFT算法,下面介绍快速傅里叶变换算法,对于快速傅里叶逆变换其优化方式非常相似,因此不做介绍

    快速傅里叶变换

    Cooley-Tukey快速傅里叶算法是常见的FFT算法,其思想是利用了DFT变换中的对称性和周期性来简化计算

    首先我们定义
    WN=ei2πNW_N=e^{-i \frac{2 \pi}{N}}
    WNW_N满足如下性质
    周期性:WNk+N=WNkW_N^{k+N} = W_N^k
    对称性:WNk+N2=WNkW_N^{k+\frac {N}{2}} = -W_N^k
    mmNN的约数:WNmkn=WNmknW_N^{mkn} = W_{\frac{N}{m}}^{kn}

    我们只需几行代码就可验证上述特性

    def Wn(k,N):
            return np.exp(-1j*2*np.pi*k/N)
    

    定义如下一些变量:

    N = 8
    k = 3
    m = 2
    n = 2
    

    验证周期性:

    print(np.allclose(Wn(k,N),-Wn(k+N,,N)))
    

    验证对称性:

    print(np.allclose(Wn(k,N),-Wn(k+N//2,,N)))
    

    验证可约性:

    print(np.allclose(Wn(m*k*n,N),Wn(k*n,N//m)))
    

    结果:

    True
    True
    True
    

    基2 FFT

    根据上面的对称性,我们可以将DFT计算分为两个较小的部分

    Xk=n=0N1xnWNknX_k = \sum\limits_{n=0}^{N-1} x_n \cdot W_N^{kn}
    =m=0N/21x2mWN2mk+m=0N/21x2m+1WN(2m+1)k\quad = \sum\limits_{m=0}^{N/2 - 1} x_{2m} \cdot W_N^{2mk} + \sum_{m=0}^{N/2 - 1} x_{2m + 1} \cdot W_N^{(2m+1)k}
    =m=0N/21x2mWN2km+WNkm=0N/21x2m+1WN2km\quad = \sum\limits_{m=0}^{N/2 - 1} x_{2m} \cdot W_{\frac{N}{2}}^{km} +W_N^k \sum\limits_{m=0}^{N/2 - 1} x_{2m + 1} \cdot W_{\frac{N}{2}}^{km}
    =F1(k)+WNkF2(k)\quad = F_1(k)+W_N^kF_2(k)
    这样一个N点变换就分解为了两个N/2点变换,这里F1(k)F_1(k)F2(k)F_2(k)分别是序列x中的奇数号和偶数号序列的N/2N/2点DFT变换,根据以上公式我们也能很快写出python代码:

    def R2FFT(x):
       N = len(x)
       N2 = N // 2
       k,n = np.meshgrid(np.arange(N2),np.arange(N2))
       W = np.exp(-1j*2*np.pi*k*n/N2)
       G = np.exp(-2j * np.pi * np.arange(N2) / N)
       X_even = np.dot(x[::2],W)
       X_odd = G*np.dot(x[1::2],W)
       return np.concatenate([X_even+X_odd,X_even-X_odd])
    

    同样计算2048点DFT,速度如下:

    compute 2048 points dft using R2FFT cost: 0.081140 ms
    

    对于N=2rN=2^r,很显然两个N/2点的DFT变换还可以继续分解下去,分解为4个N/4的更短的序列,N/4的序列还可以将序列继续分解下去,直到分解为N/2个2点的DFT变换,2点的DFT变换只需要复数加法和减法就能实现,复数乘法计算量减小至(N/2)logN(N/2)logN,复数加法计算量减小至NlogNNlogN,算法复杂度为O[NlogN]\mathcal{O}[NlogN],大大减少了DFT的计算量,这就是Cooley-Tukey快速傅里叶变换的基本原理,我们将一个DFT变换分解为两个较小的DFT变换,即基2FFT,我们可以通过递归来实现该算法:

    def RecursiveR2FFT(x):
        N = x.shape[0]
        if N <= 2:
            return [x[0]+x[1],x[0]-x[1]]
        else:
            X_even = RecursiveR2FFT(x[::2])
            X_odd = RecursiveR2FFT(x[1::2])
            factor = np.exp(-2j * np.pi * np.arange(N//2) / N)
            return np.concatenate([X_even + factor * X_odd,
                                   X_even - factor * X_odd])
    

    计算2048点DFT速度如下:

    compute 2048 points dft using RecursiveR2FFT cost: 0.081140 ms
    

    相比前面的版本速度并没有提升,是因为python的递归版本并不高效,并且没有进行并行化的计算,因此,通过观察基2fft的规律我们可以将递归调用的向量乘法转换为并行计算的矩阵乘法以删除递归调用以及并行计算,python代码如下:

    def NonRecursiveR2FFT(x):
        L = len(x)
        N_base = 2
        base = L//N_base
        X = np.reshape(x,(-1,base))
        X = np.vstack([X[0]+X[1],X[0]-X[1]]).T
        for n in range(int(np.log2(base))):
            N = X.shape[1]
            W = np.exp(-1j * np.pi * np.arange(N) / N)
            X_even = X[:X.shape[0]//2]
            X_odd = W*X[X.shape[0]//2:]
            X = np.concatenate([X_even+X_odd,X_even-X_odd],axis=-1)
        return X.ravel()
    

    计算2048点DFT速度如下:

    compute 2048 points dft using NonRecursiveR2FFT cost: 0.000327 ms
    

    可以看到,速度又提高了一个数量级,相比numpy的fft只差了1倍

    基4FFT

    当DFT点数N为4的幂时,我们当然可以使用基2FFT算法进行计算,但对于这种情况使用基4FFT算法更为高效,基4FFT的原理与基2FFT类似,只不过是将N点DFT序列拆分成4个N/4的子序列:
    Xk=n=0N1xnWNknX_k = \sum\limits_{n=0}^{N-1} x_n \cdot W_N^{kn}
    =WN0F0(k)+WNkF1(k)+WN2kF2(k)+WN3kF3(k)\quad = W_N^0F_0(k)+W_N^kF_1(k)+W_N^{2k}F_2(k)+W_N^{3k}F_3(k)
    在这里直接给出非递归的基4FFT代码:

    def NonRecursiveR4FFT(x):
        L = len(x)
        N_base = 2
        base = L//N_base
        X = np.reshape(x,(-1,base))
        X = np.vstack([X[0]+X[1],X[0]-X[1]]).T
        butterfly_matrix = np.array([[1,1,1,1],[1,-1j,-1,1j],[1,-1,1,-1],[1,1j,-1,-1j]])
        for n in range(int(np.log(base)/np.log(4))):
            N4 = X.shape[1]
            N = N4*4
            G = np.exp(-2j * np.pi * (np.arange(N4).reshape(1,-1)*np.arange(4).reshape(-1,1)) / N)
            X = G*X.reshape((4,-1,N4)).transpose([1,0,2])
            X = np.dot(butterfly_matrix,X).transpose([1,0,2]).reshape((-1,N))
        return X.ravel()
    

    注意,由于N=2048=245N=2048=2\cdot 4^5,因此我们最后分成了1024个2点FFT,如果NN是4的幂例如N=1024,那么最后会得到2个512点的结果,并不满足基4FFT的条件,那么我们可以将这2个512点序列按照基2FFT原理进行计算,最终得到1024个FFT点的计算结果,这实际上是一个混合了基2FFt和基4FFT的混合基FFT算法。

    结论

    由于我们的python算法涉及较多的内存分配和复制等操作,因此内存利用率和效率方面还是较numpy的fft差了一些,在实际的FFT算法中会使用C或者汇编来编写高度优化的FFT算法,例如会使用蝶形运算结构来优化计算流程,以最小化内存的使用,对于其他长度的序列,有更复杂的混合基FFT以及一些其他算法。

    使用python实现该算法的目的也并不是为了效率,而是因为python代码在可读性方面远优于使用C或者汇编。

    尽管在工程中我们有许多开源的FFT库可以使用,但是我相信理解FFT这种基础的算法以及其工程实现是非常有必要的,在必要的情况下,我们也可以按照自己的需求实现特定的FFT算法,这对于我们的技术能力的提升以及知识的积累和沉淀是有非常大的帮助的。

    欢迎感兴趣的同学加微信交流:xu19315519614

    展开全文
  • 混合3780点FFT

    2012-09-04 17:12:17
    自己在matlab中写的3780点FFT的混合算法。分为63x60.63又分为7x9.60分为3x4x5。
  • 混合基fft变换

    2013-12-05 15:38:58
    适用于LTE标准的混合机快速傅里叶变换,非常实用
  • DIT和DIF的2FFT算法

    万次阅读 2012-12-08 16:52:30
    根据课本上分析的DIT和DIF的步骤以及特点,写了两个DIF和DIT的2fft算法。 DIT和DIF,为了方便编程,对于前者将输入按倒位序重新排列,输出几位自然顺序排列;后者的话,输入为自然顺序,输出为倒位序。 难点...

    根据课本上分析的DIT和DIF的步骤以及特点,写了两个DIF和DIT的基2fft算法。


    DIT和DIF,为了方便编程,对于前者将输入按倒位序重新排列,输出几位自然顺序排列;后者的话,输入为自然顺序,输出为倒位序。

    难点就是倒位序的算法,以及FFT的循环算法。


    下面是最简单的序列长度为2的整数幂

    1.倒位序算法,两者都是通用的。

    (1)折断一半再拼上,再折断一半再拼上,知道列向量变成行向量

                    

    这就反应了倒位序的实质。

    我只能想到用循环:

    x=[1:8]';
    L=length(x);
    for i=1:1:log2(L)
          L=L/2;
          y=mat2cell(x,[L,L]);
         y=reshape(y,1,2);
         x=cell2mat(y);
    end

    询问了大牛,不用循环,直接用matlab的函数有:

    第一种方法:

    N = 3;
    a = (1:8)';
    y = permute(reshape(a,2*ones(1,N)),N:-1:1);
    y(:).'

    第二种方法:

    N=4;
    p(1:2^N,N:-1:1)=dec2bin(0:(2^N-1));
    re=bin2dec(char(p))+1;
    re
     
    第二种方法就是说,每一次折断一半再拼上,就想等于下标数循环右移

    即000=>000,001=>100,010=>010,011=>101....类似这样。



    (2)至于循环,就是利用相互之间的关系。

    p与第几级与k之间的关系,书上说,考虑到硬件的设备成本问题,发现每一级每两个相互乘积求和的下标与其他两者无关,故不需要在引进其他的硬件设备。

    x(k)=x(k)+x(k+b)e^(-j*2*π/N*p);

    x(k+b)=x(k)-x(k+b)e^(-j*2*π/N*p);


    在编程的时候注意x(k)需要重新赋给其他字母,不然会引起混淆错误,即

     mid=x(k)

     x(k)=mid+x(k+b)e^(-j*2*π/N*p);

     x(k+b)=mid-x(k+b)e^(-j*2*π/N*p)

    上述三句语句就是核心,其他的循环无非就是找各种k,p,n等等参数的关系!


    (3)最后强调下,一直调了好久的问题,复数和实数的转置(复数转置不可用x‘,一定要用x.')切记!

    展开全文
  • (二):FFT

    千次阅读 2013-08-05 17:31:34
    FFT: //首先思路上来说,是分治,递归: /*RECURSIVE-FFT(a) *n=a.length *if n==1 return a//一个数的DFT是自身 w_n=e^(2PIi/n) w=1 a_0=(a0,a2,....an-2) a_1=(a1,a3,....an-1) y_0=RECURISIVE-FFT(a_0)...
  • 2FFT时间抽取和频域抽取算法比较

    千次阅读 2013-10-30 15:30:09
    *FFT算法*/ #include "math.h" #include "stdio.h" struct compx { double real; double imag; } compx ; struct compx EE(struct compx b1,struct compx b2)//复数相乘 { struct compx b3; b3.real=b1.real*...
  • 快速傅里叶变换的2FFT算法的C++实现 2011-01-19 05:26 快速傅里叶变换的基本原理由于公式不好显示请读者参考其它文章或书籍,本文重点给出了时域抽取法FFT的C++实现。 下面是算法的流程图 倒序的流程图...
  • 确知离散信号互相关函数Rxy(m) = 确知离散信号互相关卷积定理:Rxy(m) = x(-m)*y(m) 时域卷积=频域相乘 设x(n)长度为N1,y(n)长度为N2 ...(2)用FFT计算x(n)和y(n)的N点DFT FFT[x(n)] = X(k) ...
  • Radix-2 FFT 网上的资源很多,但是Radix-4 FFT的资源很少,我只找到一个C++版本的,而且网上几乎没有 IFFT 的代码。 我实现了一个python版本的Radix-4,包括正变换,和反变换,方便理解和学习。 其中第一个版本的...
  • -4FFT算法 当混合基FFT算法中 时 即为-4FFT算法nk都为4进制数 个 点DFT 乘N个旋转因子 个 点DFT 乘N个旋转因子 个 点DFT 1) 的4点DFT 3) 的4点 DFT
  • 2频率抽取定点FFT算法

    热门讨论 2010-12-08 15:04:45
    自己做的2FFT的定点程序 内含定点fft同时还有浮点fft程序 对定点和浮点的fft结果进行了误差比较 希望能对大家有用
  • 傅里叶变换 ~ 2 时间抽取 FFT 算法

    千次阅读 2019-10-21 22:21:59
    文章目录1、2时间抽取FFT算法原理2、2时间抽取FFT算法流图2.1、示例1 ~ 4点的序列表示成两个2点的DFT2.2、示例2 ~ 8点的序列表示成两个2点的DFT2.3、实例演示32时间抽取FFT算法流图特点3.1、蝶形图的关系3.2...
  • 2与4时分FFT算法浅析及其比较

    万次阅读 多人点赞 2019-02-28 16:43:39
    FFT 算法的实质是把一长序列的 DFT 计算分割为较短序列的 DFT 计算,对于基2算法而言,是把序列每次一分为二,最后分割成两点 DFT,也可以采用别的分割法,每次一分为三,四,五等,就得到了基3,基4,基5等算法,...
  • 3.3 频率抽取2-FFT

    2017-06-10 12:11:39
    3.3 频率抽取2-FFT
  • 3.2 时间抽取2-FFT算法

    千次阅读 2017-06-10 12:09:53
    3.2 时间抽取2-FFT算法
  • FFT Algorithms

    2009-10-24 13:10:51
    -2,-3,-5混合基FFT算法介绍,作者Brian Gough。
  • DIT2点fft的编程实现

    千次阅读 2011-11-03 12:34:28
    找了很多资料,发现这个...下一次贴出4FFT编程实现。 蝶形公式: X(K) = X’(K) + X’(K+B)W PN , X(K+B) = X’(K) - X’(K+B) W PN 其中W PN= cos(2πP/N)- jsin(2πP/N)。 设 X(K+B) = XR(K+B) + jXI
  • FFT

    千次阅读 2017-05-22 14:03:19
    快速傅立叶变换(FFT) 函数FFT和IFFT用于快速傅立叶变换和逆变换。下面介绍这些函数。 ...函数的一种调用格式为 y=fft(x) ...其中,x是序列,y是序列的FFT,x...如果x长度是2的幂次方,函数fft执行高速-2FFT算法;
  • from scipy.fftpack import fft def my_fft(x, N): """ 实现FFT x: 离散序列 N: x的长度 """ # 反序排列 F = np.array(x) LH = int(N / 2) J = LH # print(J) M = N - 2 for I in range(1, M+1

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,915
精华内容 1,166
关键字:

基3fft