离散高斯核是DFT的本征函数吗？

#1 楼

由于DFT可通过与傅立叶矩阵相乘来表示，因此您的问题等同于询问傅立叶矩阵的特征向量是什么。 org / wiki / Discrete_Fourier_transform＃Eigenvalues_and_eigenvectors）。然而，由于特征值（$ 1，-1，i，-i $）并不简单，因此特征向量也不是唯一的（即线性组合也是特征向量） ）。也没有简单的封闭公式。

维基百科提供了一个与本征矢量接近的特征向量公式。

$$ F_m = \ sum_ {k =-\ infty} ^ \ infty \ exp \ left（-\ frac {\ pi \ cdot（m + N \ cdot k）^ 2} {N} \ right）\ quad \ quad \ quad m = 0，\ ldots，N-1 $ $

因此得出结论，高斯函数本身不是特征向量，而是无限的高斯和。
当我们从FT转换到DFT时，无限和可能可以解释为等效于频域和时域的离散化。因此，它不像截断离散高斯那样简单。

#2 楼

在过去的几周中，我在这个问题上取得了巨大的进步。

窗口函数的清零正弦函数族

我发现了以前无法识别的窗口函数。这些特别令人惊奇的是，DFT的特征向量似乎是该族中的规则成员。这不仅为计算精确的最小特征向量提供了数值解决方案，而且为族提供了简洁的定义，可以证明它们是特征向量，还可以为更高阶的离散厄米-高斯函数构造新的生成器。
我仍在研究所有大家庭如何融合在一起的细节，并从初始条件中寻找差异方程式来定义特征向量。 。
这是清零正弦窗口的定义：
$$ W_L [n] = \ prod_ {l = 1} ^ {L} \ sin \ left（\ frac {n + l} { N} \ pi \ right）$$
这是凸起正弦窗的定义：
$$ W_L [n] = \ sin ^ L \ left（\ frac {n} {N} \ pi \ right）$$
前者是“出厂”版本，后者是“权力”版本。它们开始相似。....

更新：已发现用于生成单个成员的递归函数：
\ begin {aligned}
W_L [0]＆ = 0.0001 \\
W_L [n]＆= \ frac {\ sin \ left（\ frac {n + L} {N} \ pi \ right）} {\ sin \ left（\ frac {n} { N} \ pi \ right）} W_L [n-1] \\
\ end {aligned}
它直接来自于定义。
常数并不重要不为零。结果将是成比例的。
$ N $为奇数，而$ L =（N-1）/ 2 $，这将是DFT的最小宽度特征向量。
偶数请参见我的文章大小写修复。

对于任何怀疑的人，这里有一个快速的示例程序。只有形式为$ 4m + 1 $的奇数N会具有一个特征向量，在峰值处有一个样本。 $ L =（N-1）/ 2 $，峰值出现在$ L / 2 $。其他N需要稍微按摩一下才能将它们置于中心特征向量形式。

import numpy as np

#========================================================================
def main():

#---- Set parameters

        m = 4

        N = m * 4 + 1
        L = ( N - 1 ) >> 1

#---- Build the Zeroing Sine Family member

        w = np.zeros( N )

        w[0] = 2**(-L)
        
        for n in range( 1 , N ):
          w[n] = w[n-1] * np.sin( (n+L) * np.pi / N ) / np.sin( n * np.pi / N )
          print( " %3d %12.8f" % ( n, w[n] ) )

#---- Shift the function to be zero centered

        r = np.roll( w, N - (L >> 1) )

#---- Take the DFT and compare

        sr = np.fft.fft( r ) / np.sqrt( N )
        
        for n in range( N ):
          if r[n] == 0 :
             print( " %3d %12.8f %12.8f %12.8f     Undefined" % \
                    ( n, r[n], sr[n].real, sr[n].imag ) )
          else:
             print( " %3d %12.8f %12.8f %12.8f  %12.8f" % \
                    ( n, r[n], sr[n].real, sr[n].imag, sr[n].real / r[n] ) )

#========================================================================
main()

