C：数值计算 (FFT)

Question

这个问题是针对任何 Numeric Recipes 的粉丝或任何熟悉 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)。我可以想到更多一些简单地将余弦和-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

Answer 1

起始：

• 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²(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) ] + i sin(φ) [ 1 - 2 sin ²(δ/2) ] + i sin(δ) [ i * sin(φ) + cos(φ) ]

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

= e^{i φ} + e^{i φ [ - 2 sin2(δ/2) + i sin(δ)]}

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

对于任何y，y^(a+b) = y^a × y^b。因此：

e^{i (φ+δ)}

= e^{i φ} e^{i δ}

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

= e^{i φ} [ 1 - 2 sin²(δ/2) + i sin(δ) ]

Answer 2

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

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

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

Answer 3

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

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。我的实验表明，对于 float cos(2*PI/N) == 1 来说，对于 float （即 32 位精度）来说，N >= 25736 正是如此。 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

Answer 4

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

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

并根据 sin(u)*sin(v) 的第一个“乘积求和”恒等式，我们有

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

这有帮助吗？

Answer 5

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

Answer 6

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

顺便说一句，众所周知，NR 的实现速度非常慢。

保罗

Answer 7

好的，这是三角恒等式。它不是半 cos(theta) 三角恒等式的原因是因为依赖欧拉和虚数。数学仍然超出我的能力......

_{（来源：librow.com)}