精华内容
参与话题
问答
  • FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号...

    转自:https://blog.csdn.net/zhaopeizhaopeipei/article/details/53908238

    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结果的对称性,通常我们只使用前半部分的结果,即小于采样频率一半的结果。

        好了,说了半天,看着公式也晕,下面圈圈以一个实际的信号来做说明。

        假设我们有一个信号,它含有2V的直流分量,频率为50Hz、相位为-30度、幅度为3V的交流信号,以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下:

    S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

        式中cos参数为弧度,所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样,总共采样256点。按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个点之间的间距就是1Hz,第n个点的频率就是n-1。我们的信号有3个频率:0Hz、50Hz、75Hz,应该分别在第1个点、第51个点、第76个点上出现峰值,其它各点应该接近0。实际情况如何呢?我们来看看FFT的结果的模值如图所示。




                          图1 FFT结果
        从图中我们可以看到,在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看:
    1点: 512+0i
    2点: -2.6195E-14 - 1.4162E-13i 
    3点: -2.8586E-14 - 1.1898E-13i

    50点:-6.2076E-13 - 2.1713E-12i
    51点:332.55 - 192i
    52点:-1.6707E-12 - 1.5241E-12i

    75点:-2.2199E-13 -1.0076E-12i
    76点:3.4315E-12 + 192i
    77点:-3.0263E-14 +7.5609E-13i
       
        很明显,1点、51点、76点的值都比较大,它附近的点值都很小,可以认为是0,即在那些频率点上的信号幅度为0。接着,我们来计算各点的幅度值。分别计算这三个点的模值,结果如下:
    1点: 512
    51点:384
    76点:192
        按照公式,可以计算出直流分量为:512/N=512/256=2;50Hz信号的幅度为:384/(N/2)=384/(256/2)=3;75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。可见,从频谱分析出来的幅度是正确的。
        然后再来计算相位信息。直流信号没有相位可言,不用管它。先计算50Hz信号的相位,atan2(-192, 332.55)=-0.5236,结果是弧度,换算为角度就是180*(-0.5236)/pi=-30.0001。再计算75Hz信号的相位,atan2(192, 3.4315E-12)=1.5708弧度,换算成角度就是180*1.5708/pi=90.0002。可见,相位也是对的。
    根据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-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我来了

    展开全文
  • FFT算法讲解——麻麻我终于会FFT了!

    万次阅读 多人点赞 2018-03-12 17:15:19
    FFT——快速傅里叶变换 这块不写东西空荡荡的,我决定还是把FFT的定义给贴上吧 FFT(Fast Fourier Transformation)是离散傅氏变换(DFT)的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚...

    FFT——快速傅里叶变换

    这块不写东西空荡荡的,我决定还是把FFT的定义给贴上吧

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

    这三段话其实一点用也没有

    FFT是干什么的

    FFT在算法竞赛中就有一个用途:加速多项式乘法(暴言)

    简单来说,形如 a0X0+a1X1+a2X2++anXna_0X^0+a_1X^1+a_2X^2+⋯+a_nX^n 的代数表达式叫做多项式,可以记作f(X)=a0X0+a1X1+a2X2++anXnf(X)=a_0X^0+a_1X^1+a_2X^2+⋯+a_nX^n,其中a0,a1,,ana_0,a_1,⋯,a_n叫做多项式的系数XX是一个不定元(就是不可以合并),不表示任何值,不定元在多项式中最大项的次数称作多项式的次数

    如果我们当前有两个多项式f(X),g(X)f(X),g(X),现在要把他们乘起来(求卷积),最朴素的做法就是

    i=02n1(j=0iajbij)xi\sum \limits _{i=0}^{2n-1} (\sum \limits _{j=0} ^{i} a_j*b_{i-j})*x^i

    这样的复杂度是Θ(n2)\varTheta(n^2)的,十分不美观,FFT就是要将这个过程优化为Θ(nlogn)\varTheta(n \log n)

    前置技能

    多项式

    见上文

    复数

    复数形如a+bia+bi,其中i=1i=\sqrt {-1}

    aa叫作复数的实部,bibi叫做复数的虚部

    复数(a1+b1i)(a2+b2i)(a_1+b_1i)*(a_2+b_2i)相乘的值,即a1a2b1b2+(a1b2+a2b1)ia_1a_2-b_1b_2+(a_1b_2+a_2b_1)i也是一个复数,同时我们也得到了复数的乘法法则

    复数c+dic+di可以用这种方式表示出来

    (c+di)

    复数乘法的在复平面中表现为辐角相加,模长相乘

    单位根

    复数ww满足wn=1w^n=1称作wwnn次单位根,下图包含了所有的88次单位根(图中圆的半径是1)

    8次单位根

    同样的,下图是所有的4次单位根

    4次单位根

    聪明的你也许已经发现了单位根的些许性质,即

    w2n2m=wnmw_{2n}^{2m}=w_{n}^{m}

    wnm=wnm+n2w_n^m=-w_n^{m+\frac{n}{2}}

    这两个要记住,一会很有用

    多项式的系数表达法

    我们有多项式f(X)=a0X0+a1X1+a2X2++anXnf(X)=a_0X^0+a_1X^1+a_2X^2+⋯+a_nX^n,令a=(a0,a1,a2,...,an)\vec{a}=(a_0,a_1,a_2,...,a_n),则称A(X)A(X)为多项式f(X)f(X)的系数表示法

    在系数表示法下,计算多项式乘法是Θ(n2)\varTheta (n^2)

    多项式的点值表达法

    任取n+1n+1互不相同S={p1,p2,...,pn+1}S=\{p_1,p_2,...,p_{n+1}\},对f(X)f(X)分别求值得到f(p1),f(p2),...,f(pn+1)f(p_1),f(p_2),...,f(p_{n+1}),那么称A(X)={(p1,f(p1)),(p2,f(p2)),...,(pn+1,f(pn+1))}A(X)=\{(p_1,f(p_1)),(p_2,f(p_2)),...,(p_{n+1},f(p_{n+1}))\}为多项式f(X)f(X)SS下的点值表示法

    可以把多项式想象成一个nn次函数,点值表示法就是取SS下每一个横坐标时对应的点,因为nn次函数可以由n+1n+1个点确定下来(可以将每一个点列一个nn次方程),所以nn维点值与nn维系数一一对应

    更重要的一点,点值表示法下的乘法运算获得了简化

    两个多项式PP,QQ分别取点(x,y1)(x,y_1)(x,y2)(x,y_2)PQP*Q就会取到点(x,y1y2)(x,y_1*y_2)

    C=PQC=P*Q,因为C(X)=P(X)Q(X)C(X)=P(X)*Q(X),所以C(x)=P(x)Q(x)C(x)=P(x)*Q(x),即C(x)=y1y2C(x)=y_1*y_2

    所以在点值表示法下,计算多项式乘法是Θ(n)\varTheta(n)

    FFT的具体过程

    FFT就是将系数表示法转化成点值表示法相乘,再由点值表示法转化为系数表示法的过程,第一个过程叫做求值(DFT),第二个过程叫做插值(IDFT)

    求值

    还记得我们之前提到的单位根吗?再回顾一下:

    w2n2m=wnmw_{2n}^{2m}=w_{n}^{m}

    wnm=wnm+n2w_n^m=-w_n^{m+\frac{n}{2}}

    想要求出一个多项式的点值表示法,需要选出n+1n+1个数分别带入到多项式里面,带入一个数的复杂度是Θ(n)\varTheta(n)的,那么总复杂度就是Θ(n2)\varTheta(n^2)的,因为单位根有上面两个优美的性质,所以我们尝试可以取nn次单位根组成SS,看看能不能加速我们的运算

    A0(X)A_0(X)A(X)A(X)偶次项的和,设A1(X)A_1(X)A(X)A(X)奇次项的和,即

    A0(X)=a0x0+a2x1+...+an1xn/2A_0(X)=a_0x^0+a_2x^1+...+a_{n-1}x^{n/2}

    A1(X)=a1x0+a3x1+...+an2xn/2A_1(X)=a_1x^0+a_3x^1+...+a_{n-2}x^{n/2}

    因为A(wnm)=a0wn0+a1wnm+a2wn2m+a3wn3m+...+an1wn(n1)m+anwnnmA(w_n^m)=a_0w_n^0+a_1w_n^m+a_2w_n^{2m}+a_3w_n^{3m}+...+a_{n-1}w_n^{(n-1)*m}+a_nw_n^{nm}

    所以有

    ​在这里插入图片描述

    在这里插入图片描述

    也就是说,只要有了A0(X)A_0(X)A1(X)A_1(X)的点值表示,就能在Θ(n)\varTheta(n)时间算出A(X)A(X)的点值表示,对于当前层确定的位置ii,就可以用下一层的两个值更新当前的值,我们称这个操作为“蝴蝶变换”

    因为这个过程一定要求每层都可以分成两大小相等的部分,所以多项式最高次项一定是2p2^p(pNp\in N)次方,如果不是的话,直接在最高次项补零就可以啦

    于是我们有了递归的写法:

    void FFT(Complex* a,int len){
        if(len==1) return;
        Complex* a0=new Complex[len/2];
        Complex* a1=new Complex[len/2];
        for(int i=0;i<len;i+=2){
            a0[i/2]=a[i];
            a1[i/2]=a[i+1];
        }
        FFT(a0,len/2);FFT(a1,len/2);
        Complex wn(cos(2*Pi/len),sin(2*Pi/len));
        Complex w(1,0);
        for(int i=0;i<(len/2);i++){
            a[i]=a0[i]+w*a1[i];
            a[i+len/2]=a0[i]-w*a1[i];
            w=w*wn;
        }
        return;
    }
    

    但递归版的FFT常数巨大,实现起来比较复杂,于是又有了迭代的写法

    这里写图片描述

    重新考虑下递归FFT的过程,在第ii次求解中,我们将所有元素二进制ii位为00的放在了左面,ii位为11的放在了右面,事实上,每个元素最终到的是他二进制颠倒过来的位置

    这里写图片描述

    再拿一张别人的图

    这里写图片描述

    迭代写法:

    inline void DFT(Complex a[]){
        for(int i=0;i<len;i++)///pos[i]代表反转后的位置
            if(i<pos[i])
                swap(a[i],a[pos[i]]);
        for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){///len代表多项式最高次项
            ///第一层i是枚举合并到了哪一层。
            Complex wm(cos(2.0*pi/i),sin(2.0*pi/i));
            for(int j=0;j<len;j+=i){///第二层j是枚举合并区间。
                Complex w(1,0);
                for(int k=j;k<j+mid;k++,w=w*wm){///第三层k是枚举区间内的下标。
                    Complex l=a[k],r=w*a[k+mid];
                    a[k]=l+r;a[k+mid]=l-r;
                }
            }
        }
        return;
    }
    

    插值

    有人证出来插值只要将所有wnmw_n^m换成wnm+n2w_n^{m+\frac{n}{2}},也就是所有的虚部取相反数,再将最终结果除以lenlen,就是插值的过程了

    究竟为什么?我觉得人生一定要有点遗憾可以参考这里,我就不多说了

    支持插值的迭代写法:

    const double DFT=2.0,IDFT=-2.0;///进行求值,第二个参数传DFT,插值传IDFT
    inline void FFT(Complex a[],double mode){
        for(int i=0;i<len;i++)
            if(i<pos[i])
                swap(a[i],a[pos[i]]);
        for(int i=2,mid=1;i<=len;i<<=1,mid<<=1){
            Complex wm(cos(2.0*pi/i),sin(mode*pi/i));
            for(int j=0;j<len;j+=i){
                Complex w(1,0);
                for(int k=j;k<j+mid;k++,w=w*wm){
                    Complex l=a[k],r=w*a[k+mid];
                    a[k]=l+r;a[k+mid]=l-r;
                }
            }
        }
        if(mode==IDFT)
            for(int i=0;i<len;i++)
                a[i].x/=len;
        return;
    }
    

    现在,你已经会写FFT了!

    生成函数

    AAaia_i个价值为AiA_i的物品,小BBbib_i个价值为AiA_i的物品,求用两个组成价值为cic_i的方案数

    生成函数可以解决上面的这个问题,构造两个多项式,第ii项表示价值为ii的物品有多少个,对两个人分别构造,乘在一起的多项式就代表所有的方案数了

    NTT——快速数论变换

    当FFT需要取模时,复数的存在就特别尴尬了,有没有什么东西可以代替单位根呢?

    原根性质:假设gg​是素数pp​的一个原根,则g1,g2,...,g(p1)g^1,g^2,...,g^{(p-1)}​在模pp​意义下两两不同,且结果恰好为11​~p1p-1​

    wnng(p1)1(mod p)w_n^n≡g^{(p-1)}≡1(mod\ p)(中间那个是费马小定理)

    所以可以把g(p1)/ng^{(p-1)/n}看成wn1w_n^1的等价。

    但这种质数必须是NTT质数(费马质数),即(p1)(p-1)有超过序列长度的2的正整数幂因子的质数

    具体证明还是看这里

    参考资料:
    Pick’s blog
    Fenghr的博客
    Miskcoo’s Space

    展开全文
  • FFT在matlab中的使用方法

    万次阅读 多人点赞 2019-10-04 16:29:30
    FFT在matlab中的用法 一、FFT的物理意义 ​ FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析...

    FFT在matlab中的用法

    一、FFT的物理意义

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

    二、计算序列的FFT变换

    求序列{2,3,3,2}的DFT变换。

    >> N=4;
    >> n=0:N-1;
    >> xn=[2 3 3 2];
    >> xk=fft(xn)
    

    运算结果如下:

    xk =
    
          10.0000 + 0.0000i  -1.0000 - 1.0000i   0.0000 + 0.0000i  -1.0000 + 1.0000i
    

    带入公式检验:
    X[k]=n=0N1X[n]WNnk X[k]=\sum_{n=0}^{N-1}X[n]W_N^{nk}

    X[0]=2W40+3W40+3W40+2W40=10 X[0]=2W_4^{0}+3W_4^{0}+3W_4^{0}+2W_4^{0}=10

    X[1]=2W40+3W41+3W42+2W43=1i X[1]=2W_4^{0}+3W_4^{1}+3W_4^{2}+2W_4^{3}=-1-i

    X[2]=2W40+3W42+3W44+2W46=0 X[2]=2W_4^{0}+3W_4^{2}+3W_4^{4}+2W_4^{6}=0

    X[3]=2W40+3W43+3W46+2W49=1+i X[3]=2W_4^{0}+3W_4^{3}+3W_4^{6}+2W_4^{9}=-1+i

    公式运算结果与matlab仿真结果一致。

    ​ xk与xn的维数相同,共有4个元素。xk的第一个数对应于直流分量,即频率值为0,在傅里叶级数的叠加中,它仅仅影响全部波形相对于数轴整体向上或是向下而不改变波的形状。

    三、计算三角信号的FFT变换

    1、分析采样点N的数量对FFT结果的影响

    ​ 做n个点的FFT,表示在时域上对原来的信号取了n个点来做频谱分析,n点FFT变换的结果仍为n个点。

    信号一:y=0.5sin(2pi20t)+2sin(2pi40t)

    ​ 可以很容易的看出信号的由频率为20和40的两个分量组成,其中最高频率为40,根据奈科斯特定律,

    在这里我将抽样信号fs设定为100HZ,采样点N分别设定为128和2048。

    ​ N=128

    >> fs=100;
    >> N=128;
    >> n=0:N-1;
    >> t=n/fs;
    >> y=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
    >> x=fft(y,N);
    >> m=abs(x);
    >> f=n*fs/N;
    >> subplot(2,2,1),plot(f,m);
    >> xlabel('频率/Hz');
    >> ylabel('振幅');title('N=128');
    >> grid on;
    >> subplot(2,2,2),plot(f(1:N/2),m(1:N/2));
    >> grid on;
    

    ​ N=2048

    >> fs=100;
    >> N=2048;
    >> n=0:N-1;
    >> t=n/fs;
    >> y=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
    >> x=fft(y,N);
    >> m=abs(x);
    >> f=n*fs/N;
    >> subplot(2,2,3),plot(f,m);
    >> xlabel('频率/Hz');
    >> ylabel('振幅');title('N=2048');
    >> grid on;
    >> subplot(2,2,4),plot(f(1:N/2),m(1:N/2));
    >> grid on;
    

    分析结果:

    ​ 这里的fs是采样频率。而我们通常只关心0-pi中的频谱,因为根据奈科斯特定律,只有f=fs/2范围内的信号才是被采样到的有效信号。那么,在w的范围内,得到的频谱肯定是关于n/2对称的。

    ​ 用100Hz做了128个点FFT之后,因为频率分辨率是100/128(25/32)Hz,如果原来的信号在2Hz或者1Hz有分量,你在频谱上是看不见的,这就表示你越想频谱画得逼真,就必须取越多的点数来做FFT,n就越大,你在时域上就必须取更长的信号样本来做分析。但是无论如何,由于离散采样的原理,你不可能完全准确地画出原来连续时间信号的真实频谱,只能无限接近(就是n无限大的时候),这个就叫做频率泄露。在采样频率fs不变得情况下,频率泄漏可以通过取更多的点来改善,也可以通过做FFT前加窗来改善,这就是另外一个话题了。

    ​ 第k个点的实际频率的计算为f(k)=k*(fs/n)

    运行结果:

    在这里插入图片描述

    2、对FFT输出结果的探究

    ​ 从上一节“分析采样点N的数量对FFT结果的影响”我们可以看出不管采样点数目是多少,其频率为20的分量和频率为40的分量的振幅之比依然是1:4。这是为什么呢?我们知道FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值,就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢?

    信号二:y=2+0.5sin(2pi20t)+2sin(2pi40t)

    我以128HZ的采样频率对其采样,采样点N=128。我们的信号有3个频率分量0HZ、20HZ、40HZ。

    >> fs=128;
    >> N=128;
    >> n=0:N-1;
    >> t=n/fs;
    >> y=2+0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
    >> x=fft(y,N)
    

    输出结果:

    x =
    
       1.0e+02 *
    
      Columns 1 through 7
    
       2.5600 + 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i
    
      Columns 8 through 14
    
      -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i
    
      Columns 15 through 21
    
       0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   -0.0000 - 0.3200i
    
      Columns 22 through 28
    
       0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i
    
      Columns 29 through 35
    
       0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i
    
      Columns 36 through 42
    
       0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 - 1.2800i   0.0000 - 0.0000i
    
      Columns 43 through 49
    
       0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i
    
      Columns 50 through 56
    
       0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i
    
      Columns 57 through 63
    
      -0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i
    
      Columns 64 through 70
    
       0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i
    
      Columns 71 through 77
    
       0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
    
      Columns 78 through 84
    
       0.0000 - 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i
    
      Columns 85 through 91
    
      -0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 1.2800i   0.0000 - 0.0000i   0.0000 - 0.0000i
    
      Columns 92 through 98
    
      -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i
    
      Columns 99 through 105
    
       0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i
    
      Columns 106 through 112
    
       0.0000 - 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.3200i  -0.0000 - 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i
    
      Columns 113 through 119
    
      -0.0000 + 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 + 0.0000i
    
      Columns 120 through 126
    
      -0.0000 + 0.0000i  -0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i   0.0000 + 0.0000i  -0.0000 - 0.0000i   0.0000 - 0.0000i
    
      Columns 127 through 128
    
       0.0000 + 0.0000i  -0.0000 + 0.0000i
    

    可以发现:

    第1个点:2.5600 + 0.0000i 模值256

    第2个点:-0.0000 - 0.0000i

    第20个点:-0.0000 + 0.0000i

    第21个点: -0.0000 - 0.3200i 模值32

    第22个点:-0.0000 + 0.0000i

    第40个点:0.0000 + 0.0000i

    第41个点:-0.0000 - 1.2800i 模值128

    第42个点:0.0000 + 0.0000i

    可以发现只有在第1、21、41点处才有幅值,画出其频谱图观察。

    >> m=abs(x);
    >> f=n*fs/N;
    >> plot(f,m);
    >> xlabel('频率/Hz');
    >> ylabel('振幅');title('N=128');
    >> grid on;
    

    运行结果:

    在这里插入图片描述

    ​ 可以看出频率0HZ、20HZ、40HZ对应点的幅值分别为256、36、128。这跟原始信号的幅度有什么关系呢?

    结论:

    ​ 假设原始信号的峰值为A,那么FFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。 而第一个点就是直流分量,它的模值就是直流分量的N倍。因此我们可以计算出原始信号的幅值。

    根据公式:直流分量的真实幅值为256/128=2,正确

    ​ 20HZ分量的真实幅值为32/(128/2)=0.5,正确

    ​ 40HZ分量的真实幅值为128/(128/2)=2,正确

    3、信号补零对FFT的影响

    ​ 在前面我们数据个数总是2得整数次幂,当数据个数不是2的整数次幂时会自动补零,以保持数据个数为2的整数次幂,那么补零会对FFT产生什么样的影响呢?

    信号一:y=0.5sin(2pi20t)+2sin(2pi40t)

    (1)数据个数N=32,FFT所用的采样点数NFFT=32;

    >> fs=100;
    >> N=32;
    >> NFFT=32;
    >> n=0:N-1;
    >> t=n/fs;
    >> y=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
    >> x=fft(y,NFFT);
    >> m=abs(x);
    >> f=(0:NFFT-1)*fs/NFFT;
    >> subplot(2,2,1),plot(f(1:NFFT/2),m(1:NFFT/2)*2/NFFT);
    >> xlabel('频率/Hz');ylabel('振幅');
    >> title('N=32 NFFT=32');grid on;
    

    (2)数据个数N=32,FFT所用的采样点数NFFT=128;

    >> N=32;
    >> NFFT=128;
    >> n=0:N-1;
    >> t=n/fs;
    >> y=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
    >> x=fft(y,NFFT);
    >> m=abs(x);
    >> f=(0:NFFT-1)*fs/NFFT;
    >> subplot(2,2,2),plot(f(1:NFFT/2),m(1:NFFT/2)*2/NFFT);
    >> xlabel('频率/Hz');ylabel('振幅');
    >> title('N=32 NFFT=128');grid on;
    

    (3)数据个数N=300,FFT所用的采样点数NFFT=256;

    >> N=300;
    >> NFFT=256;
    >> n=0:N-1;
    >> t=n/fs;
    >> y=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
    >> x=fft(y,NFFT);
    >> m=abs(x);
    >> f=(0:NFFT-1)*fs/NFFT;
    >> subplot(2,2,3),plot(f(1:NFFT/2),m(1:NFFT/2)*2/NFFT);
    >> xlabel('频率/Hz');ylabel('振幅');
    >> title('N=300 NFFT=256');grid on;
    

    (4)数据个数N=140,FFT所用的采样点数NFFT=512。

    >> N=140;
    >> NFFT=512;
    >> n=0:N-1;
    >> t=n/fs;
    >> y=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
    >> x=fft(y,NFFT);
    >> m=abs(x);
    >> f=(0:NFFT-1)*fs/NFFT;
    >> subplot(2,2,4),plot(f(1:NFFT/2),m(1:NFFT/2)*2/NFFT);
    >> xlabel('频率/Hz');ylabel('振幅');
    >> title('N=140 NFFT=512');grid on;
    

    运行结果:

    在这里插入图片描述
    结论:

    (1)当采样点数量为32时,由于数量比较少,频率分辨率比较低,画出的频谱图不是很逼真,但没有由于添零而导致的其他频率成分。

    (2)由于采样点的数目比信号数据个数多,所以信号必须在时间域内加零。致使振幅谱中出现很多其他成分,这是加零造成的。其振幅由于加了多个零而明显减小。

    (3)当数据个数超过频率采样点的个数时,FFT会自动切断。

    (4)信号也在时域上补零了,但是由于数据的个数比较多,FFT振幅谱看起来比较逼真。

    总结:

    ​ 如果采样数据比较少,运用FFT变换时,不能很准确的分辨出其中的频率成分。添加零后可增加频谱中的数据个数,谱的密度增高了,但仍不能分辨其中的频率成分,还会导致振幅的下降,即谱的分辨率没有提高。只有数据点数足够多时才能分辨其中的频率成分。

    展开全文
  • 浅谈 FFT (终于懂一点了~~)

    万次阅读 多人点赞 2018-08-19 11:38:09
    FFT(离散傅氏变换的快速算法) FFT(Fast Fourier Transformation)是离散傅氏变换(DFT)的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。 ...
  • 本算法采用C语言编写,可自定义抽样序列,实现了对初始序列进行快速傅里叶变换(FFT)和离散傅里叶变换(DFT),通过运行窗口将运行结果输出,并且输出FFT和DFT运行时间,方便读者分析这两个算法用时效率。
  • 离散傅里叶变换-DFT(FFT基础)

    万次阅读 多人点赞 2018-08-13 11:50:03
    本文是从最基础的知识开始讲解,力求用最通俗易懂的文字将问题将的通俗易懂,大神勿喷,多多指教啊,虽然说是从零学习FFT,但是基本的数学知识还是要有的,sin,cos,等。 &nbsp;&nbsp;&nbsp;&nbsp...
  • FFT&DFT简单总结 前言 相信大家都知道大(chou)名(ming)鼎(zhao)鼎(zhu)的FFT(fake_fake_true)(fast_fast_tle),并且都有过被它各种玄学操作虐待的经历(大佬请绕路),那么希望这篇详细的FFT简介...
  • FFT与DFT的关系及应用

    万次阅读 2016-09-12 16:14:49
    看了一篇讲理解离散傅立叶变换(二.... 首先,FFT(快速傅里叶变换)是一种实现DFT(离散傅里叶变换)的快速算法,是利用复数形式的离散傅里叶变换来计算实数形式的离散傅里叶变换)。matlab中的fft函数是实
  • 对于FFT和DFT的理解

    2020-03-30 12:08:58
    FFT是最重要,也是最难懂的。简单说下原理: FFT(快速傅里叶变换)是DFT(离散傅里叶变换)的改进算法,其将DFT的N^2 步运算减少至 ( N/2 )log2(N)步。 先来讲讲DFT的原理 离散傅里叶变换(DFT)是傅里叶变换...
  • FFT详解

    万次阅读 多人点赞 2017-03-31 10:28:34
    一直想学FFT,当时由于数学基础太差,导致啥都学不懂。请教了机房里的几位学长大神,结果还是没太明白。因此下定决心写一篇关于“FFT”的文章,一篇起码我能看得懂的“FFT”。
  • FFT

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

    2011-06-10 19:15:00
    FFT是离散傅立叶变换的快速算法,可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的,但是如果变换到频域之后,就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外,FFT可以将一个信号...
  • 以前一直对MATLAB中fft()函数的使用一直存在疑惑,为什么要加一 些参数,并且如何确定这些参数,也查了许多资料,但很多都感觉只是 表面一说根本没有讲清其本质。但随着学习的推进,慢慢有所领悟,所 以打算把...
  • 【STM32】使用STM32提供的DSP库进行FFT(附详细代码)

    万次阅读 多人点赞 2019-05-12 15:42:41
    最近,因为项目需要在STM32F103系列处理器上,对采集的音频信号进行FFT运算,然而STM32F103毕竟不是STM32F4系列的处理器,对于一般的FFT运算程序还是比较缓慢的。 幸亏官方提供了针对FFT的官方库,但是去官网找了...
  • 数字信号处理公式变程序(一)——DFT、FFT

    万次阅读 多人点赞 2015-07-03 09:15:26
    陆续会更新,包括FFT、巴特沃斯滤波器(高低带通、高低带阻)、数据差值(线性、sinc、三次样条*)、数据压缩(等距、平均、峰值检测)和模仿matlab的STFT功能(spectrogram函数三维绘图)。 注:可能会有不足...
  • FFT模板以及简单的FFT入门题总结

    千次阅读 2018-05-11 14:35:52
    FFT 前言:感谢M哥,Pls,Zjq等等一系列大神的指导,让我大概了解了FFT的基本原理以及能做的事情,下面是FFT的一些简单入门题(我也只会这些题,在这里总结一下。 由于目前题量较少,对于我目前的理解,FFT主要解决...
  • java里的fft2和fft:)

    千次阅读 2005-04-26 09:38:00
    fft是改网上c代码的,fft2是自己写的,谁能不能给一个求任意长度(即不需要将输入补成2的指数倍)的fft程序 public class FftList{ //提供fft的输出 double[] fr; double[] fi; } public class Fft2List{ //提供...
  • FFT后的物理意义

    万次阅读 多人点赞 2016-11-23 17:40:45
    FFT(Fast Fourier Transform,快速傅立叶变换)是离散傅立叶变换的快速算法,也是我们在数字信号处理技术中经常会提到的一个概念。在大学的理工科课程中,在完成高等数学的课程后,数字信号处理一般会作为通信电子...
  • 《C》C语言实现FFT算法

    万次阅读 多人点赞 2019-01-29 12:31:39
    一、什么是FFT? DFT虽好,但是其计算的次数太多,不利于大数据量的计算,FFT是DFT的快速算法,可以节省大量的计算时间,快速傅里叶变换(FFT)是一种能在O(nlogn)的时间内将一个多项式转换成它的点值表示的算法。 ...
  • 说明:本文适合信号处理方面有一定的基础的人阅读,能够理解什么时候傅里叶级数和傅里叶变换,能够理解他们的核心思想以及基本原理,能够理解到底什么是“频率域”,能够从频率的角度分析信号。...
  • Java版本的FFT和Inverse FFT

    千次阅读 2009-05-04 22:24:00
    最近有些朋友在一些项目中需要用到Java版本的FFT和Inverse FFT,玄机逸士在网上找到了一个版本,供大家参考,现抄录如下:(原文地址:http://www.cs.princeton.edu/introcs/97data/FFT.java.html)/***************...
  • FFT(最详细最通俗的入门手册)

    万次阅读 多人点赞 2017-07-24 20:17:14
    声明首先,我需要声明,本文是在转载的基础上稍微修饰的,经过原创作者...我学习 FFT 是一个比较慢的过程。 期间反反复复。 我写这篇博文只是一个非常浅显的理解。同时也可以帮助初学者在学习FFT的时候。有所偏重。避免
  • 我所理解的快速傅里叶变换(FFT

    万次阅读 多人点赞 2016-10-22 18:01:58
    首先FFT是离散傅立叶变换(DFT)的快速算法,那么说到FFT,我们自然要先讲清楚傅立叶变换。先来看看傅立叶变换是从哪里来的? 傅立叶是一位法国数学家和物理学家的名字,英语原名是Jean Baptiste Joseph Fourier(1768...
  • FFT算法

    千次阅读 2012-03-20 07:58:19
    FFT算法 分类: 信号处理2008-04-19 21:03 2582人阅读 评论(1) 收藏 举报 这几天,我一直在看FFT算法,下面分享一下我这几天学到的东西    1。直接计算离散傅立叶变换具有n^2的复杂度,而cooley...
  • 频谱图的横轴表示的是 频率, 纵轴表示的是振幅 #coding=gbk import numpy as np import pandas as pd import matplotlib.pyplot as plt #依据快速傅里叶算法得到信号的频域 ... fft_size = 8...
  • FFT移植心得

    千次阅读 2008-07-18 22:09:00
    花了2个晚上, 总算把FFT从MP3编码器中提取出来并在keil上跑通, 主要还是对keil和C51不熟, 走了不少弯路 1 开始怎么都跑不过去, 总是报MAIN.C(55): error C249: 'DATA': SEGMENT TOO LARGE突然想起不是大...
  • matlab信号频谱分析FFT详解

    万次阅读 多人点赞 2019-06-12 22:27:09
    前言 做OFDM通信少不了频谱分析,基带信号DA后的频谱,以及基带数字上变频后的DA信号都要频谱分析。我觉得其实做任何工程都是这样,先规定实施方案,然后仿真成功...matlab使用FFT函数分析信号频谱 一般我使用的F...

空空如也

1 2 3 4 5 ... 20
收藏数 40,358
精华内容 16,143
关键字:

FFT