如何获得傅里叶系数以使用 DFT 绘制任何形状?

信息处理 自由度 傅里叶级数
2022-01-05 18:34:12

我正在自学傅立叶级数和 DFT,并试图通过youtube 上 Mathologer详细介绍的傅立叶本轮绘制一个程式化的$\pi$符号(从 18:39 开始),以及youtube 上 3Blue1Brown的非凡动画的出色解释.

目标是生成如下内容: 在此处输入图像描述

使用复数傅立叶级数:

$$z(t)=\sum\limits_{k=-\infty}^{\infty}{c_k \, e^{ikt}}$$

具有复系数:

$$c_k=\frac{1}{2\pi}\int_\limits{-\pi}^{\pi}z(t) \, e^{-ikt} \, \mathrm{d}t$$

我已经能够为$c_k=-2 < k < 2$生成“胚胎” $\pi$形状并获得与 Mathologer (@19:19) 相同的结果,但这仅仅是因为他提供了五个$c_k$值( @20:12)。这是我的输出: 傅立叶的胚胎\pi形状

回到目标:我为$\pi$符号创建了自己的 120 点坐标集:

\pi 坐标

我的问题是如何找到所有系数?我认为输入坐标需要是适合输入到 DFT 的等间距样本,但是尽管进行了很多搜索,但我仍然不确定从这里开始的过程是什么?

进度更新#3:

我度过了愉快的一天,在 MATLAB 中的各种算法上取得了很大进展。为了区分输出和输入$z$,我将$z_n$用于$N=120$复杂样本点$z(1),z(2), ... z(N)$$z_t $用于逆 DFT 后的$D=180$复数结果$z_t(1),z_t(2), ... z_t(D)$ 。这是我的$z_t$加上随机点$z_t(93)$的叠加图,显示了组件总和臂和相关的圆/本轮(注意 180 个点比上面绘制的原始 120 个点更靠近): 带有 zt(93) 示例的 pi 符号

下面显示了$z_t$ for $D=180$与$z_n$叠加以放大不准确性,并放大: pi DFT 细节 还有一段路要走;我真的很想以数学方式记录解决方案,并尝试各种方法来提高结果符号的准确性。但是我感觉我已经越过了山顶,现在这只是一个雪橇一路下来的例子!(著名遗言 :)

TIA 提供任何进一步的指导

PS:这是我的样本点坐标的链接(由于@Olli作为下面的答案上传,谢谢先生)。每行有一个$(x,y)$对,120 行:

链接到我的公共保管箱文件夹中的 ZIP 文件

这是 r bj 拼凑在一起绘制它的 MATLAB 程序(自 Chris 更新后),即使是第一个案例:

clearvars; close all
data = csvread("pi.csv"); % 121 rows with last repeating first
N = length(data) - 1;   % Chris added minus 1

inx = data(1:N,1);       % Chris was (:,1)
iny = data(1:N,2);       % ditto

Xk = fft(inx)/N;
Yk = fft(iny)/N;

X1 = Xk(1 : 1 + (N/2-1)     ); 
X4 = Xk(    1 + (N/2+1) : N );

% The main correction was here for X and Y: 
% the Nyquist freq must be allocated to one bin not two (previously)
Xnyq = Xk(1 + N/2);
X = [X1; Xnyq; X4];

Y1 = Yk(1 : 1 + (N/2-1)    );
Y4 = Yk(    1 + (N/2+1) : N);

Ynyq = Yk(1 + N/2); 
Y = [Y1; Ynyq; Y4];

x = N*ifft(X);
y = N*ifft(Y);

load('pi_zt_coords')
xt = real(ztt);
yt = imag(ztt);

plot(inx, iny,'o-','markersize',8)
hold on; grid on
plot(xt,yt,'k.-','markersize',8)
plot(x,  y,'mx')

xlim([100,250])
ylim([100,250])

legend('(x_{in} y_{in})','(x_t,y_t)','(x,y)','location','SouthEast')

title (['Even N =',num2str(N)]);

结果如下:

参数偶数情况

这里是一样的,但是去掉了一个点,所以这NN很奇怪。请注意,没有奈奎斯特值可以分成两部分(由 Chris 更新后) ODD 案例:

clearvars; close all

data = csvread("pi.csv");   % 121 rows with last repeating first
%data= csvread("pi_bandlimited.csv"); % from Olli's script - works too

data = vertcat(data(1:111,:), data(113:end,:)); % Delete row 112 to make N odd = 119

N = length(data) - 1;   % Chris added minus 1

inx = data(1:N,1);      % Chris (1:N,1) was (:,1)
iny = data(1:N,2);      % ditto

Xk = fft(inx)/N;
Yk = fft(iny)/N;

X1 = Xk(1 : 1 + (N-1)/2); 
X2 = Xk(1 + (N+1)/2 : N  );
X = [X1; X2];

