C: Числовые реципиенты (БПФ)

Question

Этот вопрос адресован любым любителям Числовых Рецептов или любому, кто хорошо понимает БПФ.

Может кто-нибудь объяснить, почему реальный компонент вычисляется как -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 (тэта) как реальный и -sin (тэта) как мнимый.Я могу подумать о еще нескольких ресурсах, которые просто принимают 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

Answer 1

Начать с:

• cos (A + B) = cos (A) cos (B) - sin (A) sin (B)
• грех (A + B) = грех (A), потому что (B) + cos (A), грех (B)
• cos (2A) = 1 - 2 sin ² (A)
• e ^{i & theta;} = cos (& theta;) + i sin (& theta;)

Итак:

e ^{i (& phi; + & delta;)}

= cos (& phi; + & delta;) + i sin (& phi; + & delta;)

= cos (& phi;) cos (& delta;) - грех (& phi;) sin (& delta;) + i [sin (& phi;) cos (& delta;) + cos (& phi;) sin (& delta;)]

= cos (& phi;) [1 - 2 sin ² (& delta; / 2)] + i sin (& phi;) [1 - 2 sin ² (& delta; / 2)] + i sin (& delta;) [i * sin (& phi;) + cos (& phi;)]

= [cos (& phi;) + i sin (& phi;)] [1 - 2 sin ² (& delta; / 2)] + [cos (& phi;) + i sin (& phi; )] я грешу (& delta;)

= e ^{i & phi;} + e ^{i & phi;} [- 2 sin ² (& delta; / 2) + i sin (& delta;)]

Редактировать : Это было много бесполезного форматирования с моей стороны. Это на самом деле намного проще:

y ^{(a + b)} = y ^a & times; y ^b для любого y. Итак:

e ^{i (& phi; + & delta;)}

= e ^{i & phi;} e ^{i & delta;}

= e ^{i & phi;} [cos (& delta;) + i sin (& delta;)]

= e ^{i & phi;} [1 - 2 sin ² (& delta; / 2) + i sin (& delta;)]

Answer 2

Причина в числовой точности. Если вы внимательно изучите следующий код:

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 точно для N >= 25736 для float (то есть 32-битная точность). Возможно БПФ на 25 736 точек. Таким образом, чтобы избежать этой проблемы, 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 3

Одна форма тождества половинного угла для косинусов:

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

Не уверен, отвечает ли это на ваш вопрос.

Answer 4

Что сбивает с толку, так это то, что NR использует математическую / физическую версию БПФ, которая вращает множители в противоположную сторону, что EE определяют БПФ.Таким образом, NR forward является обратным EE и наоборот.Обратите внимание, что в NR форвард имеет положительный показатель вместо отрицательного показателя EE.Метод EE преобразует время в частоту, где математическая и физическая версия играет с угловым моментом.

Answer 5

Хорошо, вот тригентная идентичность. Причина, по которой это не идентичность с половиной cos (theta) trig, заключается в зависимости от эйлеровых и мнимых чисел. Математика все еще за мной ...

_{(источник: librow.com )}

Answer 6

Они используют триггерные идентификаторы, чтобы минимизировать количество циклических функций, которые им нужно вычислить. Многие быстрые реализации просто ищут эти функции. Если вы действительно хотите знать, вам нужно просто проработать детали, развернув цикл выше и сделав соответствующие замены переменных ... утомительно да.

Кстати, реализация NR, как известно, очень медленная.

Пол

Answer 7

Я не знаю, хорошо БПФ я сделал один, но это было давно.

Итак, я посмотрел триггеры на http://www.sosmath.com/trig/Trig5/trig5/trig5.html

и из первой идентичности "Product-to-Sum" для sin (u) * sin (v) мы имеем

sin ^ 2 (тета / 2) = (1/2) [соз (ноль) - соз (тета)] = 0,5 - 0,5 соз (тета)

Помогает ли это?