Быстрое преобразование Фурье
Листинг программы на C++ (N - четное)
fft.h
/*
Fast Fourier Transformation
====================================================
Coded by Miroslav Voinarovsky, 2002
This source is freeware.
*/
#ifndef FFT_H_
#define FFT_H_
struct Complex;
struct ShortComplex;
/*
Fast Fourier Transformation
x: x - array of items
N: N - number of items in array. Must be odd
complement: false - normal (direct) transformation, true - reverse transformation
*/
void fft2(ShortComplex *x, int N, bool complement);
struct ShortComplex
{
double re, im;
inline void operator=(const Complex &y);
};
struct Complex
{
long double re, im;
inline void operator= (const Complex &y);
inline void operator= (const ShortComplex &y);
};
inline void ShortComplex::operator=(const Complex &y) { re = (double)y.re; im = (double)y.im; }
inline void Complex::operator= (const Complex &y) { re = y.re; im = y.im; }
inline void Complex::operator= (const ShortComplex &y) { re = y.re; im = y.im; }
#endif
fft.cpp
/*
Fast Fourier Transformation for even N
====================================================
Coded by Miroslav Voinarovsky, 2004
This source is freeware.
*/
#include "fft.h"
#include <math.h>
// This array contains values from 0 to 255 with reverse bit order
static unsigned char reverse256[]= {
0x00, 0x80, 0x40, 0xC0, 0x20, 0xA0, 0x60, 0xE0,
0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,
0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8,
0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,
0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,
0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,
0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,
0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,
0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,
0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,
0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,
0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,
0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,
0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,
0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,
0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,
0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,
0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,
0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,
0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,
0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,
0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,
0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED,
0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,
0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,
0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,
0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,
0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,
0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,
0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,
0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,
0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF,
};
//This is minimized version of type 'complex'. All operations is inline
static long double temp;
inline void operator+=(ShortComplex &x, const Complex &y) { x.re += (double)y.re; x.im += (double)y.im; }
inline void operator-=(ShortComplex &x, const Complex &y) { x.re -= (double)y.re; x.im -= (double)y.im; }
inline void operator*=(Complex &x, const Complex &y) { temp = x.re; x.re = temp * y.re - x.im * y.im; x.im = temp * y.im + x.im * y.re; }
inline void operator*=(Complex &x, const ShortComplex &y) { temp = x.re; x.re = temp * y.re - x.im * y.im; x.im = temp * y.im + x.im * y.re; }
inline void operator/=(ShortComplex &x, double div) { x.re /= div; x.im /= div; }
//This is array exp(-2*pi*j/2^n) for n= 1,...,32
//exp(-2*pi*j/2^n) = Complex( cos(2*pi/2^n), -sin(2*pi/2^n) )
static Complex W2n[32]={
{-1.00000000000000000000000000000000, 0.00000000000000000000000000000000}, // W2 calculator (copy/paste) : po, ps
{ 0.00000000000000000000000000000000, -1.00000000000000000000000000000000}, // W4: p/2=o, p/2=s
{ 0.70710678118654752440084436210485, -0.70710678118654752440084436210485}, // W8: p/4=o, p/4=s
{ 0.92387953251128675612818318939679, -0.38268343236508977172845998403040}, // p/8=o, p/8=s
{ 0.98078528040323044912618223613424, -0.19509032201612826784828486847702}, // p/16=
{ 0.99518472667219688624483695310948, -9.80171403295606019941955638886e-2}, // p/32=
{ 0.99879545620517239271477160475910, -4.90676743274180142549549769426e-2}, // p/64=
{ 0.99969881869620422011576564966617, -2.45412285229122880317345294592e-2}, // p/128=
{ 0.99992470183914454092164649119638, -1.22715382857199260794082619510e-2}, // p/256=
{ 0.99998117528260114265699043772857, -6.13588464915447535964023459037e-3}, // p/(2y9)=
{ 0.99999529380957617151158012570012, -3.06795676296597627014536549091e-3}, // p/(2y10)=
{ 0.99999882345170190992902571017153, -1.53398018628476561230369715026e-3}, // p/(2y11)=
{ 0.99999970586288221916022821773877, -7.66990318742704526938568357948e-4}, // p/(2y12)=
{ 0.99999992646571785114473148070739, -3.83495187571395589072461681181e-4}, // p/(2y13)=
{ 0.99999998161642929380834691540291, -1.91747597310703307439909561989e-4}, // p/(2y14)=
{ 0.99999999540410731289097193313961, -9.58737990959773458705172109764e-5}, // p/(2y15)=
{ 0.99999999885102682756267330779455, -4.79368996030668845490039904946e-5}, // p/(2y16)=
{ 0.99999999971275670684941397221864, -2.39684498084182187291865771650e-5}, // p/(2y17)=
{ 0.99999999992818917670977509588385, -1.19842249050697064215215615969e-5}, // p/(2y18)=
{ 0.99999999998204729417728262414778, -5.99211245264242784287971180889e-6}, // p/(2y19)=
{ 0.99999999999551182354431058417300, -2.99605622633466075045481280835e-6}, // p/(2y20)=
{ 0.99999999999887795588607701655175, -1.49802811316901122885427884615e-6}, // p/(2y21)=
{ 0.99999999999971948897151921479472, -7.49014056584715721130498566730e-7}, // p/(2y22)=
{ 0.99999999999992987224287980123973, -3.74507028292384123903169179084e-7}, // p/(2y23)=
{ 0.99999999999998246806071995015625, -1.87253514146195344868824576593e-7}, // p/(2y24)=
{ 0.99999999999999561701517998752946, -9.36267570730980827990672866808e-8}, // p/(2y25)=
{ 0.99999999999999890425379499688176, -4.68133785365490926951155181385e-8}, // p/(2y26)=
{ 0.99999999999999972606344874922040, -2.34066892682745527595054934190e-8}, // p/(2y27)=
{ 0.99999999999999993151586218730510, -1.17033446341372771812462135032e-8}, // p/(2y28)=
{ 0.99999999999999998287896554682627, -5.85167231706863869080979010083e-9}, // p/(2y29)=
{ 0.99999999999999999571974138670657, -2.92583615853431935792823046906e-9}, // p/(2y30)=
{ 0.99999999999999999892993534667664, -1.46291807926715968052953216186e-9}, // p/(2y31)=
};
#define M_2PI (6.283185307179586476925286766559)
inline void complex_mul(ShortComplex *z, const ShortComplex *z1, const Complex *z2)
{
z->re = (double)(z1->re * z2->re - z1->im * z2->im);
z->im = (double)(z1->re * z2->im + z1->im * z2->re);
}
static ShortComplex *createWstore(unsigned int L, bool complement)
{
unsigned int N, Skew, Skew2;
ShortComplex *Wstore, *Warray, *WstoreEnd;
Complex WN, *pWN;
Skew2 = L >> 1;
Wstore = new ShortComplex[Skew2];
WstoreEnd = Wstore + Skew2;
Wstore[0].re = 1.0;
Wstore[0].im = 0.0;
for(N = 4, pWN = W2n + 1, Skew = Skew2 >> 1; N <= L; N += N, pWN++, Skew2 = Skew, Skew >>= 1)
{
//WN = W(1, N) = exp(-2*pi*j/N)
WN= *pWN;
if (complement)
WN.im = -WN.im;
for(Warray = Wstore; Warray < WstoreEnd; Warray += Skew2)
complex_mul(Warray + Skew, Warray, &WN);
}
return Wstore;
}
bool fft_step(ShortComplex *x, unsigned int T, unsigned int M, const ShortComplex *Wstore)
{
unsigned int L, I, J, MI, MJ, ML, N, Nd2, k, m, Skew, mpNd2;
unsigned char *Ic = (unsigned char*) &I;
unsigned char *Jc = (unsigned char*) &J;
ShortComplex S;
const ShortComplex *Warray;
Complex Temp;
L = 1 << T;
ML = M * L;
//first interchanging
for(I = 1, MI = M; I < L - 1; I++, MI += M)
{
Jc[0] = reverse256[Ic[3]];
Jc[1] = reverse256[Ic[2]];
Jc[2] = reverse256[Ic[1]];
Jc[3] = reverse256[Ic[0]];
J >>= (32 - T);
if (I < J)
{
MJ = M * J;
S = x[MI];
x[MI] = x[MJ];
x[MJ] = S;
}
}
//main loop
for(Nd2 = M, N = M + M, Skew = L >> 1; N <= ML; Nd2 = N, N += N, Skew >>= 1)
{
for(Warray = Wstore, k = 0; k < Nd2; k += M, Warray += Skew)
{
for(m = k; m < ML; m += N)
{
Temp = *Warray;
mpNd2= m + Nd2;
Temp *= x[mpNd2];
x[mpNd2] = x[m];
x[mpNd2] -= Temp;
x[m] += Temp;
}
}
}
return true;
}
/*
x: x - array of items
N: N - number of items in array. Must be odd
complement: false - normal (direct) transformation, true - reverse transformation
*/
void fft2(ShortComplex *x, int N, bool complement)
{
int r, sL, m, rpsL, mprM, widx, step, step2, h, L, T, M;
//N is odd: decrease sequence
if (N & 1)
{
N--;
x[0].im = x[0].re = 0.0;
}
//find L, M and T
for(L = 1, T = 0; L < N && N % L == 0; L += L, T++)
{
}
//return to nearest good L, T
if (L != N)
{
L >>= 1;
T--;
}
//find M
M = N / L;
//find rotation multipliers
ShortComplex *Wstore= createWstore(L, complement);
//make usual FFT
for (h = 0; h < M; h++)
fft_step(x + h, T, M, Wstore);
//remove multipliers
delete [] Wstore;
if (M != 1)
{
//allocate temporary array
ShortComplex *X = new ShortComplex[N];
Complex one;
long double arg;
//rotation multiplier array allocation
ShortComplex *mult = new ShortComplex[N];
ShortComplex *multEnd = mult + N;
ShortComplex *multPtr;
mult[0].re= 1.0;
mult[0].im= 0.0;
//step2 is
for(step2 = 1; step2 < N; step2 += step2)
{
}
//rotation multiplier array evaluation
for(step = step2 >> 1; step > 0; step2 = step, step >>= 1)
{
arg= (complement ? M_2PI : -M_2PI) * step / N;
one.re= cosl(arg);
one.im= sinl(arg);
for(multPtr = mult + step; multPtr < multEnd; multPtr += step2)
complex_mul(multPtr, multPtr - step, &one);
}
ShortComplex *pX;
for(pX = X, rpsL = sL = 0; sL < N; sL += L)
{
for(r = 0; r < L; r++, rpsL++, pX++)
{
mprM = r*M;
*pX = x[mprM++];
widx = rpsL;
for(m = 1; m < M; m++, mprM++, widx += rpsL)
{
if (widx >= N)
widx -= N;
one = mult[widx];
one *= x[mprM];
*pX += one;
}
}
}
delete mult;
//copy result from temporary array and free it
for (h = 0; h < N; h++)
x[h] = X[h];
delete [] X;
}
//adjust values for complement transformation
if (complement)
{
for( h= 0; h < N; h++)
x[h] /= N;
}
}
Знаете ли Вы, что cогласно релятивистской мифологии "гравитационное линзирование - это физическое явление, связанное с отклонением лучей света в поле тяжести. Гравитационные линзы обясняют образование кратных изображений одного и того же астрономического объекта (квазаров, галактик), когда на луч зрения от источника к наблюдателю попадает другая галактика или скопление галактик (собственно линза). В некоторых изображениях происходит усиление яркости оригинального источника." (Релятивисты приводят примеры искажения изображений галактик в качестве подтверждения ОТО - воздействия гравитации на свет)
При этом они забывают, что поле действия эффекта ОТО - это малые углы вблизи поверхности звезд, где на самом деле этот эффект не наблюдается (затменные двойные). Разница в шкалах явлений реального искажения изображений галактик и мифического отклонения вблизи звезд - 10
11 раз. Приведу аналогию. Можно говорить о воздействии поверхностного натяжения на форму капель, но нельзя серьезно говорить о силе поверхностного натяжения, как о причине океанских приливов.
Эфирная физика находит ответ на наблюдаемое явление искажения изображений галактик. Это результат нагрева эфира вблизи галактик, изменения его плотности и, следовательно, изменения скорости света на галактических расстояниях вследствие преломления света в эфире различной плотности. Подтверждением термической природы искажения изображений галактик является прямая связь этого искажения с радиоизлучением пространства, то есть эфира в этом месте, смещение спектра CMB (космическое микроволновое излучение) в данном направлении в высокочастотную область. Подробнее читайте в
FAQ по эфирной физике .
НОВОСТИ ФОРУМА Рыцари теории эфира 10.11.2021 - 12:37: ПЕРСОНАЛИИ - Personalias -> WHO IS WHO - КТО ЕСТЬ КТО - Карим_Хайдаров.СОВЕСТЬ - Conscience -> РАСЧЕЛОВЕЧИВАНИЕ ЧЕЛОВЕКА. КОМУ ЭТО НАДО? - Карим_Хайдаров.ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от д.м.н. Александра Алексеевича Редько - Карим_Хайдаров.ЭКОЛОГИЯ - Ecology -> Биологическая безопасность населения - Карим_Хайдаров.ВОЙНА, ПОЛИТИКА И НАУКА - War, Politics and Science -> Проблема государственного терроризма - Карим_Хайдаров.ВОЙНА, ПОЛИТИКА И НАУКА - War, Politics and Science -> ПРАВОСУДИЯ.НЕТ - Карим_Хайдаров.ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от Вадима Глогера, США - Карим_Хайдаров.НОВЫЕ ТЕХНОЛОГИИ - New Technologies -> Волновая генетика Петра Гаряева, 5G-контроль и управление - Карим_Хайдаров.ЭКОЛОГИЯ - Ecology -> ЭКОЛОГИЯ ДЛЯ ВСЕХ - Карим_Хайдаров.ЭКОЛОГИЯ - Ecology -> ПРОБЛЕМЫ МЕДИЦИНЫ - Карим_Хайдаров.ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от Екатерины Коваленко - Карим_Хайдаров.ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от Вильгельма Варкентина - Карим_Хайдаров.