这是结果。

#3 楼

此答案是@CedronDawg的答案的补充，后者介绍了本征矢量族。更具体地说，该答案提出了三种算法和一种混合算法，用于为给定长度$ N $生成具有最大连续零个数的离散傅里叶变换（DFT）的唯一（非归一化）非负特征向量。这是正常的DFT，没有时间或频率上的半采样偏移。
迭代算法，$ \ sim O \ big（\ exp（N）\ big）$ time
这里是迭代数值算法查找满足所述约束的特征向量。算法原理是在每个域中重复执行定义的约束，并对结果求平均。在八度中：

function retval = constrain(x, N)
  zerocount = N-floor((N+2)/4)*2-1;
  if (zerocount > 0)
    firstzero = floor((N+10)/4);
    x(firstzero:firstzero+zerocount-1) = 0;
  endif
  x = abs(x);
  retval = x / x(1);
endfunction

function retval = geteig_iterative(N, precision = 0, x = ones(1, N))
  x = constrain(x, N);
  good = x;
  X = fft(x)/sqrt(N);
  error = sumsq(x-X);
  iter = 0;
  while (max(abs(X-x)) > precision)
    lasterror = error;
    X = constrain(X, N);
    temp = X;
    x = (x + X)/2;
    X = fft(x)/sqrt(N);
    error = sumsq(x-X);
    iter++;
    if (error >= lasterror)
      break;
    endif
    good = temp;
  endwhile
  printf("Iterations: %d\n", iter);
  retval = good;
endfunction

N = 21; x = geteig_iterative(N); plot([0:N-1], x, "x", [0:N-1], real(fft(x)/sqrt(N)), "+");

可以给该算法一个起点向量。否则，它以1的向量开始。参数precision将所需的精度定义为计算出的特征向量与其DFT之间的最大绝对差值，例如precision = 0.000000000001或默认的precision = 0，以在迭代不再改善解时停止算法。确定零的位置和数量。他们的公式属于经验数学领域，但也受符号化求解$ N $较小值的问题的经验指导，并部分地通过测试随机向量初始化是否会改变结果或是否仍然可以增加firstzero来提供指导。
对于$ N = 21 $，在zerocount迭代中获得结果：
图1.对于$ N = 21 $，迭代算法（✕）及其DFT（+）的输出。
不幸的是，该算法需要大量与$ N $成指数比例的迭代，因此对于较大的$ N $无效，并且对于某些不是很大的$ N $值也无法收敛$ 18 $和$ 22 $：
图2. zerocount所需的迭代次数，作为$ N $的函数。
尽管该算法在$ N $大的缺点，但在$ N $小的结果还是有用的，显示出适合算法构造的清晰结构。通过查看展开的特征向量非零部分的频域和其他z平面零点（以倍频程为2605来查看结构最好）： ）展开特征向量的非负子序列。在这些Z平面图中，极点（✕）可以忽略。不太容易看到，但是当$ \ operatorname {mod}（N，4）\ in \ {0，2 \} $中时，在$ z = -1 $处有一个双零。不包括$ N <3 $，因为如何拆分展开的特征向量以获得非零子序列是不明确的。
分析显示零分布中的$ 4 $周期模式是$ N的函数$。当$ \ operatorname {mod}（N，4）\ in \ {0，2 \} $中时，可以看到$ z = -1 $的双零。 $ \ operatorname {mod}（N，4）\ in \ {0，1 \} $中的前两列显示单位圆上零的规则间距，表示频域中的实零，对应于DFT箱。对于$ \ operatorname {mod}（N，4）\ in \ {2，3 \} $，有两个必要但有点随意的实零，一个在单位圆的内侧，一个在单位圆的外侧，其值唯一地定义为a。给定$ N $，这是通过对随机向量初始化进行测试来确认的。对称性约束要求额外的零对是倒数对，从而将未知数减少为一个，这是一个很容易用数字求解的问题，如以下几节中的快速算法所述。
用$ \ operatorname {mod}（N，4）= 1 $看到的零模式的简单性很有吸引力。选择$ \ operatorname {mod}（N，4）$来选择DFT的时移和/或频移变体可能足以吸引所有$ N $，以确保相同的简单性。但是，包含一个零频率的bin可能很重要，并且由于问题是关于通常的DFT类型，因此我将坚持此答案。
快速卷积算法，$ O \ big（N \ log ^ 2（N）\ big）$ time
@CedronDawg的方法是将特征向量构造为实值向量的乘积，每个向量最多创建两个零。这里提出了一种使用类似原理的快速算法，但是它采用了分而治之的方法，而不是将这些向量一一相乘。在这种卷积算法中，通过时域卷积创建频域零。 precision = 0.000000000001的卷积由快速傅里叶变换（FFT）加速，总时间复杂度为$ O \ big（N \ log ^ 2（N）\ big）$。
该算法的主要工作是创建一长列频域零。此类工作分为递归地创建和组合一对半长列车，直到要创建的列车长度只有3个或更小。可以简单而直接地表示与这样的短列车相对应的时域序列。对于奇长列车，中间零点是单独创建的。无需单独创建半列车，而是通过复杂的正弦调制将它们克隆并转移到适当的位置。这样，递归不会创建调用树，而只会创建长度为$ O \ log（N）$的短调用链。在Octave中：

