4

这个问题是针对任何数字食谱的粉丝或任何理解 FFT 的人。

谁能解释为什么实际分量是由 -2*(sin(theta/2))^2 计算的？我似乎无法绕过它。我看过其他示例，例如http://www.dspdimension.com/admin/dft-a-pied/教程，它只是将 cos(theta) 视为实数，将 -sin(theta) 视为虚数。我在基本的http://www.dspguide.com/ch12/3.htm中也看到过，它将 cos(theta) 列为实数，将 -sin(theta) 列为虚数。我可以想到更多的资源，它们只是将 cos 和 -sin 视为真实和想象。

cos(theta) = 1-2*(sin(theta/2))^2

如果上述三角恒等式是真的，那为什么不遵循呢？

theta=isign*(6.28318530717959/mmax);
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);

我假设数字配方必须使用一些三角标识？我似乎无法弄清楚，这本书根本没有解释。

代码在这里找到：http ://ronispc.chem.mcgill.ca/ronis/chem593/sinfft.c.html

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(double *data,unsigned long nn,int isign)
{
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        double tempr,tempi;

        n=nn << 1;
        j=1;
        for (i=1;i<n;i+=2) {
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        mmax=2;
        while (n > mmax) {
                istep=mmax << 1;
                theta=isign*(6.28318530717959/mmax);
                wtemp=sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=sin(theta);
                wr=1.0;
                wi=0.0;
                for (m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*data[j]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
}
#undef SWAP
4

7 回答 7

3

从...开始：

  • cos(A+B) = cos(A) cos(B) - sin(A) sin(B)
  • sin(A+B) = sin(A) cos(B) + cos(A) sin(B)
  • cos(2A) = 1 - 2 sin 2 (A)
  • e i θ = cos(θ) + i sin(θ)

所以：

e i (φ+δ)

= cos(φ + δ) + i sin(φ + δ)

= cos(φ) cos(δ) - sin(φ) sin(δ) + i [sin(φ) cos(δ) + cos(φ) sin(δ)]

= cos(φ) [ 1 - 2 sin 2 (δ/2) ] + i sin(φ) [ 1 - 2 sin 2 (δ/2) ] + i sin(δ) [ i * sin(φ) + cos (φ)]

= [ cos(φ) + i sin(φ) ] [ 1 - 2 sin 2 (δ/2) ] + [ cos(φ) + i sin(φ) ] i sin(δ)

= e i φ + e i φ [ - 2 sin 2 (δ/2) + i sin(δ)]

编辑：对我来说，这是很多无用的格式。它实际上更简单：

y (a+b) = y a × y b对于任何y。所以：

e i (φ+δ)

= e我 φ e我 δ

= e i φ [ cos(δ) + i sin(δ) ]

= e i φ [ 1 - 2 sin 2 (δ/2) + i sin(δ) ]

于 2010-02-11T01:24:41.927 回答
1

余弦半角恒等式的一种形式是：

cos(theta) = 1 - 2*(sin(theta/2)^2)

不确定这是否回答了你的问题。

于 2010-02-08T10:30:47.643 回答
1

原因是数值精度。如果您仔细检查以下代码：

wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);

和

wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;

它们旨在协同工作以产生正确的预期结果。一种天真的方法将按如下方式实施：

wpr = cos(theta);
wpi = sin(theta);

和

wtemp = wr;
wr =wr*wpr - wi*wpi;
wi =wi*wpr + wtemp*wpi;

并且具有无限精度在功能上是等效的。

然而，当theta接近于零（即很多样本点或大nn）时，cos(theta)就会出现问题，因为小角度cos(theta)接近 1 并且斜率接近 0。并且在某个角度cos(theta) == 1。我的实验表明 for floatcos(2*PI/N) == 1完全适用于N >= 25736for float （即 32 位精度）。一个 25,736 点的 FFT 是可以想象的。因此，为了避免这个问题wpr，计算为cos(theta)-1使用三角学中的半角公式。对于小角度，它的计算sin是非常精确的，因此对于小角度来说，两者wpr都是wpi小而精确的。然后通过在复数乘法后重新添加 1 在更新代码中撤消此操作。用数学表示，我们得到：

w_p = cos(theta) - 1    + i*sin(theta) 
    = -2*sin^2(theta/2) + i*sin(theta)

更新规则是：

w = w * (w_p + 1) 
  = w + w*w_p
于 2017-10-08T18:14:14.517 回答
0

我不太了解 FFT，我做过一个，但已经很久了。

所以我在http://www.sosmath.com/trig/Trig5/trig5/trig5.html查找了三角身份

从 sin(u)*sin(v) 的第一个“Product-To-Sum”身份，我们有

sin^2(theta/2) = (1/2) [cos(零) - cos(theta)] = 0.5 - 0.5 cos(theta)

这有帮助吗？

于 2010-02-08T10:35:45.510 回答
0

他们正在使用三角恒等式来最小化他们需要计算的循环函数的数量。许多快速实现只是查找这些函数。如果你真的想知道，你只需要通过展开上面的循环并进行适当的变量替换来计算细节......很乏味是的。

顺便说一句，众所周知，NR 的实施非常缓慢。

保罗

于 2010-02-10T03:44:42.990 回答
0

好的，这是三角标识。它不是半 cos(theta) 三角恒等式的原因是因为依赖于欧拉和虚数。数学还是超出了我...

链接文本
（来源：librow.com）

于 2010-02-11T00:55:57.647 回答
0

令人困惑的是，NR 使用 FFT 的数学/物理版本，它以与 EE 定义 FFT 的方式相反的方式旋转旋转因子。所以 NR 前向是 EE 逆向，反之亦然。请注意，在 NR 中，前向具有正指数而不是 EE 负指数。EE 方法将时间转换为频率，其中数学和物理版本使用角动量。

于 2014-01-18T19:57:29.773 回答