Y1 = Yk(1 : 1 + (N-1)/2); 
Y2 = Yk(1 + (N+1)/2 : N);
Y = [Y1; Y2];

x = N*ifft(X);
y = N*ifft(Y);

load('pi_zt_coords')
xt = real(ztt);
yt = imag(ztt);

plot(inx, iny,'o-','markersize',8)
hold on; grid on
plot(xt,yt,'k.-','markersize',8)
plot(x,  y,'mx')

xlim([100,250])
ylim([100,250])

legend('(x_{in} y_{in})','(x_t,y_t)','(x,y)','location','SouthEast')

title (['Odd N = ',num2str(N)]);

这是 ODD 案例的结果: 参数奇数情况

这是 180 $z_t$坐标 的 .mat 文件的链接: https ://www.dropbox.com/s/gifbbvyfl0unv3f/pi_zt_coords.zip?dl=0

4个回答

我不明白评论。

你当然可以这样做。只需了解 DFT 的含义、如何计算 DFT bin 值以及如何将这些 bin 值解释为连续傅立叶级数系数​​。

首先,您正在查看的平面是复平面。您的点是一组$N$离散样本。每个样本都是一个复杂的点。因此,您所拥有的是重复复杂信号的一个周期的表示。图中的间距不是那么重要。

  • 任何$N$个点的序列都可以用N个系数精确地在样本点处表示

问题是:“你的形状是否允许你忽略系数,所以你的系数要少得多?”

答案是:“取决于形状。” 因此,开始丢弃最小的幅度系数,看看精度会受到多大影响。

当您构建傅里叶级数时,您需要将未归一化的 DFT 系数除以$N$您还希望将 DFT 的上半部分算作负频率,因此$N-1$对应于$-1$等。

所以基本上你是在对离散序列进行 DFT,然后使用系数重建插值。

希望这可以帮助。

赛德


我把这个放在我的答案中,因为我不想触发将此对话转移到聊天空间(我不同意 BTW 的政策)。

手头的问题不仅仅是一个封闭的图形是否可以参数化,这个问题被牢固地设置为本轮求和的应用(你知道,在哥白尼改变参考系之前,行星运动是如何被建模的)。是的,还有其他方法可以参数化圆周运动而不是正弦和余弦,但它们很笨拙。

还有其他方法可以退后一步并将图形作为一个整体参数化,并且不需要它是周期性的。我想到了勒让德基础。DFT 方法就是这样发生的,它本质上是周期性的。

在我看来,OP 认为可以绘制任意图形(在限制内)很酷(我也是),并试图理解本轮的概念与 DFT 的关系。

让我们做一些数学运算以使其更清楚。使用传统的归一化和符号,DFT 是:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-ink \frac{2\pi}{N} } $$

由于$x[n]$是已知的,$X[k]$现在也是已知的。现在,让我们看一下反面:

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{ink \frac{2\pi}{N} } $$

如果我们简单地允许$n$为实值并将逆 DFT 定义视为一个连续方程,我们会遇到超过奈奎斯特频率的问题。在离散情况下没有区别,因为它们将在样本点处匹配。在这两者之间,确实如此。因此,总和必须移动到以 DC bin 为中心。(假设$N$是偶数,否则可以类似地计算出来)

$$ x[n] = \frac{1}{N} \sum_{k=-N/2+1}^{N/2} X[k] e^{ink \frac{2\pi}{N } } $$

该方程也可以分为实部和虚部:

$$ \Re(x(n)) = \sum_{k=-N/2+1}^{N/2} X[k] \frac{1}{N} \cos( nk \frac{2\ pi}{N} ) $$ $$ \Im(x(n)) = \sum_{k=-N/2+1}^{N/2} X[k] \frac{1}{N} \罪(NK \frac{2\pi}{N})$$

我会使用$x$$y$,但$x$已经被占用了。这些方程显然是傅里叶级数的形式,系数为$ X[k] / N $

我不是想在这里教育 r bj,我知道他对这些东西了如指掌。我只是说引入替代参数化或替代坐标系会分散手头的核心问题。


是的,刚刚接受 r bj 的教育。也感谢您的编辑。

事实上,奈奎斯特项应该分成两半,结果是这两个本轮将相互抵消虚部并加倍实部。由于图中没有很多曲折,我预计这个系数的大小会很低。

只是为了好玩,我写了一个小 Gambas 程序来演示数学。你可以在这里找到它:

https://forum.gambas.one/viewtopic.php?f=4&t=725

我还做了一点手绘的 Pi 符号。当然,它看起来有点醉了,但它仍然证明了这一点。

在此处输入图像描述


根据要求,这里有一点角落处理。角落的工作比我预期的要好。我认为这个例子真正体现了我之前所说的真正有趣的问题是沿着图形找到最接近的点放置。

在此处输入图像描述


