精华内容
下载资源
问答
  • 提出了分数傅里叶变换计算全息图(FRTCGH).采用分数傅里叶变换的快速算法和罗曼Ⅲ型迂回位相编码方法设计并制作了一个物体不同分数阶的分数傅里叶变换计算全息图,用罗曼Ⅰ型分数傅里叶变换光学系统再现原物体,得到了...
  • 提出了双重分数傅里叶变换计算全息。在这种方法中,将两个图像的信息分别经不同阶的分数傅里叶变换后,记录在同一张分数傅里叶变换计算全息图上。它需要两个特定的分数傅里叶变换系统才能再现出所记录的图像信息,利用...
  • 用MATLAB制作二元傅里叶变换计算全息图.pdf
  • FFTW ( the Faster Fourier Transform in the West) 是一个快速计算离散傅里叶变换的标准C语言程序集,其由MIT的M.Frigo 和S. Johnson 开发。可计算一维或多维实和复数据以及任意规模的DFT。 FFTW 还包含对共享和...
  • 从取样定理出发,对角谱衍射公式及逆运算的快速傅里叶变换(FFT)进行了研究。基于研究结果,讨论数字全息研究中物光通过一个光学系统到达CCD探测器的波面重建问题,提出物光场的逆运算追迹重建及像空间波面重建两种方法...
  • 来源:...然而,使用快速傅里叶变换,时间复杂度可以降低到 O(N logN loglogN)。假设我们要计算以下两个 N 位数字的乘积:a = (aN-1aN-2…a1a0)10 = aN-1x

    来源:http://www.cnblogs.com/skyivben/archive/2008/07/23/1248413.html

    我们知道,两个 N 位数字的整数的乘法,如果使用常规的算法,时间复杂度是 O(N2)。然而,使用快速傅里叶变换,时间复杂度可以降低到 O(N logN loglogN)。

    假设我们要计算以下两个 N 位数字的乘积:

    a = (aN-1aN-2…a1a0)10 = aN-1x10N-1 + aN-2x10N-2 + … + a1x101 + a0x100

    b = (bN-1bN-2…b1b0)10 = bN-1x10N-1 + bN-2x10N-2 + … + b1x101 + b0x100

    将上面两个式子相乘,得到以下公式 (共 2N - 1 项):

    c = a x b = c2N-2x102N-2 + c2N-3x102N-3 + … + c1x101 + c0x100

    非常容易验证,上式中的 ck ( 0 ≤ k ≤ 2N-2 ) 满足以下公式:

    ck = a0xbk + a1xbk-1 + … + ak-2xb2 + ak-1xb1
    + akxb0 + ak+1xb-1 + … + aN-2xb-(N-2-k) + aN-1xb-(N-1-k)

    上式共有 N 项,ai 和 bj 的下标 i 和 j 满足 i + j = k。若不满足 0 ≤ i, j ≤ N-1 时,则令 ai = bj = 0。

    我们以两个 3 ( N = 3 ) 位数 a = 678 和 b = 432 的乘积来说明以上过程吧。

    a = (678)10 = 6x102 + 7x101 + 8x100

    b = (432)10 = 4x102 + 3x101 + 2x100

    由此:
    c0 = a0xb0 + a1xb-1 + a2xb-2 = 8x2 + 7x0 + 6x0 = 16 + 0 + 0 = 16
    c1 = a0xb1 + a1xb0 + a2xb-1 = 8x3 + 7x2 + 6x0 = 24 + 14 + 0 = 38
    c2 = a0xb2 + a1xb1 + a2xb0 = 8x4 + 7x3 +6x2 = 32 + 21 + 12 = 65
    c3 = a0xb3 + a1xb2 + a2xb1 = 8x0 + 7x4 + 6x3 = 0 + 28 + 18 = 46
    c4 = a0xb4 + a1xb3 + a2xb2 = 8x0 + 7x0 + 6x4 = 0 + 0 + 24 = 24

    最后:

    c = a x b = 104xc4 + 103xc3 + 102xc2 + 101xc1 + 100xc0
    = 10000x24 + 1000x46 + 100x65 + 10x38 + 1x16
    = 292896

    如果按以上方法计算大整数的乘法,时间复杂度是 O(N2)。

    但是,我们注意到,向量 {ck} 是向量 {ai} 和向量 {bj} 的卷积。根据卷积定理,向量卷积的离散傅里叶变换是向量离散傅里叶变换的乘积。于是,我们可以按照以下步骤来计算大整数乘法:

    分别求出向量 {ai} 和向量 {bj} 的离散傅里叶变换 {Ai} 和 {Bj}。
    将 {Ai} 和 {Bj} 逐项相乘得到向量 {Ck}。
    对 {Ck} 求离散傅里叶逆变换,得到的向量 {ck} 就是向量 {ai} 和向量 {bj} 的卷积。
    对的向量 {ck} 进行适当的进位就得到了大整数 a 和 b 的乘积 c。
    对于复数向量 { xN-1, …, x1, x0 },离散傅里叶变换公式为:

    离散傅里叶逆变换公式为:

    注意到离散傅里叶逆变换除了指数的符号相反以及结果需要乘以归一化因子 1/N 外,与离散傅里叶变换是相同的。所以计算离散傅里叶变换的程序稍做修改也可以用于计算逆变换。

    在我们的例子中,乘积 c = 292896,共 6 位数字,N 需要扩展到 23 = 8。那么,向量 {ai} 和向量 {bj} 如下所示:

    { a7, a6, a5, a4, a3, a2, a1, a0 } = { 0, 0, 0, 0, 0, 6, 7, 8 }

    { b7, b6, b5, b4, b3, b2, b1, b0 } = { 0, 0, 0, 0, 0, 4, 3, 2 }

    为了求出以上向量的离散傅里叶变换,我们令

    ω = e-2πi/N = e-2πi/8 = e-πi/4 = cos(-π/4) + i sin(-π/4) = √2 / 2 - i √2 / 2 ≈ 0.7-0.7i

    为了方便计算,我们预先求出 ω 的各次方,如下:

    ω8 = ω0 = e0 = 1

    ω9 = ω1 = e-πi/4 = cos(-π/4) + i sin(-π/4) ≈ 0.7-0.7i

    ω10 = ω2 = e-πi/2 = cos(-π/2) + i sin(-π/2) = -i

    ω11 = ω3 = e-3πi/4 = cos(-3π/4) + i sin(-3π/4) ≈ -0.7-0.7i

    ω12 = ω4 = e-πi = cos(-π) + i sin(-π) = -1

    ω13 = ω5 = e-5πi/4 = cos(-5π/4) + i sin(-5π/4) ≈ -0.7+0.7i

    ω14 = ω6 = e-3πi/2 = cos(-3π/2) + i sin(-3π/2) = i

    ω15 = ω7 = e-7πi/4 = cos(-7π/4) + i sin(-7π/4) ≈ 0.7+0.7i

    注意到当 n > 2 时,an = 0,于是:
    A0 = a0xω0x0 + a1xω1x0 + a2xω2x0 = 8xω0 + 7xω0 + 6xω0 = 8x1 + 7x1 + 6x1 = 21
    A1 = a0xω0x1 + a1xω1x1 + a2xω2x1 = 8xω0 + 7xω1 + 6xω2 ≈ 8x1 + 7x(0.7 - 0.7i) + 6x(-i) = 12.9-10.9i
    A2 = a0xω0x2 + a1xω1x2 + a2xω2x2 = 8xω0 + 7xω2 + 6xω4 = 8x1 + 7x(-i) + 6x(-1) = 2-7i
    A3 = a0xω0x3 + a1xω1x3 + a2xω2x3 = 8xω0 + 7xω3 + 6xω6 ≈ 8x1 + 7x(-0.7 - 0.7i) + 6xi = 3.1+1.1i

    A4 = a0xω0x4 + a1xω1x4 + a2xω2x4 = 8xω0 + 7xω4 + 6xω8 = 8x1 + 7x(-1) + 6x1 = 7
    A5 = a0xω0x5 + a1xω1x5 + a2xω2x5 = 8xω0 + 7xω5 + 6xω10 ≈ 8x1 + 7x(-0.7 + 0.7i) + 6x(-i) = 3.1-1.1i

    A6 = a0xω0x6 + a1xω1x6 + a2xω2x6 = 8xω0 + 7xω6 + 6xω12 = 8x1 + 7xi + 6x(-1) = 2+7i
    A7 = a0xω0x7 + a1xω1x7 + a2xω2x7 = 8xω0 + 7xω7 + 6xω14 ≈ 8x1 + 7x(0.7 + 0.7i) + 6xi = 12.9+10.9i
    同样,当 n > 2 时,bn = 0,于是:
    B0 = b0xω0x0 + b1xω1x0 + b2xω2x0 = 2xω0 + 3xω0 + 4xω0 = 2x1 + 3x1 + 4x1 = 9
    B1 = b0xω0x1 + b1xω1x1 + b2xω2x1 = 2xω0 + 3xω1 + 4xω2 ≈ 2x1 + 3x(0.7 - 0.7i) + 4x(-i) = 4.1-6.1i
    B2 = b0xω0x2 + b1xω1x2 + b2xω2x2 = 2xω0 + 3xω2 + 4xω4 = 2x1 + 3x(-i) + 4x(-1) = -2-3i
    B3 = b0xω0x3 + b1xω1x3 + b2xω2x3 = 2xω0 + 3xω3 + 4xω6 ≈ 2x1 + 3x(-0.7 - 0.7i) + 4xi = -0.1+1.9i

    B4 = b0xω0x4 + b1xω1x4 + b2xω2x4 = 2xω0 + 3xω4 + 4xω8 = 2x1 + 3x(-1) + 4x1 = 3
    B5 = b0xω0x5 + b1xω1x5 + b2xω2x5 = 2xω0 + 3xω5 + 4xω10 ≈ 2x1 + 3x(-0.7 + 0.7i) + 4x(-i) = -0.1-1.9i

    B6 = b0xω0x6 + b1xω1x6 + b2xω2x6 = 2xω0 + 3xω6 + 4xω12 = 2x1 + 3xi + 4x(-1) = -2+3i
    B7 = b0xω0x7 + b1xω1x7 + b2xω2x7 = 2xω0 + 3xω7 + 4xω14 ≈ 2x1 + 3x(0.7 + 0.7i) + 4xi = 4.1+6.1i

    这样,向量 {ai} 和向量 {bj} 的离散傅里叶变换 {Ai} 和 {Bj} 如下所示:

    { A7, A6, A5, A4, A3, A2, A1, A0 } = { 12.9+10.9i, 2+7i, 3.1-1.1i, 7, 3.1+1.1i, 2-7i, 12.9-10.9i, 21 }

    { B7, B6, B5, B4, B3, B2, B1, B0 } = { 4.1+6.1i, -2+3i, -0.1-1.9i, 3, -0.1+1.9i, -2-3i, 4.1-6.1i, 9 }

    现在,将她们逐项相乘得到向量 {Ck},即 { C7, C6, C5, C4, C3, C2, C1, C0 }

    = { -13.6+123.4i, -25-8i, -2.4-5.8i, 21, -2.4+5.8i, -25+8i, -13.6-123.4i, 189 }

    为了求出向量 {Ck} 的离散傅里叶逆变换,我们令

    ω = e2πi/N = e2πi/8 = eπi/4 = cos(π/4) + i sin(π/4) = √2 / 2 + i √2 / 2 ≈ 0.7+0.7i

    为了方便计算,我们预先求出 ω 的各次方(注意 ωk+8 = ωk),如下:

    ω0 = e0 = 1

    ω1 = eπi/4 = cos(π/4) + i sin(π/4) ≈ 0.7+0.7i

    ω2 = eπi/2 = cos(π/2) + i sin(π/2) = i

    ω3 = e3πi/4 = cos(3π/4) + i sin(3π/4) ≈ -0.7+0.7i

    ω4 = eπi = cos(π) + i sin(π) = -1

    ω5 = e5πi/4 = cos(5π/4) + i sin(5π/4) ≈ -0.7-0.7i

    ω6 = e3πi/2 = cos(3π/2) + i sin(3π/2) = -i

    ω7 = e7πi/4 = cos(7π/4) + i sin(7π/4) ≈ 0.7-0.7i

    于是:
    c0 = (1/N) x ( C0xω0x0 + C1xω1x0 + C2xω2x0 + C3xω3x0
    + C4xω4x0 + C5xω5x0 + C6xω6x0 + C7xω7x0 )
    = (1/8) x ( 189xω0 + (-13.6-123.4i)xω0 + (-25+8i)xω0 + (-2.4+5.8i)xω0
    + 21xω0 + (-2.4-5.8i)xω0 + (-25-8i)xω0 + (-13.6+123.4i)xω0 )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x1 + (-25+8i)x1 + (-2.4+5.8i)x1
    + 21x1 + (-2.4-5.8i)x1 + (-25-8i)x1 + (-13.6+123.4i)x1 )
    = 0.125 x 128 = 16
    c1 = (1/N) x ( 8xc1 = C0xω0x1 + C1xω1x1 + C2xω2x1 + C3xω3x1
    + C4xω4x1 + C5xω5x1 + C6xω6x1 + C7xω7x1 )
    = (1/8) x ( 189xω0 + ( -13.6-123.4i)xω1 + (-25+8i)xω2 + (-2.4+5.8i)xω3
    + 21xω4 + (-2.4-5.8i)xω5 + (-25-8i)xω6 + (-13.6+123.4i)xω7 )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(0.7+0.7i) + (-25+8i)x(i) + (-2.4+5.8i)x(-0.7+0.7i)
    + 21x(-1) + (-2.4-5.8i)x(-0.7-0.7i) + (-25-8i)x(-i) + (-13.6+123.4i)x(0.7-0.7i) )
    = 0.125 x 300.96 = 37.62 ≈ 38

    c2 = (1/N) x ( C0xω0x2 + C1xω1x2 + C2xω2x2 + C3xω3x2
    + C4xω4x2 + C5xω5x2 + C6xω6x2 + C7xω7x2 )
    = (1/8) x ( 189xω0 + (-13.6-123.4i)xω2 + (-25+8i)xω4 + (-2.4+5.8i)xω6
    + 21xω8 + (-2.4-5.8i)xω10 + (-25-8i)xω12 + (-13.6+123.4i)xω14 )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x(i) + (-25+8i)x(-1) + (-2.4+5.8i)x(-i)
    + 21x1 + (-2.4-5.8i)x(i) + (-25-8i)x(-1) + (-13.6+123.4i)x(-i) )
    ≈ 0.125 x 518.4 = 64.8 ≈ 65
    c3 = (1/N) x ( C0xω0x3 + C1xω1x3 + C2xω2x3 + C3xω3x3
    + C4xω4x3 + C5xω5x3 + C6xω6x3 + C7xω7x3 )
    = (1/8) x ( 189xω0 + (-13.6-123.4i)xω3 + (-25+8i)xω6 + (-2.4+5.8i)xω9
    + 21xω12 + (-2.4-5.8i)xω15 + (-25-8i)xω18 + (-13.6+123.4i)xω21 )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(-0.7+0.7i) + (-25+8i)x(-i) + (-2.4+5.8i)x(0.7+0.7i)
    + 21x(-1) + (-2.4-5.8i)x(0.7-0.7i) + (-25-8i)x(i) + (-13.6+123.4i)x(-0.7-0.7i) )
    = 0.125 x 364.32 = 45.54 ≈ 46

    c4 = (1/N) x ( C0xω0x4 + C1xω1x4 + C2xω2x4 + C3xω3x4
    + C4xω4x4 + C5xω5x4 + C6xω6x4 + C7xω7x4 )
    = (1/8) x ( 189xω0 + (-13.6-123.4i)xω4 + (-25+8i)xω8 + (-2.4+5.8i)xω12
    + 21xω16 + (-2.4-5.8i)xω20 + (-25-8i)xω24 + (-13.6+123.4i)xω28 )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x(-1) + (-25+8i)x1 + (-2.4+5.8i)x(-1)
    + 21x1 + (-2.4-5.8i)x(-1) + (-25-8i)x1 + (-13.6+123.4i)x(-1) )
    = 0.125 x 192 = 24

    c5 = (1/N) x ( C0xω0x5 + C1xω1x5 + C2xω2x5 + C3xω3x5
    + C4xω4x5 + C5xω5x5 + C6xω6x5 + C7xω7x5 )
    = (1/8) x ( 189xω0 + (-13.6-123.4i)xω5 + (-25+8i)xω10 + (-2.4+5.8i)xω15
    + 21xω20 + (-2.4-5.8i)xω25 + (-25-8i)xω30 + (-13.6+123.4i)xω35 )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(-0.7-0.7i) + (-25+8i)x(i) + (-2.4+5.8i)x(0.7-0.7i)
    + 21x(-1) + (-2.4-5.8i)x(0.7+0.7i) + (-25-8i)x(-i) + (-13.6+123.4i)x(-0.7+0.7i) )
    = 0.125 x 3.04 = 0.38 ≈ 0
    c6 = (1/N) x ( C0xω0x6 + C1xω1x6 + C2xω2x6 + C3xω3x6
    + C4xω4x6 + C5xω5x6 + C6xω6x6 + C7xω7x6 )
    = (1/8) x ( 189xω0 + (-13.6-123.4i)xω6 + (-25+8i)xω12 + (-2.4+5.8i)xω18
    + 21xω24 + (-2.4-5.8i)xω30 + (-25-8i)xω36 + (-13.6+123.4i)xω42 )
    = 0.125 x ( 189x1 + (-13.6-123.4i)x(-i) + (-25+8i)x(-1) + (-2.4+5.8i)x(i)
    + 21x1 + (-2.4-5.8i)x(-i) + (-25-8i)x(-1) + (-13.6+123.4i)x(i) )
    = 0.125 x 1.6 = 0.2 ≈ 0
    c7 = (1/N) x ( C0xω0x7 + C1xω1x7 + C2xω2x7 + C3xω3x7
    + C4xω4x7 + C5xω5x7 + C6xω6x7 + C7xω7x7 )
    = (1/8) x ( 189xω0 + (-13.6-123.4i)xω7 + (-25+8i)xω14 + (-2.4+5.8i)xω21
    + 21xω28 + (-2.4-5.8i)xω35 + (-25-8i)xω42 + (-13.6+123.4i)xω49 )
    ≈ 0.125 x ( 189x1 + (-13.6-123.4i)x(0.7-0.7i) + (-25+8i)x(-i) + (-2.4+5.8i)x(-0.7-0.7i)
    + 21x(-1) + (-2.4-5.8i)x(-0.7+0.7i) + (-25-8i)x(i) + (-13.6+123.4i)x(0.7+0.7i) )
    = 0.125 x 3.68 = 0.46 ≈ 0
    这样,我们就使用离散傅里叶变换和逆变换计算出了向量 {ai} 和向量 {bj} 的卷积向量 {ck},如下所示:
    { c7, c6, c5, c4, c3, c2, c1, c0 } = { 0, 0, 0, 0, 24, 46, 65, 38, 16 }
    这和我们在前面直接使用向量 {ai} 和向量 {bj} 来计算卷积的结果是一样的。

    但是,这个算法的时间复杂度还是 O(N2)。我们绕了这么一大圈,不是白费劲了吗?

    现在就到了关键时刻,关键在于:直接进行离散傅里叶变换的计算复杂度是 O(N2)。快速傅里叶变换可以计算出与直接计算相同的结果,但只需要 O(N logN) 的计算复杂度。 N logN 和 N2 之间的差别是巨大的。例如,当 N = 106 时,在一个每秒运算百万次的计算机上,粗略地说,它们之间就是占用 30 秒 CPU 时间和两星期 CPU 时间的差别。
    快速傅里叶变换的要点如下:一个界长为 N 的离散傅里叶变换可以重新写成两个界长各为 N/2 的离散傅里叶变换之和。其中一个变换由原来 N 个点中的偶数点构成,另一个变换由奇数点构成。这个过程可以递归地进行下去,直到我们将全部数据细分为界长为 1 的变换。什么是界长为 1 的傅里叶变换呢?它正是把一个输入值复制成它的一个输出值的恒等运算。要实现以上算法,最容易的情况是原始的 N 为 2 的整幂次项,如果数据集的界长不是 2 的幂次时,则可添上一些零值,直到 2 的下一幂次。在这个算法中,每递归一次需 N 阶运算,共需要 log N 次递归,所以快速傅里叶变换算法的时间复杂度是 O(N logN)。

    由于快速傅里叶变换是采用了浮点运算,因此我们需要足够的精度,以使在出现舍入误差时,结果中每个组成部分的准确整数值仍是可辨认的。长度为 N 的 B 进制数可产生大到 B2N 阶的卷积分量。我们知道,双精度浮点数的尾数是 53 个二进位,所以:

    2 x log2B + log2N + 几个 x log2log2N < 53

    上式中左边最后一项是为了快速傅里叶变换的舍入误差。

    所以,为了能够计算尽量大的整数,一般 B 不会取得太大。在计算机程序中经常使用 256 进制进行运算。但是如果经常需要将计算结果和十进制互相转换,则往往使用 100 进制进行运算。

    关于快速傅里叶变换以及卷积定理的更深入的知识,请参阅文末的参考文献。这一篇随笔主要是讲述相关的原理,在下一篇随笔中,我将给出一个使用快速傅里叶变换进行任意精度的算术运算的 C# 程序。

    顺便说一句,我在准备正文的例题的时候,是使用 google 计算器来进行复杂的复数运算的。发现她非常好用。以计算 c2 为例, 只要将要计算的表达式复制到 goole 搜索栏,然后按回车,就能得到计算结果:
    (189 x 1) + (((-13.6) - (123.4 * i)) x i) + (((-25) + (8 * i)) x (-1)) + (((-2.4) + (5.8 * i)) x (-i)) + (21 x 1) + (((-2.4) - (5.8 * i)) x i) + (((-25) - (8 * i)) x (-1)) + (((-13.6) + (123.4 * i)) x (-i)) = 518.4 - 1.77635684 × 10-15 i
    Google 计算器详情

    找不到和您的查询 “189x1 + (-13.6-123.4i)x(i) + (-25+8i)x(-1) + (-2.4+5.8i)x(-i) + 21x1 + (-2.4-5.8i)x(i) + (-25-8i)x(-1) + (-13.6+123.4i)x(-i)” 相符的网页。
    参考文献:
    Multiplication algorithm
    Convolution
    Convolution theorem
    Discrete Fourier transform
    Fast Fourier transform

    展开全文
  • 计算傅里叶变换 傅里叶变换 鉴于这种想法,任何信号,当然任何周期性信号,都可以由一系列正弦曲线组成,我们将开始从级数(Series)的概念转向连续信号的概念。我们将要讨论一些东西,这些东西可以让我们知道图像...

    目录

    傅里叶变换

    计算傅里叶变换


    傅里叶变换

    鉴于这种想法,任何信号,当然任何周期性信号,都可以由一系列正弦曲线组成,我们将开始从级数(Series)的概念转向连续信号的概念。我们将要讨论一些东西,这些东西可以让我们知道图像中任意给定频率的功率有多大。也就是说,我们想把图像从某样东西上变换过来,该东西是时间的函数,或者仅仅是空间的函数,从而知道它的频率是多少。这种变换被称为什么?它叫做傅里叶变换(Fourier Transform)

    现在我们做了傅里叶级数(Fourier series),然后我们进入傅里叶变换(Fourier Transform),再然后我们将进行离散傅里叶变换(Discrete Fourier Transform, DFT),这将使我们得到那个图像的东西。

    所以我们想了解频率,频率通常被写成\omega(Omega),我们的信号。我们想要参数化,我们想要变换信号,而不是x,在空间中是通过\omega(Omega)。我们通过所谓的傅里叶变换(Fourier Transform)来实现这一点,傅里叶变换会拿一些空间信号,

    并且它将为我们提供一些特性,这些特性可以给每个给定\omega(Omega)的相位\varphi 和 幅度A进行编码。

    好的,所以我们需要给定每个 \omega (Omega)的相位\varphi 和 幅度A。特别是每个\omega (Omega)从0到无穷大,实际上它是 从负无穷大到无穷大(-\infty+\infty)。这个想法是f(\omega)保存了正弦曲线相应的幅度A 和 相位\varphi。因此,这个f 必须拥有幅度A 和 相位\varphi

    f 如何保留振幅和相位?复数的把戏(Complex number trick!)。注意,我没有说过一个单独的数字。因为请记住,f(x)只是F的值,所以它是一个数字,f(\omega)将是幅度和相位,我们该怎么做?非常好,实际上,大写字母F实际上是一个复数。现在,希望您能记住关于复数的一些知识。如果没有,我会给你一个非常简短的回顾。所以基本上, f(\omega)由两部分组成:一个实部(Real Part),另一个是虚部(Imaginary Part)

    所以请记住:A+Bi ,在这种情况下,  f(\omega)的实部加上虚部,也许你还记得复数的大小就是两个元素之和的平方根。通常是\sqrt{A^2+B^2},这里是 \sqrt{R^2+I^2}

    R代表实部,I是虚部,相位之间的关系将以这种方式编写,相位\varphi是 虚部 比 实部 的反正切。

    这里遗漏的点是:实部归因于 偶函数(Even),虚部归因于 奇函数(Odd),

    所有这意味着:如果我在这里有一个函数,是我的cosign函数,这是关于原点对称的。

    正弦曲线,我要把它们拧紧,因为它必须到这里,

    然后在那里,

    再然后,所以它必须是0(zero),在这里,这是我的正弦曲线,因为当那个东西高时,它必须降到0,奇怪的是,Sin (- x) 等于 -Sin (x) ,而 Cos (- x) 等于 Cos (x)。

    所以最后的这个部分(如图),

    实部是余弦部分(Cosin),

    虚部是正弦部分(Sine)。

    计算傅里叶变换

    当我们进行傅里叶变换时,我们所做的只是计算基集。我会告诉你我的意思。看到这个丑陋,丑陋的整体吗?

    好吧,我得到的是,我已经得到了一个函数的正弦值,另一个函数的正弦值。,

    两个不同的频率,a 和 b

    我将这一点从负无穷大加到正无穷大,

    我声称如果 a 不等于 b,它等于0。

    为什么这是正确的呢?有猜测吗? 让我们直观地思考一下,如果有一个正弦曲线,

    现在我必须做一个不同的频率,这很容易,

    天啊,真漂亮!哈哈? 在本质上,在某个点上,当这个是正的这个是相同的正值,我会试着找到,希望有一些是正的,

    我必须补上一个点。对,让我们假装红色的东西在这里下来就像那样 一点点。

    当红色的点是正数时,黄色的点是负数,

    所以这两个数的乘积会相互抵消,这是一种波浪的方式表示整个积分是0,只要 a \neq b

    但是你可能会问当a=b时会发生什么?

    我们先假设它们处于相同的相位,完全相同的相位。好吧,它就是Sin²。Sin²是正的,在无穷远处求和,得到什么?你得到无穷大。所以,如果你有两个频率相同的正弦波,你会得到一个无穷大的值。除非它们完全不相合,在这种情况下是0。换句话说,Sin乘以Cos继续下去。但最基本的思路是,如果我整合一个函数,如果它是由一个与\omega不同的正弦曲线组成的,当我采用\omega的正弦并且我完全间隔时,我什么也得不到。但是如果它是等于那个正弦,我会得到无穷大。而且这将是这样写的。

    所以让我们做一个简单的例子。假设有一个简单的函数:f(x)=cos(2 \pi \omega x)

    然后我们选一些频率,我们把它叫做u,简单点。

    如果我取这个积分,Okay。那么如果 u 等于那个相同的 \omega 那将是无穷的,否则它将是零。

    所以看起来像这样, Okay。我们只有这两个成员脉冲。这里只有这两个无穷大的尖峰。

    这被称为与余弦对应的脉冲。你可以看到它们是正的,这是因为它是余弦。

    如果我们有一个正弦,因为sin (- x) = - sin (x),它看起来就像在这里的向上,

    和另一个将会在这里向下。

    这就是正弦,Okay。这就是虚部部分。两个公式如下:

    我们只是在计算一个基集来说明这个正弦曲线有多少。我们可以为所有频率做到这一点。但是你可能会问,我们不是也必须为所有相位做到这一点吗?

    答案是否定的。最酷的事情之一是,如果我有cos(\omega x) 和 sin(\omega x)。所以1是90度 或 两相相变的功率。我可以通过它们的线性组合得到任意相位。你可以证明,这是真的,我们可以做一个小演示。但基本上,如果我有一个任意相位的正弦曲线,如果你告诉我它对余弦的积分是多少、它对正弦的积分是多少,我就能告诉你这个正弦曲线的相位是多少。


    复数的基础:

    https://blog.csdn.net/sw3300255/article/details/83149483


    ——学会编写自己的代码,才能练出真功夫。

    展开全文
  • 离散傅里叶变换与相关性计算

    千次阅读 2019-11-03 11:41:37
    文章目录离散傅里叶变换(DFT)与相关性计算基本定义时域的卷积对应频域的相乘循环平移的傅里叶变换相关与卷积通过傅里叶变换计算相关 本文主要介绍通过离散傅里叶变换的方法计算两个时间序列的相关性,不仅给出...

    离散傅里叶变换(DFT)与相关性计算


    本文主要介绍通过离散傅里叶变换的方法计算两个时间序列的相关性,不仅给出数学公式的推导,同时给出对应的 Matlab程序的验证。文中涉及到相关,卷积等基本概念,本文假设读者已经具备相关的数学基础。

    更改说明:由于CSDN的MD语法在解析部分公式时出现错误,所以讲部分公式改为图片。
    补充:关于循环卷积与频域的相乘关系

    基本定义

    首先,我们给出离散傅里叶变换,卷积与相关的定义:

    离散傅里叶变换

    长度为N的离散时间序列 x ( n ) x(n) x(n)的离散傅里叶变换的定义为:
    X ( k ) = ∑ n = 0 N − 1 x ( n ) e − j 2 π n k / N X(k) = \sum_{n=0}^{N-1}x(n)e^{-j2\pi nk/N} X(k)=n=0N1x(n)ej2πnk/N
    对应的Matlab函数为fft(x)

    卷积

    长度为N的离散时间序列 x ( n ) x(n) x(n) y ( n ) y(n) y(n)的卷积定义为:
    z ( n ) = ∑ m = 0 N − 1 x ( m ) y ( n − m ) z(n) = \sum^{N-1}_{m=0} x(m)y(n-m) z(n)=m=0N1x(m)y(nm)
    得到的卷积 z ( n ) z(n) z(n)的长度为 2 N − 1 2N-1 2N1,n的取值范围是0到2N-2,卷积对应的Matlab函数为conv(x, y)。举个简单的例子:

    x = [1, 2, 3]
    y = [1, -1]
    z = conv(x, y) %  [1, 1, 1, -3]
    

    相关

    长度为N的离散时间序列 x ( n ) x(n) x(n)的相关定义为:
    z ( n ) = ∑ m = 0 N − 1 x ( m ) y ( m − n ) z(n) = \sum^{N-1}_{m=0} x(m)y(m-n) z(n)=m=0N1x(m)y(mn)
    同样,得到的相关 z ( n ) z(n) z(n)的长度为 2 N − 1 2N-1 2N1,卷积对应的Matlab函数为xcorr(x, y), 注意n的范围是 − ( N − 1 ) -(N-1) (N1) ( N − 1 ) (N-1) (N1)。举个简单的例子:

    x = [1, 2, 3]
    y = [1, -1, 1]
    z = xcorr(x, y) %  [1, 1, 2, -1, 3],对应的n是[-2, -1, 0, 1, 2]
    

    注意xcorr中要求两个时间序列长度相等,如果不相等,会自动将较短的时间序列补0。

    时域的卷积对应频域的相乘

    z ( n ) z(n) z(n)是时间序列 x ( n ) x(n) x(n) y ( n ) y(n) y(n)的卷积,长度为 2 N − 1 2N-1 2N1,则对应的离散傅里叶变换为:

    Z ( k ) = ∑ n = 0 2 N − 2 z ( n ) e − j 2 π n k / 2 N − 1 = ∑ n = 0 2 N − 2 ∑ m = 0 N − 1 x ( m ) y ( n − m ) e − j 2 π n k / 2 N − 1 \begin{aligned} Z(k) & = \sum^{2N-2}_{n=0} z(n)e^{-j2\pi nk/{2N-1}} \\ & = \sum^{2N-2}_{n=0} \sum^{N-1}_{m=0} x(m)y(n-m)e^{-j2\pi nk/{2N-1}} \end{aligned} Z(k)=n=02N2z(n)ej2πnk/2N1=n=02N2m=0N1x(m)y(nm)ej2πnk/2N1
    交换求和顺序得
    Z ( k ) = ∑ m = 0 N − 1 x ( m ) ∑ n = 0 2 N − 2 y ( n − m ) e − j 2 π n k / 2 N − 1 = ∑ m = 0 N − 1 x ( m ) e − j 2 π m k / 2 N − 1 ∑ n = 0 2 N − 2 y ( n − m ) e − j 2 π ( n − m ) k / 2 N − 1 = ∑ m = 0 N − 1 x ( m ) e − j 2 π m k / 2 N − 1 ∑ n = 0 N − 1 y ( n ) e − j 2 π n k / 2 N − 1 \begin{aligned} Z(k) & = \sum^{N-1}_{m=0} x(m) \sum^{2N-2}_{n=0} y(n-m)e^{-j2\pi nk/{2N-1}} \\ & = \sum^{N-1}_{m=0} x(m)e^{-j2\pi mk/{2N-1}} \sum^{2N-2}_{n=0} y(n-m)e^{-j2\pi (n-m)k/{2N-1}} \\ & = \sum^{N-1}_{m=0} x(m)e^{-j2\pi mk/{2N-1}} \sum^{N-1}_{n=0} y(n)e^{-j2\pi nk/{2N-1}} \end{aligned} Z(k)=m=0N1x(m)n=02N2y(nm)ej2πnk/2N1=m=0N1x(m)ej2πmk/2N1n=02N2y(nm)ej2π(nm)k/2N1=m=0N1x(m)ej2πmk/2N1n=0N1y(n)ej2πnk/2N1
    注意等式(7)中第二项的变量替换。由于 Z ( k ) Z(k) Z(k)序列的长度为 2 N − 1 2N-1 2N1,而 X ( k ) X(k) X(k) Y ( k ) Y(k) Y(k)的长度都为 N N N,为了长度上一致,我们将 x ( n ) x(n) x(n) y ( n ) y(n) y(n)分别补0,使得最终的序列长度为 2 N − 1 2N -1 2N1,分别记作 x p ( n ) x_p(n) xp(n) y p ( n ) y_p(n) yp(n),以 x p ( n ) x_p(n) xp(n)为例对应的数学描述为:

    x p ( n ) = { x ( n ) n ≤ N − 1 0 n > N − 1 x_p(n) = \left\{ \begin{array}{rl} x(n) & n \leq N-1\\ 0 & n>N-1 \end{array} \right. xp(n)={x(n)0nN1n>N1
    则其对应的傅里叶变换为:
    X p ( k ) = ∑ n = 0 N − 1 x p ( n ) e − j 2 π n k / 2 N − 1 X_p(k) = \sum_{n=0}^{N-1}x_p(n)e^{-j2\pi nk/2N-1} Xp(k)=n=0N1xp(n)ej2πnk/2N1
    联系等式(8)可以得到
    Z ( k ) = X p ( k ) Y p ( k ) Z(k) = X_p(k)Y_p(k) Z(k)=Xp(k)Yp(k)
    我们使用一个简单的Matlab程序进行验证。

    N = 20;
    x = randn(1, N);
    y = randn(1, N);
    
    z = conv(x, y);
    
    xp = [x, zeros(1, N-1)]; %补0
    yp = [y, zeros(1, N-1)]; %补0
    
    zk = fft(z);
    
    xpk = fft(xp);
    ypk = fft(yp);
    zke = xpk.*ypk;%频域相乘
    
    real_err = max(abs(real(zk)-real(zke)));
    imag_err = max(abs(imag(zk) - imag(zke)));
    
    fprintf('实补误差 %e\n虚部误差 %e\n',real_err, imag_err);
    
    %% 程序输出结果
    % 实部误差 1.421085e-14
    % 虚部误差 8.881784e-15
    

    其他的文献补0的个数都选择N而非N-1,背后的道理是一样的。选择N的好处是 X p ( k ) X_p(k) Xp(k) X ( k ) X(k) X(k)有一个一致的对应关系,即 X p ( 2 k ) = X ( k ) X_p(2k) = X(k) Xp(2k)=X(k)

    上面的定理可以简单描述为时域的卷积对应频域的乘积,一个典型应用加速卷积的计算。我们知道,如果直接套用卷积的定义,那么时间复杂度是 O ( n 2 ) O(n^2) O(n2),而快速傅里叶变换的时间复杂度是 O ( n log ⁡ n ) O(n\log n) O(nlogn),一个直接的思路就是利用傅里叶变换计算卷积,即首先根据 X p ( k ) X_p(k) Xp(k) Y p ( k ) Y_p(k) Yp(k)得到时间 Z ( k ) Z(k) Z(k),然后进行傅里叶反变换得到卷积 z ( n ) z(n) z(n)

    循环卷积与相乘的关系

    在图像处理中,图像与卷积核的运算通常需要保持图像大小不变,这种情况下时域的卷积与频域的乘积的对应关系是否仍然保持呢?
    Matlab对应的卷积函数为imfilter(u, v,'circular', 'conv'),给出简单的Matlab函数进行验证。

    a=[1,2,3,4,5,6,7];
    b = [1,2,3];
    d = imfilter(a,b,'circular','conv');
    % 结果 d => [ 11     8    12    16    20    24    21]
    

    了解了循环卷积的运算过程,我们进行公式的推导。首先,对循环卷积进行表示。对于长度为N的序列 x ( n ) x(n) x(n)其中,n=0, 1, 2, …N-1。以及长度为M的序列 h ( n ) h(n) h(n),这里对于 h ( n ) h(n) h(n),如果长度为奇数2M+1,则对应的n的取值范围为-M, -(M-1), … 0, …, (M-1), M。如果长度为偶数2M,则对应的n的范围表示为-M, …M-1。则循环卷积y(n)可以表示为:
    y ( n ) = ∑ m = 0 N − 1 x ( m ) h ) ( n − m ) , n = 0 , 1 , 2 , . . . N − 1 y(n) = \sum_{m=0}^{N-1}x(m)h)(n-m), n=0,1,2,...N-1 y(n)=m=0N1x(m)h)(nm),n=0,1,2,...N1
    根据上述公式,进行形同推导,则可以得到:
    Y ( k ) = ∑ m = 0 N − 1 x ( m ) ∑ n = 0 N − 1 h ( n − m ) e − j 2 π n k / N = ∑ m = 0 N − 1 x ( m ) e − j 2 π m k / N ∑ n = 0 N − 1 h ( n − m ) e − j 2 π ( n − m ) k / N \begin{aligned} Y(k) & = \sum^{N-1}_{m=0} x(m) \sum^{N-1}_{n=0} h(n-m)e^{-j2\pi nk/N} \\ & = \sum^{N-1}_{m=0} x(m)e^{-j2\pi mk/N} \sum^{N-1}_{n=0} h(n-m)e^{-j2\pi (n-m)k/N} \\ \end{aligned} Y(k)=m=0N1x(m)n=0N1h(nm)ej2πnk/N=m=0N1x(m)ej2πmk/Nn=0N1h(nm)ej2π(nm)k/N
    定义 h p ( n ) h_p(n) hp(n)为将 h ( n ) h(n) h(n)补0到长度为N,然后将向左循环平移M位。即:
    h p ( n ) = { h ( n ) 0 ≤ n ≤ M − 1 x ( n − m + N ) n ≥ n − M 0 o t h e r w i s e h_p(n) = \left\{ \begin{array}{rl} h(n) & 0\leq n \leq M-1\\ x(n-m+N) & n \geq n-M\\ 0 & otherwise \end{array} \right. hp(n)=h(n)x(nm+N)00nM1nnMotherwise
    Y ( k ) = X ( k ) H p ( k ) Y(k)=X(k)H_p(k) Y(k)=X(k)Hp(k),我们通过一个Matlab程序进行验证。

    clc
    clear
    x=[1,2,3,4,5,6,7];
    h = [1,2,-1];
    
    y = imfilter(x, h, 'circular', 'conv'); %-3     6     8    10    12    14     9
    
    
    hp = [h, zeros(1, length(x)-length(h))]; %补0
    m = floor(length(h)/2);
    hp = circshift(hp, -m); %向左循环平移
    Yk = fft(x).*fft(hp);
    
    ye = ifft(Yk);  %  -3.0000    6.0000    8.0000   10.0000   12.0000   14.0000    9.0000
    % 计算结果可知,y与ye的结果一致(忽略舍入误差)
    

    循环平移的傅里叶变换

    循环平移的本质是将序列 x ( n ) x(n) x(n)看做周期为N的周期序列,也就是N到2N-1的取值与0到N-1的取值有一一对应的关系。这里我们给出证明:

    将序列 x ( n ) x(n) x(n)向右做循环平移m位,得到y(n),则:
    y ( n ) = { x ( n − m ) n ≥ m x ( n − m + N ) n < m y(n) = \left\{ \begin{array}{rl} x(n-m) & n \geq m\\ x(n-m+N) & n < m \end{array} \right. y(n)={x(nm)x(nm+N)nmn<m
    则对应的傅里叶变换为:
    Y ( k ) = ∑ n = 0 N − 1 y ( n ) e − j 2 π n k / N = ∑ n = m N − 1 x ( n − m ) e − j 2 π n k / N + ∑ n = 0 m − 1 x ( n − m + N ) e − j 2 π n k / N = ∑ n = m N − 1 x ( n − m ) e − j 2 π ( n − m ) k / N e − j 2 π m k / N + ∑ n = 0 m − 1 x ( n − m + N ) e − j 2 π ( n − m + N ) k / N e − j 2 π ( m − N ) k / N = ∑ n = 0 N − 1 x ( n ) e − j 2 π n k / N e − j 2 π m k / N = X ( k ) e − j 2 π m k / N \begin{aligned} Y(k) &= \sum_{n=0}^{N-1} y(n)e^{-j2\pi nk/N} \\ &= \sum_{n=m}^{N-1} x(n-m)e^{-j2\pi nk/N} + \sum_{n=0}^{m-1} x(n-m+N)e^{-j2\pi nk/N} \\ &= \sum_{n=m}^{N-1} x(n-m)e^{-j2\pi (n-m)k/N}e^{-j2\pi mk/N} + \sum_{n=0}^{m-1} x(n-m+N)e^{-j2\pi (n-m+N)k/N}e^{-j2\pi (m-N)k/N} \\ &= \sum_{n=0}^{N-1}x(n)e^{-j2\pi nk/N}e^{-j2\pi mk/N} \\ &= X(k)e^{-j2\pi mk/N} \end{aligned} Y(k)=n=0N1y(n)ej2πnk/N=n=mN1x(nm)ej2πnk/N+n=0m1x(nm+N)ej2πnk/N=n=mN1x(nm)ej2π(nm)k/Nej2πmk/N+n=0m1x(nm+N)ej2π(nm+N)k/Nej2π(mN)k/N=n=0N1x(n)ej2πnk/Nej2πmk/N=X(k)ej2πmk/N

    相关与卷积

    熟悉卷积的读者应该知道卷积的计算过程分为两步:翻转和平移,与相关的计算相比,只是多了一步翻转操作。如果直接序列进行翻转,是否能用卷积代替相关操作呢?答案是是显然的,我们用一个Matlab程序做证明。

    x = [1,2,1];
    y = [1,2,3];
    
    z1 = xcorr(x, y); %相关
    z2 = conv(x, flip(y)); %先翻转,然后卷积
    
    z1 % => [3, 8, 8, 4, 1]
    z2 % => [3, 8, 8, 4, 1]
    

    这里不做数学的证明,有兴趣的读者可以自己证明。证明的过程中注意相关计算结果的第一项对应的n是 − ( N − 1 ) -(N-1) (N1),而卷积操作的第一项对应的n是0。

    通过傅里叶变换计算相关

    如果直接利用相关的公式,时间复杂度是 O ( n 2 ) O(n^2) O(n2),而快速傅里叶变换的时间复杂度是 O ( n log ⁡ n ) O(n\log n) O(nlogn),那么能否使用fft方法来计算相关,从而提高运算速度呢?上一节我们已经知道相关与翻转后做卷积是等价的,那么一个直接的思路是:

    1. 将序列 y ( n ) y(n) y(n)翻转,得到 u ( n ) u(n) u(n),然后补0,得到 u p ( n ) u_p(n) up(n),进行傅里叶变换得到 U p ( k ) U_p(k) Up(k)
    2. 将序列 x ( n ) x(n) x(n)补0,得到 x p ( n ) x_p(n) xp(n),进行傅里叶变换得到 X p ( k ) X_p(k) Xp(k)
    3. 计算 Z ( k ) = X ( p ) U p ( k ) Z(k) = X_(p)U_p(k) Z(k)=X(p)Up(k),通过傅里叶反变换得到 z ( n ) z(n) z(n)

    上述计算的缺点是多了一个中间变量 u ( n ) u(n) u(n),那么能否用 Y p ( k ) Y_p(k) Yp(k) 来计算 U p ( k ) U_p(k) Up(k)呢?我们进行一个简单的推导。
    U p ( k ) = ∑ n = 0 N − 1 u ( n ) e − j 2 π n k / 2 N − 1 = ∑ n = 0 N − 1 y ( N − 1 − n ) e − j 2 π n k / 2 N − 1 = ∑ n = 0 N − 1 y ( N − 1 − n ) e − j 2 π ( n + 1 − N ) k / 2 N − 1 e − j 2 π ( N − 1 ) k / 2 N − 1 = { ∑ n = 0 N − 1 y ( N − 1 − n ) e − j 2 π ( N − 1 − n ) k / 2 N − 1 } ∗ e − j 2 π ( N − 1 ) k / 2 N − 1 = Y p ∗ ( k ) e − j 2 π ( N − 1 ) k / 2 N − 1 \begin{aligned} U_p(k) & = \sum_{n=0}^{N-1}u(n)e^{-j2\pi nk/2N-1} \\ & = \sum_{n=0}^{N-1}y(N-1-n)e^{-j2\pi nk/2N-1} \\ & = \sum_{n=0}^{N-1}y(N-1-n)e^{-j2\pi (n+1-N)k/2N-1}e^{-j2\pi (N-1)k/2N-1} \\ & = \left\{\sum_{n=0}^{N-1}y(N-1-n)e^{-j2\pi (N-1-n)k/2N-1}\right\}^*e^{-j2\pi (N-1)k/2N-1} \\ & = Y^*_p(k)e^{-j2\pi (N-1)k/2N-1} \end{aligned} Up(k)=n=0N1u(n)ej2πnk/2N1=n=0N1y(N1n)ej2πnk/2N1=n=0N1y(N1n)ej2π(n+1N)k/2N1ej2π(N1)k/2N1={n=0N1y(N1n)ej2π(N1n)k/2N1}ej2π(N1)k/2N1=Yp(k)ej2π(N1)k/2N1
    其中 ∗ * 表示共轭运算,对应的matlab函数为conj(X)

    也就是说如果令 Z 1 ( k ) = X p ( k ) Y p ∗ ( k ) Z_1(k) = X_p(k)Y_p^*(k) Z1(k)=Xp(k)Yp(k),对应的傅里叶反变换为 z 1 ( n ) z_1(n) z1(n),则 Z ( k ) = Z 1 ( k ) e − j 2 π ( N − 1 ) k / 2 N − 1 Z(k) = Z_1(k)e^{-j2\pi (N-1)k/2N-1} Z(k)=Z1(k)ej2π(N1)k/2N1,根据循环平移的傅里叶变换易知 z ( n ) = z 1 ( n − ( N − 1 ) ) z(n) = z1(n-(N-1)) z(n)=z1(n(N1))

    因此,通过傅里叶变换计算相关运算的思路为:

    1. x ( n ) x(n) x(n) y ( n ) y(n) y(n)分别补N-1个0,得到 x p ( n ) x_p(n) xp(n) y p ( n ) y_p(n) yp(n),进行快速傅里叶变换得到 X p ( k ) X_p(k) Xp(k) Y p ( k ) Y_p(k) Yp(k)
    2. 计算 Z 1 ( k ) = X ( k ) Y ∗ ( k ) Z_1(k) = X(k)Y^*(k) Z1(k)=X(k)Y(k),并进行傅里叶反变换得到 z 1 ( n ) z_1(n) z1(n)
    3. 进行循环平移得到 z ( n ) z(n) z(n)

    同之前讨论过的一样,补0的个数只要大于等于N-1不影响计算结果,很多资料中都补N个0.

    最后我们通过一个Matlab程序对比一下4种计算相关性的方法。

    % 通过4种方式计算相关
    %% 数据生成
    N = 30;
    x = randn(1, N);
    y = [randn(1, 4), x(1:N-4)];
    %y = randn(1, N);
    
    %% 方法1, 通过Matlab自带函数xcorr进行计算
    
    z1 = xcorr(x, y);
    
    %% 方法2,通过翻转,求卷积
    
    z2 = conv(x, flip(y));
    
    %% 方法3,通过翻转,傅里叶变换计算
    u = flip(y);
    xp = [x, zeros(1, N-1)];
    up = [u, zeros(1, N-1)];
    
    zk = fft(xp).*fft(up);
    z3 = ifft(zk);
    
    %% 方法4,傅里叶变换方法
    xp = [x, zeros(1, N-1)];
    
    yp = [y, zeros(1, N-1)];
    
    z1k = fft(xp).*conj(fft(yp));
    z1n =  ifft(z1k);
    z4 = circshift(z1n, N-1); %需要平移操作
    
    %% 计算结果的可视化
    lags = -(N-1):N-1;
    
    
    scalediff = 5;
    plot(lags, z1);
    hold on;
    plot(lags, z2 + scalediff*1);
    plot(lags, z3 + scalediff*2);
    plot(lags, z4 + scalediff*3);
    
    plot([-4, -4], ylim, 'r:','LineWidth', 1.5); % 平移量
    hold off
    legend('方法一','方法二','方法三', '方法四');
    
    

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

    展开全文
  • 傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换傅里叶变换
  • N点离散傅里叶变换同时计算两个N点实序列的离散傅里叶变换
  • 采用坐标扩展变换,突破了快速傅里叶变换计算过程中输入屏和衍射屏空间尺度必须相同、抽样格点必须等间隔的限制,使上述问题得到解决,可以得到更加详细的聚焦光束光场分布。同时,采用分两步计算的思想,避免了计算...
  • 十分简明易懂的FFT(快速傅里叶变换

    万次阅读 多人点赞 2018-08-07 11:38:56
    快速傅里叶变换 (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 ( n 2 ) O(n^2) O(n2),但 F F T FFT FFT O ( n log ⁡ 2 n ) O(n\log_2 n) O(nlog2n)的时间解决
    F F T FFT FFT名字逼格高,也难懂,其他教程写得让人看不太懂,于是自己随便写一下

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

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


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

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

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

    系数表示法

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

    点值表示法

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

    • n n n个不同的 x x x代入,会得出 n n n个不同的 y y y,在坐标系内就是 n n n个不同的点

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

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

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


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

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

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

    什么意思呢?

    设两个点值多项式分别为 f ( x ) = { ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , ( x 2 , f ( x 2 ) ) , . . . , ( x n − 1 , f ( x n − 1 ) ) } f(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2)),...,(x_{n-1},f(x_{n-1}))\} f(x)={(x0,f(x0)),(x1,f(x1)),(x2,f(x2)),...,(xn1,f(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 ) ) } g(x)=\{(x_0,g(x_0)),(x_1,g(x_1)),(x_2,g(x_2)),...,(x_{n-1},g(x_{n-1}))\} g(x)={(x0,g(x0)),(x1,g(x1)),(x2,g(x2)),...,(xn1,g(xn1))}
    设它们的乘积是 h ( x ) h(x) h(x),那么 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 ) ) } 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}))\} h(x)={(x0,f(x0)g(x0)),(x1,f(x1)g(x1)),...,(xn1,f(xn1)g(xn1))}

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

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

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

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

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


    DFT前置知识&技能

    复数

    毕竟高中有所以不多说

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

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

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

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


    x x x轴就是实数集中的坐标轴, y y y轴就是虚数单位 i i i

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

    复数的运算

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

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


    这张是从网上盗的

    A B + A D = A C AB+AD=AC AB+AD=AC

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


    DFT(离散傅里叶变换)

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

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

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

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

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

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

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


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

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


    单位根的一些性质

    F F T FFT FFT的过程中需要用到 ω \omega ω的一些性质

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

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

    • 证明

    • ω n k = c o s k n 2 π + i s i n k n 2 π = c o s 2 k 2 n 2 π + i s i n 2 k 2 n 2 π = ω 2 n 2 k \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=cosnk2π+isinnk2π=cos2n2k2π+isin2n2k2π=ω2n2k

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

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

    • 证明

    • ω n n 2 = c o s n 2 n 2 π + i s i n n 2 n 2 π = c o s π + i s i n π = − 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 ωn2n=cosn2n2π+isinn2n2π=cosπ+isinπ=1

    • (这个东西和 e i x = c o s x + i s i n x e^{ix}=cosx+isinx eix=cosx+isinx e i π + 1 = 0 e^{i\pi}+1=0 eiπ+1=0有点关系,我不会就不讲了

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

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

    FFT(快速傅里叶变换)

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

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

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

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

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

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

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

    A 1 ( x ) = a 0 + a 2 x + a 4 x 2 + . . . + a n − 2 x n 2 − 1 A_1(x)=a_0+a_2x+a_4x^2+...+a_{n-2}x^{{n\over 2}-1} A1(x)=a0+a2x+a4x2+...+an2x2n1 A 2 ( x ) = a 1 + a 3 x + a 5 x 2 + . . . + a n − 1 x n 2 − 1 A_2(x)=a_1+a_3x+a_5x^2+...+a_{n-1}x^{{n \over 2}-1} A2(x)=a1+a3x+a5x2+...+an1x2n1

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

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

    A ( ω n k ) = A 1 ( ( ω n k ) 2 ) + ω n k A 2 ( ( ω n k ) 2 ) A(\omega^k_n)=A_1((\omega^k_n)^2)+\omega^k_nA_2((\omega^k_n)^2) A(ωnk)=A1((ωnk)2)+ωnkA2((ωnk)2) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) = A 1 ( ω n 2 k ) + ω n k A 2 ( ω n 2 k ) =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}) =A1(ωn2k)+ωnkA2(ωn2k)=A1(ω2nk)+ωnkA2(ω2nk)

    那么对于 A ( ω n k + n 2 ) A(\omega^{k+{n\over2}}_n) A(ωnk+2n)的话,代入 ω n k + n 2 \omega^{k+{n \over 2}}_n ωnk+2n
    A ( ω n k + n 2 ) = A 1 ( ω n 2 k + n ) + ω n k + n 2 A 2 ( ω n 2 k + 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) A(ωnk+2n)=A1(ωn2k+n)+ωnk+2nA2(ωn2k+n) = A 1 ( ω n 2 k ω n n ) − ω n k A 2 ( ω n 2 k ω n n ) =A_1(\omega^{2k}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^n_n) =A1(ωn2kωnn)ωnkA2(ωn2kωnn) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k ) = A 1 ( ω n 2 k ) − ω n k A 2 ( ω n 2 k ) =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}) =A1(ωn2k)ωnkA2(ωn2k)=A1(ω2nk)ωnkA2(ω2nk)

    • 发现了什么?

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

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


    IFFT(快速傅里叶逆变换)

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

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

    • I D F T IDFT IDFT暴力 O ( n 2 ) O(n^2) O(n2)做?其实也可以用 F F T FFT FFT O ( n log ⁡ 2 n ) O(n\log_2n) O(nlog2n)的时间搞

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

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

    意思就是说 F F T FFT FFT I F F T IFFT IFFT可以一起搞


    朴素版FFT板子

    c + + c++ c++有自带的复数模板 c o m p l e x complex complex
    a . r e a l ( ) a.real() a.real()即表示复数 a a a的实部

    #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 , b a,b a,b相乘再转系数多项式 c c c,通常只用打这么一小段

    	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只能处理 n n n 2 2 2的整数次幂的多项式
    所以在 F F T FFT FFT前一定要把 n n n调到 2 2 2的次幂

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

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


    FFTの优化——迭代版FFT

    这个图也是盗的

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

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

    这样的话我们可以先把原序列变换好,把每个数放在最终的位置上,再一步一步向上合并
    一句话就可以 O ( n ) O(n) O(n)预处理第 i i i位最终的位置 r e v [ i ] 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后记

    本人版权意识薄弱……

    N T T NTT NTT我来了

    展开全文
  • 用python实现傅里叶的快速变换,实现与传统傅里叶变换不同的快速傅里叶变换,也使用了一些图形界面设计
  • 离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换离散傅里叶变换
  • 傅里叶变换.pdf傅里叶变换.pdf傅里叶变换.pdf傅里叶变换.pdf
  • 实现图像的傅里叶变换通过OPENCV,应用经典的傅里叶变换,实现图像变换
  • 傅里叶变换(二维离散傅里叶变换)

    万次阅读 多人点赞 2018-06-15 22:22:35
    离散二维傅里叶变换一常用性质: 可分离性、周期性和共轭对称性、平移性、旋转...根据快速傅里叶变换计算要求,需要图像的行数、列数均满足2的n次方,如果不满足,在计算FFT之前先要对图像补零以满足2的n次。 ...
  • 傅里叶变换

    2018-06-13 20:40:41
    傅里叶变换深入了解傅里叶变换深入了解傅里叶变换深入了解傅里叶变换深入了解
  • 条纹投影轮廓测量中不连续点的相位计算傅里叶变换,开窗傅里叶变换和小波变换方法的比较
  • 离散傅里叶级数,离散傅里叶变换及逆傅里叶变换的实现。
  • 精通matlab傅里叶变换及逆变换和快速傅里叶变换
  • 傅里叶变换到加窗傅里叶变换

    万次阅读 多人点赞 2017-12-10 14:00:49
    傅里叶变换到加窗傅里叶变换我们都做了什么
  • 傅里叶变换 一维离散傅里叶变换

    万次阅读 热门讨论 2019-11-06 21:08:43
    DFT:(Discrete Fourier Transform)离散傅里叶变换傅里叶变换在时域和频域上都呈离散的形式,将信号的时域采样变换为其DTFT的频域采样。在形式上,变换两端(时域和频域上)的序列是有限长的,而实际上这两组序列...
  • 傅里叶变换&短时傅里叶变换&小波变换

    万次阅读 多人点赞 2016-05-03 11:26:36
    傅里叶变换&短时傅里叶变换&小波变换
  • 计算信号的傅里叶变换,信号的反傅里叶变换,周期信号的展开,
  • 离散傅里叶变换和快速傅里叶变换,希望对大家有帮助
  • C语言快速红丝线傅里叶变换的方法,仅供参考
  • OpenCV中的图像变换——傅里叶变换

    万次阅读 多人点赞 2021-07-22 20:11:24
    这篇博客将介绍OpenCV中的图像变换,包括用Numpy、OpenCV计算图像的傅里叶变换,以及傅里叶变换的一些应用;
  • 连续傅里叶变换 傅里叶变换的性质 离散傅里叶变换(DFT) 从前面我们已经知道,非周期连续函数傅里叶变换如下 F(ω)=∫−∞+∞f(t)e−iωtdt F(\omega)=\int ^{+\infty}_{-\infty}f(t)e^{-i\omega t}dt F(ω)=∫−...
  • 傅里叶变换(一)——认识傅里叶变换

    万次阅读 多人点赞 2018-05-29 22:36:07
    注:本文为博主参考书籍和他人文章并加上自己的理解所编,作为学习笔记使用并将其分享出去供大家学习。若涉及到引用您的文章内容请评论区告知!...一、什么是傅里叶变换   时域及频域  在讲...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 43,398
精华内容 17,359
关键字:

傅里叶变换如何计算