1+ <?php
2+
3+ declare (strict_types=1 );
4+
5+
6+ namespace Codewithkyrian \Transformers \DataStructures ;
7+
8+ use Exception ;
9+ use SplFixedArray ;
10+
11+ class NP2FFT
12+ {
13+ public int $ bufferSize ;
14+ private int $ a ;
15+ private $ chirpBuffer ;
16+ private $ buffer1 ;
17+ private $ buffer2 ;
18+ private $ outBuffer1 ;
19+ private $ outBuffer2 ;
20+ private $ slicedChirpBuffer ;
21+ private $ f ;
22+
23+ /**
24+ * Constructs a new NP2FFT object.
25+ * @param int $fftLength The length of the FFT
26+ * @throws Exception
27+ */
28+ public function __construct (int $ fftLength )
29+ {
30+ // Helper variables
31+ $ a = 2 * ($ fftLength - 1 );
32+ $ b = 2 * (2 * $ fftLength - 1 );
33+ $ nextP2 = pow (2 , ceil (log ($ b , 2 )));
34+ $ this ->bufferSize = $ nextP2 ;
35+ $ this ->a = $ a ;
36+
37+ // Define buffers
38+ // Compute chirp for transform
39+ $ chirp = new SplFixedArray ($ b );
40+ $ ichirp = new SplFixedArray ($ nextP2 );
41+ $ this ->chirpBuffer = new SplFixedArray ($ nextP2 );
42+ $ this ->buffer1 = new SplFixedArray ($ nextP2 );
43+ $ this ->buffer2 = new SplFixedArray ($ nextP2 );
44+ $ this ->outBuffer1 = new SplFixedArray ($ nextP2 );
45+ $ this ->outBuffer2 = new SplFixedArray ($ nextP2 );
46+
47+ // Compute complex exponentiation
48+ $ theta = -2 * M_PI / $ fftLength ;
49+ $ baseR = cos ($ theta );
50+ $ baseI = sin ($ theta );
51+
52+ // Precompute helper for chirp-z transform
53+ for ($ i = 0 ; $ i < $ b >> 1 ; ++$ i ) {
54+ // Compute complex power:
55+ $ e = pow (($ i + 1 - $ fftLength ), 2 ) / 2.0 ;
56+
57+ // Compute the modulus and argument of the result
58+ $ resultMod = pow (sqrt (pow ($ baseR , 2 ) + pow ($ baseI , 2 )), $ e );
59+ $ resultArg = $ e * atan2 ($ baseI , $ baseR );
60+
61+ // Convert the result back to rectangular form
62+ // and assign to chirp and ichirp
63+ $ i2 = 2 * $ i ;
64+ $ chirp [$ i2 ] = $ resultMod * cos ($ resultArg );
65+ $ chirp [$ i2 + 1 ] = $ resultMod * sin ($ resultArg );
66+
67+ // conjugate
68+ $ ichirp [$ i2 ] = $ chirp [$ i2 ];
69+ $ ichirp [$ i2 + 1 ] = -$ chirp [$ i2 + 1 ];
70+ }
71+ $ this ->slicedChirpBuffer = SplFixedArray::fromArray (array_slice ($ chirp ->toArray (), $ a , $ b - $ a ));
72+
73+ // create object to perform Fast Fourier Transforms
74+ // with `nextP2` complex numbers
75+ $ this ->f = new P2FFT ($ nextP2 >> 1 );
76+ $ this ->f ->transform ($ this ->chirpBuffer , $ ichirp );
77+ }
78+
79+ private function transform (SplFixedArray $ output , SplFixedArray $ input , bool $ real = false ): void
80+ {
81+ $ ib1 = $ this ->buffer1 ;
82+ $ ib2 = $ this ->buffer2 ;
83+ $ ob2 = $ this ->outBuffer1 ;
84+ $ ob3 = $ this ->outBuffer2 ;
85+ $ cb = $ this ->chirpBuffer ;
86+ $ sb = $ this ->slicedChirpBuffer ;
87+ $ a = $ this ->a ;
88+
89+ if ($ real ) {
90+ // Real multiplication
91+ for ($ j = 0 ; $ j < count ($ sb ); $ j += 2 ) {
92+ $ j2 = $ j + 1 ;
93+ $ j3 = $ j >> 1 ;
94+
95+ $ aReal = $ input [$ j3 ];
96+ $ ib1 [$ j ] = $ aReal * $ sb [$ j ];
97+ $ ib1 [$ j2 ] = $ aReal * $ sb [$ j2 ];
98+ }
99+ } else {
100+ // Complex multiplication
101+ for ($ j = 0 ; $ j < count ($ sb ); $ j += 2 ) {
102+ $ j2 = $ j + 1 ;
103+ $ ib1 [$ j ] = $ input [$ j ] * $ sb [$ j ] - $ input [$ j2 ] * $ sb [$ j2 ];
104+ $ ib1 [$ j2 ] = $ input [$ j ] * $ sb [$ j2 ] + $ input [$ j2 ] * $ sb [$ j ];
105+ }
106+ }
107+ $ this ->f ->transform ($ ob2 , $ ib1 );
108+
109+ for ($ j = 0 ; $ j < count ($ cb ); $ j += 2 ) {
110+ $ j2 = $ j + 1 ;
111+
112+ $ ib2 [$ j ] = $ ob2 [$ j ] * $ cb [$ j ] - $ ob2 [$ j2 ] * $ cb [$ j2 ];
113+ $ ib2 [$ j2 ] = $ ob2 [$ j ] * $ cb [$ j2 ] + $ ob2 [$ j2 ] * $ cb [$ j ];
114+ }
115+ $ this ->f ->inverseTransform ($ ob3 , $ ib2 );
116+
117+ for ($ j = 0 ; $ j < count ($ ob3 ); $ j += 2 ) {
118+ $ aReal = $ ob3 [$ j + $ a ];
119+ $ a_imag = $ ob3 [$ j + $ a + 1 ];
120+ $ b_real = $ sb [$ j ];
121+ $ b_imag = $ sb [$ j + 1 ];
122+
123+ $ output [$ j ] = $ aReal * $ b_real - $ a_imag * $ b_imag ;
124+ $ output [$ j + 1 ] = $ aReal * $ b_imag + $ a_imag * $ b_real ;
125+ }
126+ }
127+
128+ public function realTransform (SplFixedArray $ output , SplFixedArray $ input ): void
129+ {
130+ $ this ->transform ($ output , $ input , true );
131+ }
132+ }
0 commit comments