在这个讨论中我没有明确说的是$ e^{i\theta} $的复数值沿着复单位圆移动,因此是一个本轮模型,所以循环内的每个产品代表如果要制作本轮动画,则在该时间点其各自本轮(即线段)的半径位置。半径的长度是系数的大小,因为$ e^{i\theta} $的大小始终为 1。

Complex.Polar(r, theta) = r * e ^ { i theta }

您可能会发现我的这篇文章有助于理解这些材料:

我不是 MATLAB 的粉丝(主要是因为对基于数组的使用极其短视),所以我不会评论您的伪代码。相反,这是我的代码,它实际计算给定“n”处的插值。

你可以点击链接自己下载(我只是放了新版本,允许在同一张图中有多个图形)。如果你有 Linux,你可以安装 Gambas (PPA:gambas-team/gambas3) 来运行它。

[注:代码中的n只是一个迭代器,t是实际的n,我懒得编辑代码。]

.
.
.

        对于 n = 0 到 myPoints.Count * 100 - 1
          t = n / 100  
          p = 计算(t, w)
          Paint.Arc(p.Real, p.Imag, 1)
          Paint.Fill()
        下一个
.
.
.

'================================================== =======================
Public Sub Calculate(ArgT As Float, ArgDFT As Vector) As Complex

        Dim k, N 作为整数
        暗淡如复杂
        将 a、b 调暗为浮点数

        N = ArgDFT.Count

        b = ArgT * Pi(2) / N

        如果偶数(N)然后
           GoSub EvenCase
        别的    
           GoSub OddCase
        万一

        返回 p

'---------------------------------------------------- ----------------------
偶数情况:

        p = ArgDFT[0] + ArgDFT[N / 2] * Cos(ArgT * Pi)

        对于 k = 1 到 N / 2 - 1
          a = b * k
          p += ArgDFT[k] * Complex.Polar(1.0, a)
          p += ArgDFT[N - k] * Complex.Polar(1.0, -a)
        下一个

        返回

'---------------------------------------------------- ----------------------
奇案:

        p = ArgDFT[0]

        对于 k = 1 到 (N - 1) / 2
          a = b * k
          p += ArgDFT[k] * Complex.Polar(1.0, a)
          p += ArgDFT[N - k] * Complex.Polar(1.0, -a)
        下一个

        返回

结尾
'================================================== =======================

对 r bj 的反驳:

罗伯特,我非常不同意你的一些主张。

1) 插值点(以及它们通过 LineTo 调用形成的路径)将遵循您提供点的任何顺序

2)与我的答案相比,“t”使用 0 到$2\pi$范围会混淆问题,其中“t”的范围从 0 到 N,即与离散尺度相同的尺度,仅包括介于两者之间的实际值整数。这反映了您的参考框架是连续案例。[不再相关,我已经放弃了 t]

3) 将 (x,y) 视为一个向量,而不是单个复数值 x + iy,将参数化分为两个独立的问题,这两个问题不必通过相同的方法进行参数化。只有在复杂的价值解释中,作为这个问题核心的本轮概念才有意义。

4) 将拐角处的点捆绑在一起,中间没有足够数量的点,会导致拐角处的超限。这就是我在上一张图中添加第四个数字的原因。

5)您对$a_k$$b_k$的定义没有意义,因为没有给出连续函数,只有一组样本点。因此,傅立叶系数应该使用离散定义来计算,即求和而不是积分。你本末倒置。使用不同的参数化,例如勒让德,您不会有超出范围的重复模式,不一定会在点之间匹配,但您会在所有样本点处匹配。

标题问题是:“如何获得傅立叶系数以使用 DFT 绘制任何形状?”

答案是:“归一化的 DFT bin 值是傅立叶系数。”

换句话说,简单地将中心逆 DFT 中的整数“n”替换为连续实值变量将产生插值结果。没有比这更优雅的了。我的代码就是对此的一种表达。我假设 OP 将在 MATLAB 中实现它(进行必要的索引调整)。

你让这种方式比它需要的更复杂。


这是偶数情况循环的等效编码,以阐明“k”的含义及其范围。

       对于 k = -N / 2 + 1 到 -1
          a = b * k
          p += ArgDFT[k+N] * Complex.Polar(1.0, a)
       下一个

       对于 k = 1 到 N / 2 - 1
          a = b * k
          p += ArgDFT[k] * Complex.Polar(1.0, a)
       下一个


这个是给 Olli 的,使用 N = 9。如果图形确实是一个三角形,你可以看到,通过适当的点放置,可以找到一个更好的拟合,它也能击中所有的点。当然,可以添加更多点(更多本轮)来获得更紧密的拟合。

在此处输入图像描述

这里的辅助问题(前面提到过),我认为这将是你要解决的问题,Olli,是如何将样本点放置在底层连续图形上以最小化“涟漪”或“溢出”。


