精华内容
参与话题
问答
  • FFT

    千次阅读 多人点赞 2018-04-11 10:04:56
    FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号...

    FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如
    果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。 
        虽然很多人都知道FFT是什么,可以用来做什么,怎么去做,但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。
        现在圈圈就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号,经过ADC采样之后,就变成了数字信号。采样定理告诉我们,采样频率要大于信号频率的两倍,这些我就不在此罗嗦了。

        采样得到的数字信号,就可以做FFT变换了。N个采样点,经过FFT之后,就可以得到N个点的FFT结果。为了方便进行FFT运算,通常N取2的整数次方。

        假设采样频率为Fs,信号频率F,采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值,就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢?假设原始信号的峰值为A,那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。而第一个点就是直流分量,它的模值就是直流分量的N倍。而每个点的相位呢,就是在该频率下的信号的相位。第一个点表示直流分量(即0Hz),而最后一个点N的再下一个点(实际上这个点是不存在的,这里是假设的第N+1个点,也可以看做是将第一个点分做两半分,另一半移到最后)则表示采样频率Fs,这中间被N-1个点平均分成N等份,每个点的频率依次增加。例如某点n所表示的频率为:Fn=(n-1)*Fs/N。由上面的公式可以看出,Fn所能分辨到频率为为Fs/N,如果采样频率Fs为1024Hz,采样点数为1024点,则可以分辨到1Hz。1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时间的信号并做FFT,则结果可以分析到0.5Hz。如果要提高频率分辨力,则必须增加采样点数,也即采样时间。频率分辨率和采样时间是倒数关系。
         假设FFT之后某点n用复数a+bi表示,那么这个复数的模就是An=根号a*a+b*b,相位就是Pn=atan2(b,a)。根据以上的结果,就可以计算出n点(n≠1,且n<=N/2)对应的信号的表达式为:An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号,是直流分量,幅度即为A1/N。
         由于FFT结果的对称性,通常我们只使用前半部分的结果,即小于采样频率一半的结果。
         总结:假设采样频率为Fs,采样点数为N,做FFT之后,某一点n(n从1开始)表示的频率为:Fn=(n-1)*Fs/N;该点的模值除以N/2就是对应该频率下的信号的幅度(对于直流信号是除以N);该点的相位即是对应该频率下的信号的相位。相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值,范围从-pi到pi。要精确到xHz,则需要采样长度为1/x秒的信号,并做FFT。要提高频率分辨率,就需要增加采样点数,这在一些实际的应用中是不现实的,需要在较短的时间内完成分析。解决这个问题的方法有频率细分法,比较简单的方法是采样比较短时间的信号,然后在后面补充一定数量的0,使其长度达到需要的点数,再做FFT,这在一定程度上能够提高频率分辨力。
    具体的频率细分法可参考相关文献。

    [附录:本测试数据使用的matlab程序]
    close all; %先关闭所有图片
    Adc=2;  %直流分量幅度
    A1=3;   %频率F1信号的幅度
    A2=1.5; %频率F2信号的幅度
    F1=50;  %信号1频率(Hz)
    F2=75;  %信号2频率(Hz)
    Fs=256; %采样频率(Hz)
    P1=-30; %信号1相位(度)
    P2=90;  %信号相位(度)
    N=256;  %采样点数
    t=[0:1/Fs:N/Fs]; %采样时刻

    %信号
    S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
    %显示原始信号
    plot(S);
    title('原始信号');

    figure;
    Y = fft(S,N); %做FFT变换
    Ayy = (abs(Y)); %取模
    plot(Ayy(1:N)); %显示原始的FFT模值结果
    title('FFT 模值');

    figure;
    Ayy=Ayy/(N/2);   %换算成实际的幅度
    Ayy(1)=Ayy(1)/2;
    F=([1:N]-1)*Fs/N; %换算成实际的频率值
    plot(F(1:N/2),Ayy(1:N/2));   %显示换算后的FFT模值结果
    title('幅度-频率曲线图');

    figure;
    Pyy=[1:N/2];
    for i="1:N/2"
    Pyy(i)=phase(Y(i)); %计算相位
    Pyy(i)=Pyy(i)*180/pi; %换算为角度
    end;
    plot(F(1:N/2),Pyy(1:N/2));   %显示相位图
    title('相位-频率曲线图');

    看完这个你就明白谐波分析了。
    展开全文
  • fft

    2018-02-22 20:44:23
    资料 http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform ...两道板子题 【模板】A*B Problem升级版(FFT快速傅里叶) #include&lt;cstdio...

    资料
    http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
    https://blog.csdn.net/iamzky/article/details/22712347
    两道板子题
    【模板】A*B Problem升级版(FFT快速傅里叶)

    #include<cstdio>
    #include<cmath>
    #include<iostream>
    const double PI=acos(-1);
    const int N=2*610000;
    using namespace std;
    struct complex{
        double x,y;
        complex(){}
        complex(double xx,double yx) {x=xx,y=yx;} 
        complex operator *(complex &a)const{
            return complex(x*a.x-y*a.y,x*a.y+y*a.x);
        }
        complex operator +(complex &a)const{
            return complex(x+a.x,y+a.y);
        }
        complex operator -(complex &a)const{
            return complex(x-a.x,y-a.y);
        }
    }a[N],b[N];
    int n,m,c[N],R[N];
    void fft(complex *a,int f){
        for(int i=0,j=0;i<n;i++){
            if(i<R[i]) swap(a[i],a[R[i]]);
            //for(int l=n>>1;(j^=1)<l;l>>=1);
        }
        for(int i=1;i<n;i<<=1) {complex w(cos(PI/i),f*sin(PI/i));
            for(int j=0;j<n;j+=(i<<1)){complex e(1,0);
                for(int k=0;k<i;k++,e=e*w){ 
                complex x=a[j+k],y=e*a[j+k+i];a[j+k]=x+y,a[j+k+i]=x-y;
                }   
            }
        }
        if(f==-1) for(int i=0;i<n;i++) a[i].x/=n;
    }
    char x[N],y[N];
    int main(){
        scanf("%d%s%s",&n,x,y);n--;
        for(int i=0;i<=n;i++) a[i].x=x[n-i]-'0',b[i].x=y[n-i]-'0';
        m=n*2;int L=0; 
        for(n=1;n<=m;n<<=1)L++;
        for(int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
        fft(a,1),fft(b,1);
        for(int i=0;i<=n;i++) a[i]=a[i]*b[i]; 
        fft(a,-1);
        for(int i=0;i<=m;i++) c[i]=(int)(a[i].x+0.1); 
        for(int i=0;i<=m;i++){
    
            if(c[i]>=10){
                c[i+1]+=c[i]/10;c[i]%=10;
                if(i==m) m++;
            }
        }
        while(m) if(c[m]) break; else m--;
        while(~m) printf("%d",c[m--]);
        return 0; 
    
    } 

    P3803 【模板】多项式乘法(FFT)

    #include<cstdio>
    #include<cmath>
    #include<iostream>
    const double PI=acos(-1);
    const int N=2201000;
    using namespace std;
    struct complex{
        double x,y;
        complex(){}
        complex(double xx,double yx) {x=xx,y=yx;} 
        complex operator *(complex &a)const{
            return complex(x*a.x-y*a.y,x*a.y+y*a.x);
        }
        complex operator +(complex &a)const{
            return complex(x+a.x,y+a.y);
        }
        complex operator -(complex &a)const{
            return complex(x-a.x,y-a.y);
        }
    }a[N],b[N];
    int n,m,c[N],R[N];
    void fft(complex *a,int f){
        for(int i=0,j=0;i<n;i++){
            if(i<R[i]) swap(a[i],a[R[i]]);
            //for(int l=n>>1;(j^=1)<l;l>>=1);
        }
        for(int i=1;i<n;i<<=1) {complex w(cos(PI/i),f*sin(PI/i));
            for(int j=0;j<n;j+=(i<<1)){complex e(1,0);
                for(int k=0;k<i;k++,e=e*w){ 
                complex x=a[j+k],y=e*a[j+k+i];a[j+k]=x+y,a[j+k+i]=x-y;
                }   
            }
        }
        if(f==-1) for(int i=0;i<n;i++) a[i].x/=n;
    }
    char x[N],y[N];
    int main(){
        scanf("%d%d",&n,&m);
        for(int i=0;i<=n;i++) scanf("%lf",&a[n-i].x);
        for(int j=0;j<=m;j++) scanf("%lf",&b[m-j].x);
        m=n+m+1;int L=0; 
        for(n=1;n<=m;n<<=1)L++;
        for(int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
        fft(a,1),fft(b,1);
        for(int i=0;i<=n;i++) a[i]=a[i]*b[i]; 
        fft(a,-1);
        for(int i=0;i<=m;i++) c[i]=(int)(a[i].x+0.1); 
        while(m) if(c[m]) break; else m--;
        while(~m) printf("%d ",c[m--]);
        return 0; 
    
    } 
    展开全文
  • 十分简明易懂的FFT(快速傅里叶变换)

    万次阅读 多人点赞 2018-08-07 11:38:56
    FFT有什么用 快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使...

    FFT前言

    快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称,简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。

    FFT(Fast Fourier Transformation) 是离散傅氏变换(DFT)的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。

    ——百度百科

    FFT(Fast Fourier Transformation),中文名快速傅里叶变换,用来加速多项式乘法
    朴素高精度乘法时间O(n2)O(n^2),但FFTFFTO(nlog2n)O(n\log_2 n)的时间解决
    FFTFFT名字逼格高,也难懂,其他教程写得让人看不太懂,于是自己随便写一下

    • 建议对复数三角函数相关知识有所耳闻 (不会也无所谓)

    • 下面难懂的点我会从网上盗


    多项式的系数表示法和点值表示法

    • FFTFFT其实是一个用O(nlog2n)O(n\log_2n)的时间将一个用系数表示的多项式转换成它的点值表示的算法

    • 多项式的系数表示和点值表示可以互相转换

    系数表示法

    一个n-1次n项多项式f(x)f(x)可以表示为f(x)=i=0n1aixif(x)=\sum^{n-1}_{i=0}a_ix^i
    也可以用每一项的系数来表示f(x)f(x),即f(x)={a0,a1,a2,...,an1}f(x)=\{a_0,a_1,a_2,...,a_{n-1} \}
    这就是系数表示法,也就是平时数学课上用的方法

    点值表示法

    • 把多项式放到平面直角坐标系里面,看成一个函数

    • nn个不同的xx代入,会得出nn个不同的yy,在坐标系内就是nn个不同的点

    • 那么这nn个点唯一确定该多项式,也就是有且仅有一个多项式满足k,f(xk)=yk∀k,f(x_k)=y_k

    • 理由很简单,把nn条式子联立起来成为一个有n条方程的n元方程组,每一项的系数都可以解出来

    那么f(x)f(x)还可以用f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),...,(xn1,f(xn1))}f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\}来表示
    这就是点值表示法


    高精度乘法下两种多项式表示法的区别

    对于两个用系数表示的多项式,我们把它们相乘
    设两个多项式分别为A(x),B(x)A(x),B(x)
    我们要枚举AA每一位的系数与BB每一位的系数相乘
    那么系数表示法做多项式乘法时间复杂度O(n2)O(n^2)

    但两个用点值表示的多项式相乘,只需要O(n)O(n)的时间

    什么意思呢?

    设两个点值多项式分别为f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),...,(xn1,f(xn1))}f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\}g(x)={(x0,g(x0)),(x1,g(x1)),(x2,g(x2)),...,(xn1,g(xn1))}g(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),...,(x_{n-1},g(x_{n-1}))\}
    设它们的乘积是h(x)h(x),那么h(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),...,(xn1,f(xn1)g(xn1))}h(x)=\{(x_0,f(x_0)·g(x_0)),(x_1,f(x_1)·g(x_1)),...,(x_{n-1},f(x_{n-1})·g(x_{n-1}))\}

    所以这里的时间复杂度只有一个枚举的O(n)O(n)

    • 突然感觉高精度乘法能O(n)O(n)暴艹一堆题?

    • 但是朴素的系数表示法转点值表示法的算法还是O(n2)O(n^2)的,逆操作类似

    • 朴素系数转点值的算法叫DFT(离散傅里叶变换),点值转系数叫IDFT(离散傅里叶逆变换)

    • 难道高精度乘法只能O(n2)O(n^2)了吗?


    DFT前置知识&技能

    复数

    毕竟高中有所以不多说

    我们把形如a+bi(a,b均为实数)的数称为复数,其中a称为实部,b称为虚部,i称为虚数单位。当虚部等于零时,这个复数可以视为实数;当z的虚部不等于零时,实部等于零时,常称z为纯虚数。复数域是实数域的代数闭包,也即任何复系数多项式在复数域中总有根。 复数是由意大利米兰学者卡当在十六世纪首次引入,经过达朗贝尔、棣莫弗、欧拉、高斯等人的工作,此概念逐渐为数学家所接受。
    ——百度百科

    初中数学老师会告诉你没有1\sqrt{-1},但仅限RR
    扩展至复数集CC,定义i2=1i^2=-1,一个复数zz可以表示为z=a+bi(a,bR)z=a+bi(a,b\in R)
    其中aa称为实部bb称为虚部ii称为虚数单位

    • 在复数集中就可以用ii表示负数的平方根,如7=7i\sqrt{-7}=\sqrt{7}i

    还可以把复数看成平面直角坐标系上的一个点,比如下面


    xx轴就是实数集中的坐标轴,yy轴就是虚数单位ii

    这个点(2,3)(2,3)表示的复数就是2+3i2+3i,或者想象它代表的向量(2,3)(2,3)
    其实我们还可以自己想象 (其实没有这种表达方式) 它可以表示为(13,θ)(\sqrt{13},\theta)
    一个复数zz定义为它到原点的距离,记为z=a2+b2|z|=\sqrt{a^2+b^2}
    一个复数z=a+biz=a+bi共轭复数abia-bi(虚部取反),记为z=abi\overline{z}=a-bi

    复数的运算

    复数不像点或向量,它和实数一样可以进行四则运算
    设两个复数分别为z1=a+bi,z2=c+diz_1=a+bi,z_2=c+di,那么
    z1+z2=(a+c)+(b+d)iz_1+z_2=(a+c)+(b+d)iz1z2=(acbd)+(ad+bc)iz_1z_2=(ac−bd)+(ad+bc)i

    复数相加也满足平行四边形法则


    这张是从网上盗的

    AB+AD=ACAB+AD=AC

    复数相乘还有一个值得注意的小性质
    (a1,θ1)(a2,θ2)=(a1a2,θ1+θ2)(a_1,\theta_1)*(a_2,\theta_2)=(a_1a_2,\theta_1+\theta_2)
    模长相乘,极角相加


    DFT(离散傅里叶变换)

    • 一定注意从这里开始所有的nn都默认为22的整数次幂

    对于任意系数多项式转点值,当然可以随便取任意nnxx值代入计算
    但是暴力计算xk0,xk1,...,xkn1(k[0,n))x_k^0,x_k^1,...,x_k^{n-1}(k\in[0,n))当然是O(n2)O(n^2)的时间
    其实可以代入一组神奇xx,代入以后不用做那么多的次方运算
    这些xx当然不是乱取的,而且取这些xx值应该就是 傅里叶 的主意了

    考虑一下,如果我们代入一些xx,使每个xx若干次方等于11,我们就不用做全部的次方运算了
    ±1±1是可以的,考虑虚数的话±i±i也可以,但只有这四个数远远不够

    • 傅里叶说:这个圆圈上面的点都可以做到

    以原点为圆心,画一个半径为11单位圆
    那么单位圆上所有的点都可以经过若干次次方得到11
    傅里叶说还要把它给nn等分了,比如此时n=8n=8

    橙色点即为n=8n=8时要取的点,从(1,0)(1,0)点开始,逆时针从00号开始标号,标到77
    记编号为kk的点代表的复数值为ωnk\omega_n^k,那么由模长相乘,极角相加可知(ωn1)k=ωnk(\omega_n^1)^k=\omega_n^k
    其中ωn1\omega_n^1称为nn次单位根,而且每一个ω\omega都可以求出 (我三角函数不好)
    ωnk=coskn2π+isinkn2π\omega_n^k=\cos{k\over n}2π+i\sin{k\over n} 2π

    或者说也可以这样解释这条式子


    注意sin2θ+cos2θ=1sin^2\theta+cos^2\theta=1什么的,就容易理解了

    那么ωn0,ωn1,...,ωnn1\omega^0_n,\omega^1_n,...,\omega^{n-1}_n即为我们要代入的x0,x1,...,xn1x_0,x_1,...,x_{n-1}


    单位根的一些性质

    FFTFFT的过程中需要用到ω\omega的一些性质

    ωnk=ω2n2k\omega^k_n=\omega^{2k}_{2n}

    • 它们表示的点(或向量)表示的复数是相同的

    • 证明

    • ωnk=coskn2π+isinkn2π=cos2k2n2π+isin2k2n2π=ω2n2k\omega^k_n=cos{k\over n}2π+isin{k\over n} 2π=cos{2k\over 2n}2π+isin{2k\over 2n} 2π=\omega^{2k}_{2n}

    ωnk+n2=ωnk\omega^{k+{n \over 2}}_n=-\omega_n^k

    • 它们表示的点关于原点对称,所表示的复数实部相反,所表示的向量等大反向

    • 证明

    • ωnn2=cosn2n2π+isinn2n2π=cosπ+isinπ=1\omega^{n\over 2}_n=cos{{n\over 2}\over n}2\pi+isin{{n\over 2}\over n}2\pi=cos\pi+isin\pi=-1

    • (这个东西和eix=cosx+isinxe^{ix}=cosx+isinxeiπ+1=0e^{i\pi}+1=0有点关系,我不会就不讲了

    ωn0=ωnn\omega^0_n=\omega^n_n

    • 都等于11,或1+0i1+0i

    FFT(快速傅里叶变换)

    虽然DFTDFT搞出来一堆很牛逼ω\omega作为代入多项式的xx
    但只是代入这类特殊xx值法的变换叫做DFTDFT而已,还是要代入单位根暴力计算

    • DFT还是暴力O(n2)O(n^2)

    DFTDFT可以分治来做,于是 FFT(快速傅里叶变换) 就出来了
    首先设一个多项式A(x)A(x)
    A(x)=i=0n1aixi=a0+a1x+a2x2+...+an1xn1A(x)=\sum^{n-1}_{i=0}a_ix^i=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}

    A(x)A(x)下标的奇偶性A(x)A(x)分成两半,右边再提一个xx

    A(x)=(a0+a2x2+...+an2xn2)+(a1x+a3x3+...+an1xn1)A(x)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})

    =(a0+a2x2+...+an2xn2)+x(a1+a3x2+...+an1xn2)=(a_0+a_2x^2+...+a_{n-2}x^{n-2})+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})

    两边好像非常相似,于是再设两个多项式A1(x),A2(x)A_1(x),A_2(x),令

    A1(x)=a0+a2x+a4x2+...+an2xn21A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{{n\over 2}-1}A2(x)=a1+a3x+a5x2+...+an1xn21A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{{n \over 2}-1}

    比较明显得出
    A(x)=A1(x2)+xA2(x2)A(x)=A_1(x^2)+xA_2(x^2)

    再设k&lt;n2k&lt;{n\over 2},把ωnk\omega^k_n作为xx代入A(x)A(x)(接下来几步变换要多想想单位根的性质)

    A(ωnk)=A1((ωnk)2)+ωnkA2((ωnk)2)A(\omega^k_n)=A_1((\omega^k_n)^2)+\omega^k_nA_2((\omega^k_n)^2)=A1(ωn2k)+ωnkA2(ωn2k)=A1(ωn2k)+ωnkA2(ωn2k)=A_1(\omega^{2k}_n)+\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})+\omega^k_nA_2(\omega^k_{n\over 2})

    那么对于A(ωnk+n2)A(\omega^{k+{n\over2}}_n)的话,代入ωnk+n2\omega^{k+{n \over 2}}_n
    A(ωnk+n2)=A1(ωn2k+n)+ωnk+n2A2(ωn2k+n)A(\omega^{k+{n\over 2}}_n)=A_1(\omega^{2k+n}_n)+\omega^{k+{n\over 2}}_nA_2(\omega^{2k+n}_n)=A1(ωn2kωnn)ωnkA2(ωn2kωnn)=A_1(\omega^{2k}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^n_n)=A1(ωn2k)ωnkA2(ωn2k)=A1(ωn2k)ωnkA2(ωn2k)=A_1(\omega^{2k}_n)-\omega^k_nA_2(\omega^{2k}_n)=A_1(\omega^k_{n\over2})-\omega^k_nA_2(\omega^k_{n\over2})

    • 发现了什么?

    A(ωnk)A(\omega^k_n)A(ωnk+n2)A(\omega^{k+{n\over2}}_n)两个多项式后面一坨东西只有符号不同
    就是说,如果已知A1(ωn2k)A_1(\omega^k_{n\over 2})A2(ωn2k)A_2(\omega^k_{n\over 2})的值,我们就可以同时知道A(ωnk)A(\omega^k_n)A(ωnk+n2)A(\omega^{k+{n\over2}}_n)的值
    现在我们就可以递归分治来搞FFTFFT

    每一次回溯时只扫当前前面一半的序列,即可得出后面一半序列的答案
    n==1n==1时只有一个常数项,直接returnreturn
    时间复杂度O(nlog2n)O(n\log_2n)


    IFFT(快速傅里叶逆变换)

    想一下,我们不仅要会FFTFFT,还要会IFFT(快速傅里叶逆变换)
    我们把两个多项式相乘 (也叫求卷积),做完两遍FFTFFT也知道了积的多项式的点值表示
    可我们平时用系数表示的多项式,点值表示没有意义

    • 怎么把点值表示的多项式快速转回系数表示法?

    • IDFTIDFT暴力O(n2)O(n^2)做?其实也可以用FFTFFTO(nlog2n)O(n\log_2n)的时间搞

    你有没有想过为什么傅里叶是把ωnk\omega^k_n作为xx代入而不是别的什么数?
    原因是有的但是有我也看不懂
    由于我是沙雕所以只用记住一个结论

    • 一个多项式在分治的过程中乘上单位根的共轭复数,分治完的每一项除以nn即为原多项式的每一项系数

    意思就是说FFTFFTIFFTIFFT可以一起搞


    朴素版FFT板子

    c++c++有自带的复数模板complexcomplex
    a.real()a.real()即表示复数aa的实部

    #include<complex>
    #define cp complex<double>
    
    void fft(cp *a,int n,int inv)//inv是取共轭复数的符号
    {
        if (n==1)return;
        int mid=n/2;
        static cp b[MAXN];
        fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1];
        fo(i,0,n-1)a[i]=b[i];
        fft(a,mid,inv),fft(a+mid,mid,inv);//分治
        fo(i,0,mid-1)
        {
            cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n));//inv取决是否取共轭复数
            b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid];
        }
        fo(i,0,n-1)a[i]=b[i];
    }
    

    两个多项式a,ba,b相乘再转系数多项式cc,通常只用打这么一小段

    	cp a[MAXN],b[MAXN];
    	int c[MAXN];
    	fft(a,n,1),fft(b,n,1);//1系数转点值
    	fo(i,0,n-1)a[i]*=b[i];
    	fft(a,n,-1);//-1点值转系数
    	fo(i,0,n-1)c[i]=(int)(a[i].real()/n+0.5);//注意精度
    

    很明显,FFT只能处理nn22的整数次幂的多项式
    所以在FFTFFT前一定要把nn调到22的次幂

    这个板子看着好像很优美,但是

    递归常数太大,要考虑优化…


    FFTの优化——迭代版FFT

    这个图也是盗的

    这个很容易发现点什么吧?

    • 每个位置分治后的最终位置为其二进制翻转后得到的位置

    这样的话我们可以先把原序列变换好,把每个数放在最终的位置上,再一步一步向上合并
    一句话就可以O(n)O(n)预处理第ii位最终的位置rev[i]rev[i]

    fo(i,0,n-1)rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
    

    至于蝴蝶变换它死了其实是我不会


    真·FFT板子

    void fft(cp *a,int n,int inv)
    {
        int bit=0;
        while ((1<<bit)<n)bit++;
        fo(i,0,n-1)
        {
            rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
            if (i<rev[i])swap(a[i],a[rev[i]]);//不加这条if会交换两次(就是没交换)
        }
        for (int mid=1;mid<n;mid*=2)//mid是准备合并序列的长度的二分之一
        {
        	cp temp(cos(pi/mid),inv*sin(pi/mid));//单位根,pi的系数2已经约掉了
            for (int i=0;i<n;i+=mid*2)//mid*2是准备合并序列的长度,i是合并到了哪一位
    		{
                cp omega(1,0);
                for (int j=0;j<mid;j++,omega*=temp)//只扫左半部分,得到右半部分的答案
                {
                    cp x=a[i+j],y=omega*a[i+j+mid];
                    a[i+j]=x+y,a[i+j+mid]=x-y;//这个就是蝴蝶变换什么的
                }
            }
        }
    }
    

    这个板子好像不是那么好背
    至少这个板子已经很优美


    FFT后记

    本人版权意识薄弱……

    NTTNTT我来了

    展开全文
  • <div><p>Using some known answer tests from another project I find that many of the 16 bit complex FFT functions are giving the wrong answer. I have modified the Ne10 conformance testing for 16 bit ...
  • 以前一直对MATLAB中fft()函数的使用一直存在疑惑,为什么要加一 些参数,并且如何确定这些参数,也查了许多资料,但很多都感觉只是 表面一说根本没有讲清其本质。但随着学习的推进,慢慢有所领悟,所 以打算把...
    以前一直对MATLAB中fft()函数的使用一直存在疑惑,为什么要加一
    些参数,并且如何确定这些参数,也查了许多资料,但很多都感觉只是
    表面一说根本没有讲清其本质。但随着学习的推进,慢慢有所领悟,所
    以打算把自己的一些所懂分享下,有什么问题也希望大家指正。
    本文主要先对DFT、FFT的一些概念进行介绍,然后通过MATLAB仿真
    进行fft()分析,从而解释上述参数。
    

    一、DFT与FFT
    首先是对DFT与FFT的一些概念上的介绍,其实FFT与DFT是等价的,他们实现的功能是一样的,只是FFT是DFT的算法优化,因为毕竟要用电脑来计算,DFT算的太慢了,就优化下也就成了FFT。所以此处我们对DFT与FFT的介绍是等价的。
    那么我们就来介绍DFT,它也被叫做离散傅里叶变换,其实它就是DFS离散傅里叶级数的时域频域主值序列,或者也是DTFT离散时间傅里叶变换的频域采样。至于DFS与DTFT相信大家也是明白,那么很多人好奇为什么还要DFT这玩意,主要呢,还是因为计算机,因为计算机不可能处理无限长的信号而DFS和DTFT要么时域无限长,要么频域无限长,所以就搞了个DFT。
    好的,那么我们就来直接看公式(这里我就假想大家都对其来源及原因都清楚了):
    在这里插入图片描述

    在这里插入图片描述
    注意看上图,这就是我们接下来的所有依据,X(k)是正变换,x(n)反变换,此处务必注意x(n)那个式子有个1/N,这对后面的理解很关键。
    好的那么接下来我们就可以利用MATLAB来进行分析了。

    二、MATLAB的FFT函数使用
    首先说明下什么情况下我们要用FFT,这是很简单的,但还是要说说:因为现实世界都是连续的信号,相信我们要分析它时单纯靠我们人算是很麻烦的,所有这些都是需要计算机进行计算,那么答案就有了,我们其实都是处理连续信号,也就是说我们都是想要FFT来分析连续信号。
    OK
    1、那么这样FFT第一个步骤就出现了,那就是A/D转换,也就是对连续信号采样。而这里我们就要确定一个参数采样频率Fs,说到采样,大家肯定会想到采样定理,其内容就是当对信号进行以采样频率为Fs>=2fc的采样时,信息不急丢失(fc为原始连续信号的最大截止频率),那么这样第一个参数Fs就这样确定了,其应该满足Fs>=2fc,这时采样周期为Ts=1/Fs。
    2、既然确定了采样频率Fs,那么我们就要想信号的时域应该怎么确定,换句话说,我们将采样好的信号给计算机处理时,应该给计算机多少个。这里就出现了第二个参数即序列个数L,好的,在这里由于第一个参数确定的Fs,相应的也是采样周期Ts,所以我们就可以得出时域信号的横坐标时间的变化范围,也就是
    t=(0:L-1)Ts
    那么这里就确定了N吗?不,我们这里还没确定,对于L的确立我们需要上面这个语句的帮助。那么究竟怎么确定呢,这里我们就要想到DFT的来源了:为什么DFT可以用来进行表示信号的频谱响应,因为其时域频域都是信号的相应时域频域的主值区间,所以说我们要用DTF可以表示这个信号时,应尽可能的包含这个信号的一个周期或整数倍周期也就是上述t的范围应该是信号的一个周期范围或整数倍周期范围,否则则会发生频谱泄露,至于频谱泄露其实就是出现了本来没有的频率分量。比如说,10Hz的cos(2π10t),本来只有一种频率分量f=10Hz,但分析结果却包含了与10Hz频率相近的其它频率分量。下述将会进行仿真分析。所以说此处L的确定应是LxTs=mT(T为信号的最小周期,m为正整数),常常m=1。那么或许有人问怎么平常我看到的都不是这样,L都是取和Fs一样的值(注意为了频率分辨率为1Hz,下文会讨论)。其实对于这个也有它的道理,不妨我们推到下:
    LxTs=L/Fs=1,也就是说其表达的时域t范围为1,那么平常我们的信号最小周期一般都小于1的,要是当信号最小周期大于1时其就是错误的啦,所以这个只是对大部分来说是有用的。
    那么说了这么多究竟怎么确定,要是按LxTs=mT确定会不会太麻烦,确实有时我们真的不能知道确定的信号最小周期,所以这里就放宽点,只要
    保证LxTs>=T
    就可以,当然如果可以取整数倍就最好了。
    OK
    可是上面说又会发生频率泄露,这可怎么办?那么这里就只能尽可能提高N-DFT的N点啦,这样能使尽可能更多的能量落到正确的频率点上,也就是说这样误差会更小,这样就能达到我们的目的。
    3、于是第三个参数N-DFT的N就出现了,根据上述2情况,我们知道N的提高有利于提高精确的,那么它还有其它约束吗?答案肯定有的,这里就需要些理论知识了:
    先来了解下频域采样定理:频域以N点为采样周期的采样,在时域就是原序列的以周期N的延拓
    那么这样也就是说上述2中的序列长度L必须小于或等于N,否则就会发生时域混叠,所以说L>=N这就是限制条件,常常我们取L=N,当发生频谱泄露比较严重时我们就可以考虑将L增大,这样就可以减小频谱泄露的的影响。

    综上,以上就是MATLAB相关的所有参数了,它们分别就采样频率Fs,原序列长度L,N-DFT变换的N,其主要确定关系总结如下:
    Fs>=2fc(fc为信号最高截止频率)
    LxTs>=T(T为信号的最小周期,若能取LxTs=T则应尽可能取,可以避免频谱泄露)
    N>=L(N常常应为2^m,因为FFT运算就是不断模2进行DFT,所以N为2的次方利于提高速度)
    对于其常取值,Fs=1024,N=1024,L=1024,注意满足上述条件就好

    说完了这些参数,我们就讲讲MATLAB的fft函数,其用法主要Y=fft(x,N),x为原信号序列,Y为DFT(也就是FFT)变换后的,但是大家会发现其变换后幅度变了而且变了很大,其实就是跟最初的那个式子有关,其频域幅值增大了N倍
    在这里插入图片描述
    并且发现其横轴也都是整数,所以这时也要对其横轴坐标进行变换,即要乘以相应的频率分辨率Fs/N,所以说上面为什么,常常也取Fs=L,N=L,这样为了频率分辨率为1Hz,对于这些我们需用以下MATLAB语句进行实现

    Y=fft(x,N);%N一定到大于信号序列x的长度,不过一般等于,决不能小于
    %因为若要小于时域就会发生混叠
    Y=abs(Y);
    Y=Y./N;
    %Y=Y.*2./N;Y(1)=Y(1)./2;若考虑将负频率的幅度折算到正频率时应这样处理
    %因为变换后有个N的乘积因子的影响,根据DFT公式可知,故消除其影响
    f=(0:N-1)*Fs./N;
    %频率正常化,因为变换后横坐标是每个点对应的个数,故应转化成实际的频率
    %由于频谱是对称的,且周期的故常常只画一半如下
    subplot(2,1,1);
    plot(f(1:N/2),Y(1:N./2));
    %当然也可以画平常我们见到的有对称的如下
    f=f-Fs./2;
    subplot(212);
    plot(f,fftshift(Y));
    

    上述代码也有相关解释,那么我也就对相关几个语句进行介绍解释:
    f=(0:N-1)*Fs./N;这是N-DFT变换后的横轴坐标,也就是真实频率值;
    Y=Y./N;这是对DFT变换后幅值处理,也就是要除于N才是真实的幅值
    %Y=Y.*2./N;Y(1)=Y(1)./2;而这个是鉴于无真实的负频率,其与正频率对称的原理,把负频率去掉,全变为正频率的幅值处理。

    下面语句是为了变成我们平常看到有正负频率的图像
    f=f-Fs./2;
    subplot(212);
    plot(f,fftshift(Y));

    三、仿真试验
    OK接下来就是仿真试验啦。
    我们直接考虑x(t)=cos(2π10t);其周期T=0.1s,fc=10Hz
    1、首先考虑Fs=512,L=512,N=512,满足上述三个关系吧,
    Fs=512>=2fc;LxTs=1s>=T(也是10T,故没频率泄露);N=512>=L=512;那么我们看下仿真结果是否对得到:
    在这里插入图片描述
    可以看到只有一个频率,可以完全表示,故与上述推理一样。代码如下

    Fs=512;%采样频率
    Ts=1/Fs;%采样周期
    N=512;%N—DFT
    L=512;%原信号序列长度
    t=(0:L-1).*Ts;%时域自变量
    x=cos(10*2*pi*t);%原信号
    subplot(311);
    plot(t,x);
    title("原信号")
    xlabel("t/s")
    grid on
    Y=fft(x,N);%fft变换
    Y=abs(Y)./N;%实际幅值变换
    f=(0:N-1)*Fs./N;%实际频率变换
    subplot(312);
    plot(f(1:N/2),Y(1:N./2));
    title("N-DFT变换幅频响应单边")
    xlabel("f/Hz")
    grid on
    f=f-Fs./2;%移位
    subplot(313);
    plot(f,fftshift(Y));%移位
    title("N-DFT变换幅频响应双边")
    xlabel("f/Hz")
    grid on
    

    2、那么我们再来考虑Fs=512,L=128,N=128,满足上述三个关系吧,
    Fs=512>=2fc;LxTs=0.15s>=T(不是T的整数倍,故有频率泄露);N=128>=L=128;那么我们看下仿真结果是否对得到(这次频谱用stem绘制更明显点):
    在这里插入图片描述
    可以看到出现了些不应该出现的频率,也验证了我们的推论,其代码如下:

    Fs=512;%采样频率
    Ts=1/Fs;%采样周期
    N=128;%N—DFT
    L=128;%原信号序列长度
    t=(0:L-1).*Ts;%时域自变量
    x=cos(10*2*pi*t);%原信号
    subplot(311);
    plot(t,x);
    title("原信号")
    xlabel("t/s")
    grid on
    Y=fft(x,N);%fft变换
    Y=abs(Y)./N;%实际幅值变换
    f=(0:N-1)*Fs./N;%实际频率变换
    subplot(312);
    stem(f(1:N/2),Y(1:N./2));
    title("N-DFT变换幅频响应单边")
    xlabel("f/Hz")
    grid on
    f=f-Fs./2;%移位
    subplot(313);
    stem(f,fftshift(Y));%移位
    title("N-DFT变换幅频响应双边")
    xlabel("f/Hz")
    grid on
    

    3、那么我们再来看下增大N点是否可以减小频谱泄露的影响,考虑Fs=512,L=128,N=1024,满足上述三个关系吧,
    Fs=512>=2fc;LxTs=0.15s>=T(不是T的整数倍,故有频率泄露);N=1024>=L=128;那么我们看下仿真结果是否对得到(这次频谱用stem绘制更明显点):
    在这里插入图片描述
    很明显可以看到与N=128时,其往所求正确频率处更靠近集中,故增大N可以减小频率响应的影响。代码如下:

    Fs=512;%采样频率
    Ts=1/Fs;%采样周期
    N=1024;%N—DFT
    L=128;%原信号序列长度
    t=(0:L-1).*Ts;%时域自变量
    x=cos(10*2*pi*t);%原信号
    subplot(311);
    plot(t,x);
    title("原信号")
    xlabel("t/s")
    grid on
    Y=fft(x,N);%fft变换
    Y=abs(Y)./N;%实际幅值变换
    f=(0:N-1)*Fs./N;%实际频率变换
    subplot(312);
    stem(f(1:N/2),Y(1:N./2));
    title("N-DFT变换幅频响应单边")
    xlabel("f/Hz")
    grid on
    f=f-Fs./2;%移位
    subplot(313);
    stem(f,fftshift(Y));%移位
    title("N-DFT变换幅频响应双边")
    xlabel("f/Hz")
    grid on
    

    4、最后看下它可不可ifft()变回去:
    在这里插入图片描述
    这是对于1的ifft只在后面加上下述代码:

    figure
    subplot(211)
    plot(t,x)
    title("原信号")
    xlabel("t/s")
    grid on
    Y=fft(x);
    xx=ifft(Y);
    subplot(212)
    plot(t,xx)
    title("ifft原信号")
    xlabel("t/s")
    grid on
    

    所以ifft变换时,应把一开始的N因子变回去。
    对于非周期信号,其原理差不多,关键是时域能否尽可能取到一个周期等,或者对于无限长的信号,只能采取截取出其中最主要的信息了。

    综上,差不多就这样了,可能多少有些问题,欢迎大家指正互相进步。

    展开全文
  • Matlab FFT参数设置研究

    千次阅读 2019-04-12 11:48:09
    近期要对一款高速ADC进行测试,用到Matlab的fft函数分析其动态性能,为了对Matlab 的fft有一个全方位立体的认识,对其参数进行了小实验,记录如下。 使用Matlab生成采样数据 clear; fs = 1000; ts = 1/fs; L = ...
  • FFT with fftw3

    2020-11-23 03:50:53
    <p>It is not an issue, but I would like to understand how works FFT in CImg, especially with fftw3. Since I am concerned about performance I've test how scales the time consumption with 3d images ...
  • 快速傅里叶变换(FFT

    千次阅读 2016-10-31 16:40:56
    快速傅里叶变换(FFT)说到FFT,我们就不得不介绍傅里叶同学、傅里叶级数、傅里叶变换、离散傅里叶变换。下面对他们这些玩意一一介绍。 1.傅里叶同学简介让·巴普蒂斯·约瑟夫·傅里叶(Baron Jean Baptiste Joseph...
  • Fix FFT overlap calculation

    2020-11-28 19:26:44
    <div><p>Gqrx displays the FFT overlap percentage in the FFT Settings pane. I noticed two errors in the calculation: <ol><li>The calculation doesn't take into account the fact that the FFT interval...
  • 如何正确理解FFT

    2020-11-04 07:39:10
    您可以通过周期性地收集大量的 ADC 输出转换采样来生成 FFT图。一般而言,ADC 厂商们将一种单音、满量程模拟输入信号用于其产品说明书的典型性能曲线。您从这些转换获得数据,然后绘制出一幅与图 1 相似的图。该图的...
  • 语音处理快速傅里叶变换xf = np.fft.rfft(xs,fft_size)/fft_size之后信号长度没有变为fft_size/2 +1 长度仍然是xs的长度,按理来说频率个数应为:np.linspace(0, sample_rate/2, fft_size/2+1) 但是这样就会报错: ...

空空如也

1 2 3 4 5 ... 20
收藏数 13,822
精华内容 5,528
关键字:

FFT