function retval = getzeros(N, M)
  if (M == 1)
    retval = [1, -1];
  elseif (M == 2)
    retval = [1, -2*cos(pi/N), 1];
  elseif (M == 3)
    retval = [1, -1 - 2*cos(2*pi/N), 1 + 2*cos(2*pi/N), -1];
  elseif (rem(M, 2) == 0)
    half = getzeros(N, M/2);
    modulator = exp(i*[((0:M/2)-M/4)*2*pi*M/4/N]);
    retval = real(conv(conj(modulator).*half, modulator.*half));
  else
    half = getzeros(N, (M-1)/2);
    modulator = exp(i*[((0:(M-1)/2)-(M-1)/4)*2*pi*(M+1)/4/N]);
    retval = conv(real(conv(conj(modulator).*half, modulator.*half)), getzeros(N, 1));
  endif
endfunction
#[h, w] = freqz(getzeros(10, 4), [], 500); plot(abs(h));

function retval = geteig_convolutional(N)
  M = N-floor((N+2)/4)*2-1;
  x = getzeros(N, M);
  if (rem(M, 2) == 1)
    x = conv(x, getzeros(N, 1));
  endif
  x .*= (-1).^[0:length(x)-1];
  mid = (length(x) + 1)/2;
  if (rem(N, 4) == 2 || rem(N, 4) == 3)
    s = sum(x)/sqrt(N);
    x = conv(x, [1, 2*(s - x(mid + 1))/(x(mid) - s), 1]);
  endif
  retval = circshift(horzcat(x/x(mid), zeros(1, N-length(x))), (rem(N,4) < 2)-mid);
endfunction

N = 512; x = geteig_convolutional(N); plot([0:N-1], x, "x", [0:N-1], real(fft(x)/sqrt(N)), "+");
plot([0:N-1], x - real(fft(x)/sqrt(N)), "x");

此图： > 图5.对于$ N = 512 $，快速卷积算法的输出与其DFT之间的差异。
对于大约$ N> 750 $的较大值，卷积算法开始出现数值误差，我想部分是因为前后FFT –逆FFT（IFFT）。下一节介绍的乘法算法将更加稳定。
代码包括创建$ \ operatorname {mod}（N，4）\ in \ {2，3 \} $中所需的倒数零为此，从等式$ x [-1] + a \求解了变量$ a $的卷积$ [1，a，1] $的值，该值创建了倒数零。 x [0] + x [1] = X [0] + a \，X [0] + X [0] $。该方程式要求卷积后，时域原点的值（由等式的左侧表示）等于频域原点的值（右侧），其中$ x $和$ X $代表新生特征向量及其DFT快速卷积算法。$ O \ big（N \ log（N）\ big）$ time
这里是快速算法的乘法版本。它在构造特征向量的同一域中创建零。中间结果以对数刻度存储，以适应数字的大动态范围。在八度中：