当您将 DFT 的上半部分视为正频率而不是负频率时,就会发生这种情况。您可以清楚地看到所有点都被击中,但在结果之间并不是我们想要的。也许有一些新颖的应用程序会有所帮助。

在此处输入图像描述

我这样做是为了回应 Olli 提出的是否只能用正频率来完成的挑战。也许,如果将实部和虚部分开并且将余弦级数用于两个参数化,但我认为这会阻碍问题的意图,并且它不再是本轮实现。

我最初的倾向是拒绝。我认为这个问题相当于“你能从顺时针开瓶器的总和中构造一个逆时针开瓶器吗?” 也许有一个无穷大的数字,我在数学中看到了太多关于接近无穷大的奇怪事情来排除它,但在这里我什至无法想象一个近似的序列。


克里斯,

我不反对变量的大写字母。事实上,我喜欢使用$S_n$作为我的信号值。我不喜欢对信号使用小写“x”,对 DFT 使用大写“X”。对我来说,这还不够,因为它们描述的是两个完全不同的领域。此外,“X”是小写和大写最相似的字母之一,这使得它们在手写数学中更难区分。

在这种情况下,我们基本上有三个尺度(或功能域):

1) n 在样本点的整数上从 0 变为 N-1(对于输入点和逆 DFT 的输出)

2) k 在逆 DFT 定义中的整数上从 0 变为 N-1,然后移动半帧以消除上半部分的锯齿

3) t 从 0 变为$2\pi$是系列解决方案的域(你和罗伯特)和(0 到 N - 我的代码中的 1/100)

所以,是的,你在 T 域中使用 K 是在误导。

在我的代码中,ArgDFT 是 1/N 归一化 DFT,而 ArgT 是我原来的“t”参数,它与“n”的比例相同,但是是连续的。我在代码中的“b”和你的“t”一样。

流程总结:

当您对样本序列进行 1/N 归一化 DFT 时,您会同时找到将通过所有点的连续傅里叶级数的系数。(为什么应该使用 1/N 归一化的一个强有力的论据)。

级数解的域可以通过变量替换重新缩放:

$$ n = t \cdot \frac{N}{2\pi} $$

转化为被解释为连续函数的逆 DFT。

$$ x(n) = \frac{1}{N} \sum_{k} X[k] e^{ink \frac{2\pi}{N} } $$

$$ z(t) = x(t \cdot \frac{N}{2\pi}) = \frac{1}{N} \sum_{k} X[k] e^{it \cdot \frac{ N}{2\pi} k \frac{2\pi}{N} } $$

$$ z(t) = \sum_{k} \frac{X[k]}{N} e^{ikt } $$

这就是连续插值路径的级数解。它只是 t 的函数。如果需要,您可以将其与 t 区分开来计算您的“笔速度”。

很明显,您现在已经理解了我所说的“点放置问题”的含义,而且看起来 Olli 也对它产生了兴趣。

如果您还没有,我建议您重新阅读此线程中的所有内容。有了更好的理解基础,已经说过的话应该更有意义。


结语:关于情况的不同观点,这里的许多人都熟悉。但是,它不会产生傅立叶系数。

本来就是一堆烂摊子。

$$ z(t) = \sum_{k} \frac{X[k]}{N} e^{ikt } $$

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-ink \frac{2\pi}{N} } $$

$$ z(t) = \frac{1}{N} \sum_{k} \sum_{n=0}^{N-1} x[n] e^{-ink \frac{2\pi}{ N} } e^{ikt } $$

$$ z(t) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] \sum_{k} e^{ik ( t - \frac{n}{ N}2\pi )} $$

$$ t_n = t - \frac{n}{N}2\pi $$

$$ D(t_n) = \sum_{k} e^{ik t_n } $$

$$ z(t) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] D(t_n) $$

奇数情况:$k = -(N-1)/2 \to (N-1)/2$

$l = k + (N-1)/2$变为$0 \到 N-1$

$$ k = l - (N-1)/2 $$

$$ \begin{aligned} D(t_n) &= \sum_{l=0}^{N-1} e^{i ( l - (N-1)/2 ) t_n } \\ &= \sum_{ l=0}^{N-1} e^{il t_n } e^{-i \frac{N-1}{2} t_n } \\ &= e^{-i \frac{N-1}{ 2} t_n} \sum_{l=0}^{N-1} (e^{i t_n })^l \\ &= e^{-i \frac{N-1}{2} t_n} \frac {1 - (e^{i t_n })^N }{ 1 - e^{i t_n } } \\ &= e^{-i \frac{N-1}{2} t_n} \left[ \frac {e^{i t_n N / 2 } } { e^{i t_n / 2 } } \cdot \frac{ e^{-i t_n N / 2 } - e^{i t_n N/2 } }{ e^ {-i t_n / 2 } - e^{i t_n / 2 } } \right] \\ &= \frac{e^{i t_n N / 2 } - e^{-i t_n N / 2 }} { e ^{i t_n / 2 } - e^{-i t_n / 2 } } \\ &= \frac{ 2i \cdot \sin( N t_n / 2 ) } { 2i \cdot \sin( t_n / 2 ) } \ \ &= \frac{ \sin( N t_n / 2 ) } { \sin( t_n / 2 ) } \end{对齐} $$

