C: 数值计算方法（FFT）

Question

C: 数值计算方法（FFT）

4

这个问题是针对 Numerical 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) 列为虚部。我可以想到更多的资源，它们只是将 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);

我假设Numerical Recipe必须使用某些三角恒等式？我似乎无法弄清楚，而且这本书根本没有解释。

代码在这里：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

- Mac

今天晚些时候，我将比较这个结果和DSP指南，并报告哪个更准确。 - Mac

让数据正确执行...对我来说，为什么它不是与cos（theta）等效的三角恒等式仍然是个谜。 - Mac

1

本内容适用于 Numerical Recipes 爱好者或对 FFT 有良好理解的人群。两者不可兼得。任何对 FFT 有相当了解的人都会立即意识到 Numerical Recipes 中提供的信息已过时40年，因此不会成为其爱好者。 - J D

如果你正在寻找一个不错的现代FFT实现，那么请使用FFTW。http://www.fftw.org - J D

2

@Jon Harrop：完全同意；Numerical Recipes不仅非常古老，而且特别是C版本是“在没有对C的热爱或知识的情况下编写的”。作者们非常熟悉Fortran，虽然这确实导致了合理的Fortran，但绝对不会导致好的C。 - FrankH.

7个回答

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。我的实验表明，对于浮点数（即32位精度），当 N >= 25736 时，cos(2*PI/N) == 1。可以想象一个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

- Commodore63

1

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

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

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

- Mitch Wheat

我猜这是一个打字错误，他们漏掉了1吧？所以应该是：wpr = 1 -2.0wtempwtemp; - Mac

Mac：不确定。浏览整个代码将非常耗时... - Mitch Wheat

@MAc：我的 NR 副本也有相同的代码。你是否查看了 NR 的勘误网站？ - Mitch Wheat

@Mac：这很不可能是一个打字错误；否则会对结果产生太大的影响。 - Mitch Wheat

Mitch Wheat: 谢谢...我认为这可能是一个打字错误...我读到了一些其他代码中的打字错误...很奇怪，因为如果你在谷歌上搜索它，它会完全显示而没有任何更正？它必须不会显著改变结果。+1 为引导我走上正确的道路 - Mac

0

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

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

从sin(u)*sin(v)的第一条“积和式”恒等式中我们得到：

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

这有帮助吗？

- Paul

0

他们正在使用三角恒等式来最小化需要计算的圆函数数量。许多快速实现只需查找这些函数即可。如果您真的想知道，您需要通过展开上面的循环并进行适当的变量替换来解决细节问题...是有点繁琐的。

顺便说一句，NR实现被认为非常慢。

保罗

- Paul

@Paul - 谢谢Paul。我明白NR很慢，但我需要一些严谨的东西去理解。我发现大多数可用的教程/源代码在其代码中做出了太多的假设。尽管我不明白为什么它不是完整的三角恒等式，但我已经读到三角恒等式可以避免sin/cos之间的舍入误差。这是一个有趣的话题，但至少可以说它让我感到困惑。 - Mac

我理解，高级实现可能非常复杂！祝好。 - Paul

保罗，除了NR算法，哪里可以找到更快的算法？ - FreelanceConsultant

快速傅里叶变换库：英特尔的MKL库（C，C ++，Fortran）或MIT许可的FFTW（C，C ++），或CenterSpaces的NMath库（C＃，.NET）。 - Paul

0

好的，这里是三角恒等式。它不是半余弦（cos(theta)）三角恒等式的原因是因为它依赖于欧拉和虚数。这种数学对我来说仍然超出了范围...

_{(来源：librow.com)}

- Mac

0

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

- Jeff

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Alok Singhal · Accepted Answer

从以下公式开始：

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 sin²(δ/2) + i sin(δ)]

编辑：我在格式方面做了很多无用的工作。其实它要简单得多：

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

e^{i (φ+δ)}

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

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

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