function retval = get_log2_zeros_multiplicative(N, M)
  k = [0:N-1];
  if (M == 1)
    x = log2(abs(sin(pi*k/N)));
  elseif (M == 2)
    x = log2(abs(cos(pi/N) - cos(2*pi*k/N + pi/N)));
  elseif (M == 3)
    x = log2(abs(sin(pi*k/N).*(2*cos(2*pi/N) + 1) - sin(3*pi*k/N)));
  elseif (rem(M, 4) == 0)
    half = get_log2_zeros_multiplicative(N, M/2);
    x = circshift(half, M/4) + circshift(half, -M/4);
  elseif (rem(M, 4) == 1)
    half = get_log2_zeros_multiplicative(N, (M-1)/2);
    x = circshift(half, (M-1)/4+1) + circshift(half, -(M-1)/4) + get_log2_zeros_multiplicative(N, 1);
  elseif (rem(M, 4) == 2)
    half = get_log2_zeros_multiplicative(N, M/2);
    x = circshift(half, (M-2)/4) + circshift(half, -((M-2)/4+1));
  else
    half = get_log2_zeros_multiplicative(N, (M-1)/2);
    x = circshift(half, (M-3)/4+1) + circshift(half, -((M-3)/4+1)) + get_log2_zeros_multiplicative(N, 1);
  endif
  retval = x - round(max(x));
endfunction

function retval = geteig_multiplicative(N)
  M = N-floor((N+2)/4)*2-1;
  x = get_log2_zeros_multiplicative(N, M);  
  if (rem(M, 2) == 1)
    x += get_log2_zeros_multiplicative(N, 1);
  endif
  x = ifftshift(2.^x);
  firstzero = floor((N+10)/4);
  x(firstzero:firstzero+M-1) = 0;
  if (rem(N, 4) == 2 || rem(N, 4) == 3)
    x0 = sum(x)/sqrt(N);
    x1 = sum(x.*cos(2*pi*[0:N-1]/N))/sqrt(N);
    s = x(1);
    x .*= 2*cos(2*pi*[0:N-1]/N) + 2*(s - x1)/(x0 - s);
  endif
  retval = x/x(1);
endfunction

N = 512; x = geteig_multiplicative(N); plot([0:N-1], x - real(fft(x)/sqrt(N)), "x");

在$ N = 512 $时，快速乘法算法的误差非常小：
图6.快速乘法算法的输出与其输出之间的差异对于$ N = 512 $，DFT。对于$ N = 65536 $，使用快速乘法算法计算的特征向量与其DFT之间的绝对差（未显示，类似于图6）将在zplane(nonzeros(fftshift(geteig_iterative(N)))');之下。
混合乘法-迭代算法，$ O \ big（\ log（N）N \ big）$ time
仍可以使用迭代算法完善乘法算法的输出。在八度中：使用迭代细化，对于$ N = 65536 $，可以将使用快速乘法算法计算的特征向量与其DFT之间的绝对差（未显示，类似于图7）从conv以下减少到1e-12中9以下
结论
即使对于非常大的$ N $，乘积算法1e-12和混合乘积迭代算法3e-16也非常快，是用于生成这些特征向量的推荐算法。由于Octave的FFT是在双精度浮点中完成的，因此混合算法的迭代优化不会出现任何错误，这在某些应用中可能会比用于计算迭代起点的乘法算法产生更多的错误。与混合算法的输出不同，在展开的意义上，乘法算法的输出保证是单峰的（图8）。算法（粗蓝线）及其最适合的二次趋势线（细黑线），N = 65536 $。高斯函数的对数是二次的。紧密匹配表明，当$ N \ to \ infty $时，这种类型的特征向量接近高斯函数。对于任何$ N $，第一个零将出现在这样的图中，时间为$ \ pm $ 23，在当前情况下为$ \ pm 16385 $，这意味着高斯和本征向量都将跳水到极小的值距离零已经很远了。