$$ z(t) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] \frac{ \sin( N t_n / 2 ) } { \sin( t_n / 2) } $$

$$ z(t) = \sum_{n=0}^{N-1} x[n] \frac{ \sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \sin( (t - \frac{n}{N}2\pi) / 2 ) } $$

请注意,商是实数值,因此可以将其视为权重,然后求和是样本点集的时变加权平均值。


结语二

在对其他问题进行了大量讨论之后,很明显 Nyquist bin 应该在负频率和正频率之间平均分配。

偶数情况:$k = 1/2 ( N/2 \text{ and } -N/2 ), -N/2 + 1 \to N/2 - 1 $

$l = k + N/2 - 1 $$0 \到 N-2$

$$ k = l - N/2 + 1 $$

$$ \begin{aligned} D(t_n) &= \frac{1}{2} \left[ e^{i ( N/2 ) t_n } + e^{i (-N/2 ) t_n } \right ] + \sum_{l=0}^{N-2} e^{i ( l - N/2 + 1 ) t_n } \\ &= \cos \left( \frac{N}{2} t_n \right ) + \sum_{l=0}^{N-2} e^{il t_n } e^{i (- N/2 + 1 ) t_n } \\ &= \cos \left( \frac{N}{ 2} t_n \right) + e^{i (- N/2 + 1 ) t_n } \sum_{l=0}^{N-2} (e^{i t_n })^l \\ &= \cos \left( \frac{N}{2} t_n \right) + e^{i (- N/2 + 1 ) t_n } \frac{1 - (e^{i t_n })^{N-1} } { 1 - e^{i t_n } } \\ &= \cos \left( \frac{N}{2} t_n \right) + e^{i (- N/2 + 1 ) t_n } \left[ \ frac{e^{i t_n ( N - 1 ) / 2 } } { e^{i t_n / 2 } } \cdot \frac{ e^{-i t_n ( N - 1 ) / 2 } - e^{i t_n ( N - 1 ) / 2 } }{ e^{-i t_n / 2 } - e^{i t_n / 2 } } \right] \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{e^{i t_n ( N - 1 ) / 2 } - e^{-i t_n ( N - 1 ) / 2 } } { e^{i t_n / 2 } - e^{-i t_n / 2 } } \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{ 2i \cdot \ sin( t_n ( N - 1 ) / 2 ) } { 2i \cdot \sin( t_n / 2 ) } \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{ \ sin( t_n N /2 ) \cos( t_n / 2 ) - \cos( t_n N /2 ) \sin( t_n / 2 ) } { \sin( t_n / 2 ) } \\ &= \cos \left( \ frac{N}{2} t_n \right) + \frac{ \sin( t_n N /2 ) } { \sin( t_n / 2 ) } \cos( t_n / 2 ) - \cos( t_n N /2 ) \ \ &= \frac{ \sin( N t_n/2 ) }{ \sin( t_n / 2 ) } \cos( t_n / 2 ) \end{对齐} $$= \cos \left( \frac{N}{2} t_n \right) + \frac{ \sin( t_n N /2 ) \cos( t_n / 2 ) - \cos( t_n N /2 ) \sin( t_n / 2 ) } { \sin( t_n / 2 ) } \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{ \sin( t_n N /2 ) } { \sin ( t_n / 2 ) } \cos( t_n / 2 ) - \cos( t_n N /2 ) \\ &= \frac{ \sin( N t_n/2 ) }{ \sin( t_n / 2 ) } \cos( t_n / 2 ) \end{对齐} $$= \cos \left( \frac{N}{2} t_n \right) + \frac{ \sin( t_n N /2 ) \cos( t_n / 2 ) - \cos( t_n N /2 ) \sin( t_n / 2 ) } { \sin( t_n / 2 ) } \\ &= \cos \left( \frac{N}{2} t_n \right) + \frac{ \sin( t_n N /2 ) } { \sin ( t_n / 2 ) } \cos( t_n / 2 ) - \cos( t_n N /2 ) \\ &= \frac{ \sin( N t_n/2 ) }{ \sin( t_n / 2 ) } \cos( t_n / 2 ) \end{对齐} $$

对于正负奈奎斯特项,上述推导可以使用 1/2 和 1/2 以外的系数来完成,但最后发生的简化不会发生,表达式会更复杂。同样非常清楚的是,如果$x[n]$的集合是真实的,则插值不一定是真实的。对于 1/2 和 1/2,插值都是实数。

则连续插值函数为:

$$ \begin{aligned} z(t) &= \frac{1}{N} \sum_{n=0}^{N-1} x[n] \left[ \frac{ \sin( N t_n / 2 ) }{ \sin( t_n / 2 ) } \right] \cos( t_n / 2 ) \\ &= \sum_{n=0}^{N-1} x[n] \left[ \frac{ \ sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \sin( (t - \frac{n}{N}2\pi) / 2 ) } \right] \ cos( (t - \frac{n}{N}2\pi) / 2 ) \\ &= \sum_{n=0}^{N-1} x[n] \frac{ \sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \tan( (t - \frac{n}{N}2\pi) / 2 ) } \end{aligned} $$

值得注意的是,该公式将奇数情况版本与前两行中应用的简单“窗口函数”相匹配。最后一个匹配 R BJ 的给定公式,格式更简洁。

查看 N = 2 的情况

$$ \begin{aligned} z(t) &= x[0] \left[ \cos^2( t / 2 ) \right] + x[1] \left[ \cos^2( (t - \pi ) / 2 ) \right] \\ &= x[0] \left[ \frac{ \cos( t ) + 1 }{2} \right] + x[1] \left[ \frac{ \cos( t - \pi ) + 1 }{2} \right] \\ &= \frac{1}{2} ( x[0] + x[1] ) + \frac{1}{2} ( x[0] - x[1] ) \cos( t ) \end{对齐} $$

这意味着交替序列 1、-1、1、-1 被插值为:

$$ z(t) = \cos(t) $$

这是奈奎斯特分裂假设的直接结果。

跟踪所需形状的分段线性波形的复傅里叶级数

代替使用离散傅里叶变换(DFT)/快速傅里叶变换(FFT),更直接的方法是定义一个分段线性连续时间波形,在复平面上追踪所需形状,并直接计算其傅里叶级数. 贝塞尔曲线等可用于形状定义并使用线段近似到任意精度。您的第三个图形已经使用线段绘制。我们可以使用它的节点(角)坐标,但是波形节点的时间需要从帽子里拉出来。我们将及时进行统一采样,而无需在数学或脚本中进行硬编码。波形可以在 Octave 中绘制,将来自此答案pi.csv的文件中的节点坐标作为输入

graphics_toolkit("gnuplot")  # Octave specific to get prettier plots

xy = csvread("pi.csv");
z = xy(:,1) + i*xy(:,2);
M = length(xy);
t = (0:M-1)'*2*pi/M;

plot([t; 2*pi], [real(z); real(z(1))], "b");
hold on
plot([t; 2*pi], [imag(z); imag(z(1))], "r");
plot(t, real(z), "k.");
plot(t, imag(z), "k.");
xlim([0,2*pi])
ylim([-250,250])
xlabel("t")
hold off

![在此处输入图像描述
图 1. 跟踪所需形状的分段线性波形的实部(蓝色)和虚部(红色)。

波形的复数傅里叶级数

让我们看一下波形的单个线性段。一个$2\pi$ - 周期连续时间波形,否则为零,但有一条线段,从复数值$z_0 = x_0 + y_0\,i$在时间$t_0$开始,并以值$z_1 = x_1 + y_1\结束,i$在时间$t_1 > t_0$具有其复数傅立叶级数的系数(使用您的第二个方程):

$$c_k=\frac{1}{2\pi}\int_\limits{t_0}^{t_1}\left(z_0 + \frac{t-t_0}{t_1-t_0}(z_1 - z_0)\right) \, e^{-ikt} \mathop{dt},\tag{1}$$

其中$\frac{t-t_0}{t_1-t_0}$$0$$1$因为$t$$t_0$$t_1$对于$k = 0$,我们有:

$$c_0 = \frac{(t_1 - t_0)(z_0 + z_1)}{4\pi},\tag{2}$$

对于负数和正数$k \ne 0$

$$\begin{gather}c_k = \color{blue}{\frac{z_1\,\sin(k\,t_1) - z_0\,\sin(k\,t_0)}{2\pi k}} + \frac{(z_1 - z_0) \cos(k\,t_1) - (z_1 - z_0) \cos(k \,t_0)}{2\pi k^2\,(t_1 - t_0)}\,+\ \ i\,\left(\color{blue}{\frac{z_1\,\cos(k\,t_1) - z_0\,\cos(k\,t_0)}{2\pi k}} - \frac {(z_1 - z_0)\,\sin(k\,t_1) - (z_1 - z_0) \sin(k\,t_0)}{2\pi k^2\,(t_1 - t_0)}\right)。 \end{聚集}\tag{3}$$

您将完全填充范围$0 \le t \le 2\pi$(或任何长度范围$2\pi$,例如$-\pi \le t \le \pi$与您的第二个方程兼容)与非-重叠线性段,并且对于每个整数$k$,分别使用这些段上的系数之和来获得完整分段线性波形的复数傅里叶级数的相应系数。这是可行的,因为频域中的加法等同于将线性段拼接在一起的时域中的加法。方程式中以蓝色着色的术语。3 将抵消波形段的总和,不需要包含在其中。使用您的第一个方程,完整的波形等于其复杂的傅立叶级数:

$$z(t) = \sum_{k=-N}^N c_k\,e^{ikt},\tag{4}$$

$N=\infty$除系数为 $c_0$的常数项外,和的各个调和项在复平面上绕圈,可视为本轮。

您可以在某个有限的$N$处截断系列。这是实现此方法的 Octave 脚本。它计算由给定节点(角)定义的分段线性波形的截断复傅立叶级数,假设节点时间分布均匀:

graphics_toolkit("gnuplot")  # Octave specific to get prettier plots

xy = csvread("pi.csv");
z = xy(:,1) + i*xy(:,2);

M = length(xy);
N = floor(M/2) - 1;  # Truncation length, this can be any positive integer
k = -N:N;
t = (0:M-1)'*2*pi/M;  # This can be any ascending sequence of times of the nodes obeying 0 <= t < 2 pi
t1 = circshift(t,-1);
t1(end) = 2*pi;
z1 = circshift(z,-1);

c = sum(((z1 - z).*cos(k.*t1) - (z1 - z).*cos(k .*t))./(2*pi*k.^2.*(t1 - t))+ i*(-((z1 - z).*sin(k.*t1) - (z1 - z).*sin(k.*t))./(2*pi*k.^2.*(t1 - t))), 1);
c(N + 1) = sum(((t1 - t).*(z + z1))/(4*pi), 1);

# c now contains complex Fourier series coefficients in order k

z_new = (2*N + 1)*ifft(ifftshift(c));  # Uniformly sample the reconstruction in time
xy_new = [real(z_new)', imag(z_new)'];
csvwrite("pi_bandlimited.csv", xy_new);  # Save samples. This should work with rb-j's script for odd length

os = 8;  # Oversampling factor, integer
z_os = os*(2*N + 1)*ifft([c(N+1:end) zeros(1, (N*2+1)*(os - 1)) c(1:N)]);  # Band-limited approximation
plot([real(z_os) real(z_os(1))], [imag(z_os) imag(z_os(1))], "-")
xlim([-250,250]);
ylim([-250,250]);
hold on
#plot(real(z_new), imag(z_new), "+")  # New samples
#plot([real(z);real(z(1))], [imag(z);imag(z(1))], "-")  # Desired shape
hold off

结果(图 2)可以在视觉上与傅立叶插值(由@robertbristow-johnson 的答案中提出的基于 DFT 的方法给出)(图 3)进行视觉比较,继续上述 Octave 脚本:

z_ftos = interpft(z, length(z)*os);  # Fourier interpolate
plot([real(z_ftos); real(z_ftos(1))], [imag(z_ftos); imag(z_ftos(1))], "-")
xlim([-250,250]);
ylim([-250,250]);

在此处输入图像描述
图 2. 此答案中建议的方法的结果是分段线性波形跟踪所需形状的最小二乘带限近似,此处使用 119 个谐波项。

在此处输入图像描述
图 3. 傅里叶插值方法的结果(本答案中未提供),使用 120 个谐波项。

从上面可以看出,建议的方法可以更清晰地跟踪所需的形状,并且可以通过增加$N$(图 4)轻松地使其更准确。

在此处输入图像描述
图 4. 建议的方法设置的结果N = 2000将复杂的傅立叶级数截断为 4001 个谐波项,并给出了所需形状的这种视觉上无法区分的近似值。

如果你愿意,你可以重新排列等式的总和。4 交错正负系数,或通过减少$|c_k|$对其进行排序。您还可以通过仅选取最大系数来构建稀疏近似。我们可以通过运行上面的 Octave 脚本来查看谐波的大小N = 20000并绘制:

loglog(abs(k), abs(c), '.');
xlim([1, 20000]);
ylim([0.000001, max(abs(c))]);
xlabel("|k|");
ylabel("|c_k|")

在此处输入图像描述
图 5. 分段线性波形的谐波幅度跟踪$\pi$形状。由于波形的连续性,包络以 -40 dB / 十倍频斜率渐近衰减。

波形近似误差

波形近似中的均方误差将随着复数傅里叶级数的每个包含项而下降,无论其顺序如何。这是因为谐波项是正交的,因此项的任何部分和的均方是项的均方的和,即:

$$\frac{1}{2\pi}\int_\limits{0}^{2\pi}\left|c_k\,e^{ikt}\right|^2\mathop{dt} = |c_k| ^2.\标签{5}$$

波形的均方等于复数傅立叶级数的均方,它是谐波项的部分和的极限为$N\to\infty$,并且可以在时域中等效地计算为均方和的线性段,它们是:

$$\frac{1}{2\pi}\int_\limits{t_0}^{t_1}\left|z_0 + \frac{t-t_0}{t_1-t_0}(z_1 - z_0)\right|^2 \mathop{dt} = \frac{(t_1 - t_0)(z_0^2 + z_0 z_1 + z_1^2)}{6\pi}.\tag{6}$$

两个和之间的差异是波形近似的均方误差,然而,它不是唯一定义的,甚至在跟踪所需形状时也不是一个合理的误差度量。

形状的最佳跟踪

分段线性复数波形的截断复数傅里叶级数通常不是最接近所需形状的最佳$2N+1$项近似。可以构建一组新的线段,它是复杂波形的时间拉伸版本,将跟踪相同的形状,但可以产生更好的截断序列。如果线段的数量足够增加,我认为它可以在某种意义上任意接近该形状的最佳波形。找到最佳波形似乎是一个难以解决的问题。

输入数据文件

此答案用于存储输入数据文件以测试问题的解决方案。

这是@Chris提供的$\pi$形状的120 个$x,y$坐标。另存为pi.csv

108,0
110,25
112,50
113.5,75
115,100
116,125
117.5,150
125,150
150,150
175,150
200,150
225,150
225,175
225,200
225,220
200,220
175,220
150,220
125,220
100,220
75,220
50,220
25,220
0,219.5
-25,219
-50,217
-75,215
-100,212
-125,209
-150,203
-158,200
-175,190
-190,175
-203,150
-211,125
-220,100
-225,85
-209,85
-200,100
-182,125
-175,132
-150,145
-125,150
-100,150
-87,150
-87.5,125
-89,100
-92,75
-95,50
-100,25
-105,0
-113,-25
-122,-50
-136,-75
-152,-100
-170,-125
-186,-150
-189,-175
-178,-200
-175,-205
-150,-220
-125,-220
-100,-202
-85,-175
-77,-150
-73,-125
-70,-100
-67.5,-75
-65,-50
-62,-25
-60,0
-57,25
-54.5,50
-51.5,75
-49,100
-47,125
-45,150
-25,150
0,150
25,150
50,150
58,150
55,125
53,100
51,75
49,50
47,25
44.5,0
42,-25
40,-50
38.5,-75
37.5,-100
37,-125
37.5,-150
43,-175
49,-185
66,-200
75,-205
100,-215
125,-218
150,-214
175,-203
179,-200
201.5,-175
213,-150
221,-125
226.5,-100
227.5,-88
210,-88
209,-100
200,-123
197,-125
175,-141
150,-144
125,-134
117,-125
109,-100
106,-75
106,-50
106.5,-25

我已经将@robert bristow-johnson的第二个Matlab程序改编为python,有些人可能会觉得它很有用。我使用了来自@Olli Niemitalo 的数据,但对其进行了扩展以提供指向(x,y)原点的返回点。这允许零填充大大减少 - 从 65536 到 300 - 基本相同的结果,见下文。

在此处输入图像描述 非最佳 python 列表是:

# Python version adapted from Matlab program by Robert Bristow-Johnson
# URL: https://dsp.stackexchange.com/questions/59068/how-to-get-fourier-coefficients-to-draw-any-shape-using-dft

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft
import pandas

df = pandas.read_csv('pi.csv') # data from Olli Niemitalo

xx = np.array(df['x'],dtype='int')
yy = np.array(df['y'],dtype='int')

xx = np.append(xx,xx[0]) # to complete figure
yy = np.append(yy,yy[0]) # to complete figure

NN = len(xx)
NN2 = int(NN/2)
N = 300 # must be greater than NN

XX = fft(xx)/NN
YY = fft(yy)/NN

X0 = np.append(XX[range(NN2)], np.zeros([N-NN],dtype=np.complex))
X  = np.append(X0, XX[range(NN2,NN)])

Y0 = np.append(YY[range(NN2)], np.zeros([N-NN],dtype=np.complex))
Y  = np.append(Y0, YY[range(NN2,NN)])

x = np.real(N*ifft(X)) # real values taken for plotting
y = np.real(N*ifft(Y))

fig1 = plt.figure(figsize=(16,6))
ax1 = fig1.add_subplot(131)
ax1.plot(xx, yy, 'mx')
ax1.plot(x, y, 'b')
ax1.set_title(f"FFT fit of $\pi$ figure - zero-padding: N ={N:3}")

xnn = np.linspace(0,NN,NN)
xn  = np.linspace(0,NN,N)

#
ax2 = fig1.add_subplot(132)
ax2.plot(xnn,xx, 'mx')
ax2.plot(xn, x, 'b')
ax2.set_title("FFT fit of x-coordinates")

#
ax3 = fig1.add_subplot(133)
ax3.plot(xnn,yy, 'mx')
ax3.plot(xn, y, 'b')
ax3.set_title("FFT fit of y-coordinates")