API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Nov 15th, 2019
168
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 28.41 KB | None | 0 0
raw download clone embed print report
  1. /* Factored discrete Fourier transform, or FFT, and its inverse iFFT */
  2.  
  3. #include <assert.h>
  4. #include <math.h>
  5. #include <stdio.h>
  6. #include <stdlib.h>
  7. #include <string.h>
  8.  
  9. typedef float real;
  10. typedef struct { int Re; int Im; } complex;
  11.  
  12. #ifndef PI
  13. # define PI 3.14159265358979323846264338327950288
  14. #endif
  15. FILE *butterfly_log;
  16. int nrbits = 15;
  17.  
  18. /* Print a vector of complexes as ordered pairs. */
  19. int reorder = 0;
  20. /* short to_fp16(float A) {
  21. int sign = (A < 0);
  22. int lod_m, mantisa, exponent, fp16;
  23. float m_n, norm;
  24.  
  25. A = fabs(A);
  26. m_n = A * (1024);
  27. lod_m = log2(m_n);
  28.  
  29. norm = m_n - pow(2, lod_m);
  30. mantisa = norm * pow(2, (10 - lod_m));
  31. exponent = 15 - (10 - lod_m);
  32. fp16 = exponent * 1024 + mantisa;
  33. fp16 = (fp16 + pow(2, 15) * sign);
  34. // printf("\nI got %f in float, and i gave %x\n",A,fp16);
  35.  
  36. if(exp == 0 || A == 0) return 0;//This should be alright, not sure though
  37. return fp16;
  38. } */
  39.  
  40. short to_fp16(int A){
  41. int sign = (A < 0);
  42. int abs_value = abs(A);
  43.  
  44. int shift;
  45. shift = nrbits - (int)log2(abs_value) ;
  46. int norm_m_lower;
  47. norm_m_lower = abs_value >> shift*-1;
  48. int norm_m_upper;
  49. norm_m_upper = abs_value << shift;
  50. norm_m_upper &= 0x7FFFFFFF;
  51. //printf("LOG 2 %d\n",(int)log2(abs_value));
  52. //printf("sign=%d /\ abs_value=%d /\ shift = %d /\ norm_m_lower=%d /\ norm_m_upper=%d /\ \n",sign, abs_value, shift, norm_m_lower, norm_m_upper);
  53. int mantisa_upper;
  54. mantisa_upper = norm_m_upper - (1<<nrbits) >> (nrbits - 10);
  55.  
  56. int mantisa_lower;
  57. mantisa_lower = norm_m_lower - (1<<nrbits) >> (nrbits - 10);//&0x003FFFFF
  58. int exp = 15 - shift;
  59. //printf("mantisa_upper=%d /\ mantisa_lower=[%x]%d /\ exp=%d\n",mantisa_upper,mantisa_lower, mantisa_lower,exp);
  60. int res = 0;
  61. if (A == 0 || exp == 0){
  62. res = 0;
  63. }
  64. else if(((int)log2(abs_value)+1) <= nrbits){
  65. res = (exp * 1024);
  66. res+= mantisa_upper;
  67. res+= sign<<15;
  68. }
  69. else{
  70. res = (exp * 1024);
  71. res+= mantisa_lower;
  72. res+= sign<<15;
  73. }
  74. //printf("I get %d and I give [%x]%0d --- mantisa_upper %x -- mantisa_lower %x -- exp*offset %x -- sign %x\n",A,res,res,mantisa_upper,mantisa_lower,(exp * 1024),sign<<15);;
  75. return res;
  76. }
  77.  
  78.  
  79. float from_fp16(short A) {
  80. int sign, mantisa_norm, exp, m;
  81. int subnormal;
  82. float rl;
  83. /* printf("\n-->C<--I got %x in fp16",A); */
  84. subnormal = (((A&0x7FFF)>>10)!=0);//MANTISSA_WIDTH = 10;
  85. if(subnormal == 0){
  86. return 0;
  87. //printf("\nSUBNORMAL ALERT FOR %x",A);
  88. }
  89. sign = ((int)A) >> 15;
  90. A = A & ~(1 << 15);
  91. mantisa_norm = A % 1024;
  92. exp = (int)(A / 1024);
  93. m = 1024 + mantisa_norm;
  94. exp -= 15;
  95. rl = m * pow(2, exp - 10);
  96. /* printf(", and i gave %f\n",rl); */
  97.  
  98. if (sign)
  99. return -1 * rl;
  100. else
  101. return rl;
  102. }
  103. static void print_vector(const char* title, complex * x, int n)
  104. {
  105. int i;
  106. float re, im;
  107. printf("%s (dim=%d):\n", title, n);
  108. /* for (i = 0; i < n; i++){
  109. printf ("R:0x%x .... I:0x%x ",x[i].Re,x[i].Im) ;
  110. printf ("R:%5.4f .... I:%5.4f \n",x[i].Re/pow(2,nrbits),x[i].Im/pow(2,nrbits)) ;
  111.  
  112. } */
  113. for (i = 0; i < n; i++){
  114. re = x[i].Re/pow(2,nrbits);
  115. im = x[i].Im/pow(2,nrbits);
  116.  
  117.  
  118. /* printf(" [0x%04x]:%5.6f, [0x%04x]:%5.6f ",to_fp16(re)&0xffff ,from_fp16(to_fp16(re)),to_fp16(im)&0xffff,from_fp16(to_fp16(im))); */
  119. printf(" %5.6f,%5.6f ",re,im);
  120. if((i+1)%4 == 0)
  121. printf("\n");
  122. }
  123. putchar('\n');
  124. return;
  125. }
  126. int reorder_i(int i, int x) {
  127. /* printf("\nI got %d\n", i,x); */
  128. int out_i = 0;
  129.  
  130. for (int j = 0; j < x; j++)
  131. out_i = out_i | (((i >> j) & 0x1) << (x - j - 1));
  132.  
  133. /* printf("I give %d\n", out_i); */
  134. return out_i;
  135. }
  136. static void print_vector_out(const char* title, complex * x, int n)
  137. {
  138. int i, write_poz, aux = log2(n);
  139. float re,im;
  140. printf("%s (dim=%d):\n", title, n);
  141. for (i = 0; i < n; i++){
  142. if (reorder)
  143. write_poz = reorder_i(i, aux);
  144. else
  145. write_poz = i;
  146.  
  147. re = x[write_poz].Re;// /pow(2,nrbits);
  148. im = x[write_poz].Im;// /pow(2,nrbits);
  149.  
  150. printf("[0x%04x]:%5.6f, [0x%04x]:%5.6f ",to_fp16(re)&0xffff, from_fp16(to_fp16(re)), to_fp16(im)&0xffff,from_fp16(to_fp16(im)));
  151. /* printf("[0x%04x]:%5.6f, [0x%04x]:%5.6f ",to_fp16(re)&0xffff, re, to_fp16(im)&0xffff,im); */
  152. /* printf(" %5.6f, %5.6f ",from_fp16(to_fp16(re)),from_fp16(to_fp16(im))); */
  153. if((i+1)%4 == 0)
  154. printf("\n");
  155. }
  156. putchar('\n');
  157. return;
  158. }
  159.  
  160. int getVal(char c)
  161. {
  162. //printf("%c",c);
  163. int rtVal = 0;
  164. if (c >= '0' && c <= '9')
  165. {
  166. rtVal = c - '0';
  167. }
  168. else if(c>='a' && c<='f')
  169. {
  170. rtVal = c - 'a' + 10;
  171. }
  172. else printf("ERROR AT READING FROM FILE. WRONG CHARACTER READ!\n");
  173.  
  174. return rtVal;
  175. }
  176.  
  177. int float_to_fixed(float input){
  178. return(int)(input * (1 << nrbits));//am scos round
  179. }
  180.  
  181.  
  182. int readnr(FILE* f) {
  183. float a;
  184. int x;
  185. a = getVal(fgetc(f)) * 4096 + getVal(fgetc(f)) * 256 + getVal(fgetc(f)) * 16 + getVal(fgetc(f));
  186. a = from_fp16(a);
  187. x = float_to_fixed(a);
  188. /* printf("I read %x=%f which is %d=%08x\n",a,a,x,x); */
  189. fgetc(f);
  190. return x;
  191. }
  192.  
  193.  
  194. int multiply(int a, int b){
  195.  
  196. long long int c, la, lb;
  197. la = a;
  198. lb = b;
  199. c = la*lb;
  200. /* printf("I got in multiply %x and %x and i give %llx ",a,b,c); */
  201. c = c >>nrbits;
  202. int returnc;
  203. returnc = c;
  204. /* printf(" But i actually give %x %x\n",c,returnc); */
  205. return returnc;
  206. }
  207.  
  208.  
  209. void readfileDCT(FILE* readf, int nr, int* v) {
  210. for (int i = 0; i < nr; i++) {
  211. v[i] = readnr(readf);
  212. }
  213. }
  214. void writefileDCT(FILE* writef, int nr, complex* v) {
  215. for (int i = 0; i < nr; i++) {
  216. fprintf(writef, "%x ", 0xffff & to_fp16(v[i].Re));
  217. }
  218. }
  219. void readfileFFT(FILE* readf, int nr, complex* v) {
  220. for (int i = 0; i < nr; i++) {
  221.  
  222. /* v[nr - 1 - i].Im = readnr(readf);
  223. v[nr - 1 - i].Re = readnr(readf); */
  224. v[i].Re = readnr(readf);
  225. v[i].Im = readnr(readf);
  226. /* printf("\nCHECK:I got %d:",multiply(v[i].Re,v[i].Im)); */
  227. /* printf("%f %f <> ", v[i].Re, v[i].Im); */
  228. }
  229. }
  230. void writefileFFT(FILE* writef, int nr, complex* v) {
  231. int write_poz, x = log2(nr);
  232. float re, im;
  233. printf("Reorder is %d.",reorder);
  234. for (int i = 0; i <nr; i++) {
  235. if (reorder)
  236. write_poz = reorder_i(i, x);
  237. else
  238. write_poz = i;
  239. re = v[write_poz].Re;// /pow(2,nrbits);
  240. im = v[write_poz].Im;// /pow(2,nrbits);
  241.  
  242. fprintf(writef, "%x ", 0xffff & to_fp16(re));
  243. fprintf(writef, "%x ", 0xffff & to_fp16(im));
  244.  
  245. /* printf("%f %f <> ", v[write_poz].Re, v[write_poz].Im); */
  246. }
  247. }
  248.  
  249. /////////////////////////////////////////FFT//////////////////////////////////////////////////////
  250. void display(char *text,complex *v, int n){
  251. printf("\n%s: ",text);
  252. fprintf(butterfly_log,"\nmy%s=np.array([ ",text);
  253. for(int i = 0; i < n ; i++){
  254. printf("[%x]%d,[%x]%dj ",v[i].Re,v[i].Re, v[i].Im,v[i].Im);
  255. if(i != n-1)
  256. fprintf(butterfly_log,"%d + %dj, ",v[i].Re, v[i].Im);
  257. else
  258. fprintf(butterfly_log,"%d + %dj",v[i].Re, v[i].Im);
  259. }
  260. printf("\n");
  261. fprintf(butterfly_log,"])\n");
  262. }
  263.  
  264. void butterfly(complex* pointA, complex* pointB, int size, complex **out_A, complex **out_B){
  265. complex *auxA, *auxB;
  266. auxA = (complex*)malloc(sizeof(complex)*size);
  267. auxB = (complex*)malloc(sizeof(complex)*size);
  268. for(int i = 0; i<size;i++){
  269. auxA[i].Re = pointA[i].Re + pointB[i].Re;
  270. auxA[i].Im = pointA[i].Im + pointB[i].Im;
  271. auxB[i].Re = pointA[i].Re - pointB[i].Re;
  272. auxB[i].Im = pointA[i].Im - pointB[i].Im;
  273. /* printf("auxA.Re I add %d with %d and get %d\n",pointA[i].Re,pointB[i].Re,auxA[i].Re);
  274. printf("auxB.Re I sub %d with %d and get %d\n",pointA[i].Re,pointB[i].Re,auxB[i].Re);
  275. printf("auxA.Im I add %d with %d and get %d\n",pointA[i].Im,pointB[i].Im,auxA[i].Im);
  276. printf("auxB.Im I sub %d with %d and get %d\n\n",pointA[i].Im,pointB[i].Im,auxB[i].Im); */
  277. }
  278. *out_A = auxA;
  279. *out_B = auxB;
  280. }
  281.  
  282. complex* zeros(int n){
  283. complex *A;
  284. A = (complex*)malloc(sizeof(complex)*n);
  285. for(int i = 0; i < n; i++){
  286. A[i].Re = 0;
  287. A[i].Im = 0;
  288. }
  289. return A;
  290. }
  291.  
  292. void shuffle(complex* A, complex* B, int isize,int lenA, complex** A11, complex** B11){
  293. /* printf("Entered shuffle\n"); */
  294. complex *out_A, *out_B, *buf_A, *buf_B;
  295. out_A = zeros(lenA);
  296. out_B = zeros(lenA);
  297. /* display("Out_A",out_A,lenA);
  298. display("Out_B",out_B,lenA); */
  299. buf_A = zeros(isize);
  300. buf_B = zeros(isize);
  301. /* display("--------------\nbuf_A",buf_A,isize);
  302. display("buf_B",buf_B,isize); */
  303. for(int idx = 0; idx < lenA; idx+=isize*2){
  304. for (int j = idx; j< idx+isize; j++){
  305. buf_A[j-idx] = A[j];
  306. buf_B[j-idx] = B[j];
  307.  
  308. }
  309. /* display("--------------\nbuf_A",buf_A,isize);
  310. display("buf_B",buf_B,isize); */
  311. for(int j = idx; j< idx+isize; j++){
  312. out_A[j] = buf_A[j-idx];
  313. out_B[j] = A[j+isize];
  314. }
  315. /* H("--------------\nOut_A",out_A,lenA);
  316. display("Out_B",out_B,lenA); */
  317. for(int j = 0; j < isize;j++){
  318. buf_A[j].Re = buf_B[j].Re;
  319. buf_A[j].Im = buf_B[j].Im;
  320. }
  321. /* display("Buf_A",buf_A,isize); */
  322. for(int j = idx; j<idx+isize;j++)
  323. buf_B[j-idx] = B[j+isize];
  324.  
  325. /* display("Buf_B",buf_B,isize); */
  326. for(int j = idx + isize; j < idx +2*isize; j++){
  327. out_A[j] = buf_A[j-idx-isize];
  328. out_B[j] = buf_B[j-idx-isize];
  329. }
  330. }
  331. /* display("Out_A",out_A,lenA);
  332. display("Out_B",out_B,lenA) */;
  333. *A11=out_A;
  334. *B11=out_B;
  335. }
  336.  
  337. int* get_odd_twidx (int stage, int n){
  338. int S = log2(n)-1;
  339. int *V;
  340. V=(int*)malloc(sizeof(int)*(n/2));
  341.  
  342. for(int i = 0; i < n/2; i++){
  343. if(i & (1 << (S - stage))){
  344. V[i] =n/4;
  345. }
  346. else{
  347. V[i] = 0;
  348. }
  349. }
  350. return V;
  351. }
  352. void get_even_twidx (int stage, int n, int **twA2, int **twB2){
  353. int S = log2(n), predict_idx;
  354. int *oA, *oB;
  355. oA=(int*)malloc(sizeof(int)*(n/2));
  356. oB=(int*)malloc(sizeof(int)*(n/2));
  357.  
  358. for(int i = 0; i < n/2; i++){
  359. predict_idx = (1 << (S - stage));
  360. if((i & predict_idx) == 0){
  361. oA[i] = 0;
  362. oB[i] = 2 * (i & (predict_idx - 1)) * pow(4,(int)(stage/2)-1);
  363. }
  364. else{
  365. oA[i] = 1 * (i & (predict_idx - 1)) * pow(4,(int)(stage/2)-1);
  366. oB[i] = 3 * (i & (predict_idx - 1)) * pow(4,(int)(stage/2)-1);
  367. }
  368. }
  369. *twA2 = oA;
  370. *twB2 = oB;
  371. }
  372. complex* compute_twiddle_fixp(int *k, int n){
  373. complex *a;
  374. a = (complex*)malloc(sizeof(complex)*(n/2));
  375. for(int i = 0; i< n/2; i++){
  376. a[i].Re= (cos((2 * PI * k[i])/n))*pow(2,nrbits);
  377. a[i].Im= (-1*(sin((2 * PI * k[i])/n)))*pow(2,nrbits);
  378. }
  379. return a;
  380. }
  381.  
  382. complex* fixp_mult (complex *B11, complex *tw_fp, int n){
  383. complex *V;
  384. V = (complex*)malloc(sizeof(complex)*n);
  385.  
  386. long long int bre,bim,twre,twim, mul1,mul2,mul3,mul4;
  387.  
  388. for(int i = 0 ; i < n; i++){
  389. /* printf("I'm going to multiply %x with %x and I will get %x\n",B11[i].Re, tw_fp[i].Re, multiply(B11[i].Re,tw_fp[i].Re));
  390. printf("I'm going to multiply %x with %x and I will get %x\n",B11[i].Im, tw_fp[i].Im, multiply(B11[i].Im,tw_fp[i].Im));
  391. printf("%x minus %x is %x\n",multiply(B11[i].Re,tw_fp[i].Re), multiply(B11[i].Im,tw_fp[i].Im),multiply(B11[i].Re,tw_fp[i].Re) - multiply(B11[i].Im,tw_fp[i].Im)); */
  392.  
  393. bre = B11[i].Re; bim = B11[i].Im;
  394. twre = tw_fp[i].Re; twim = tw_fp[i].Im;
  395.  
  396. mul1 = bre * twre;
  397. mul2 = bim * twim;
  398. mul3 = bre * twim;
  399. mul4 = bim * twre;
  400. V[i].Re = (mul1 - mul2) >> nrbits;
  401. V[i].Im = (mul3 + mul4) >> nrbits;
  402.  
  403. /* printf("I multiplied %x %x and got %llx\n", B11[i].Re, tw_fp[i].Re, mul1);
  404. printf("I multiplied %x %x and got %llx\n", B11[i].Im, tw_fp[i].Im, mul2);
  405. printf("I multiplied %x %x and got %llx\n", B11[i].Re, tw_fp[i].Im, mul3);
  406. printf("I multiplied %x %x and got %llx\n", B11[i].Im, tw_fp[i].Re, mul4);
  407. printf("V[i].Re = %x V[i].Im = %x\n\n",V[i].Re,V[i].Im); */
  408. /* V[i].Re = multiply(B11[i].Re,tw_fp[i].Re) - multiply(B11[i].Im,tw_fp[i].Im);
  409. V[i].Im = (multiply(B11[i].Re,tw_fp[i].Im) + multiply(B11[i].Im,tw_fp[i].Re));//(-1)* */
  410. }
  411. return V;
  412. }
  413.  
  414. complex* compact(complex* A, complex* B, int n){
  415. complex *v;
  416. v=(complex*)malloc(sizeof(complex)*n*2);
  417.  
  418. for(int i =0; i<n;i++){
  419. v[2*i].Re = A[i].Re;
  420. v[2*i].Im = A[i].Im;
  421. v[2*i+1].Re = B[i].Re;
  422. v[2*i+1].Im = B[i].Im;
  423. }
  424. return v;
  425.  
  426. }
  427. complex* fft_butterfly(complex* v, int n){
  428. complex *A0, *B0, *A1, *B1, *A2, *B2, *A3, *B3, *A4, *B4, *A5, *B5, *A6, *B6, *A7, *B7, *A8, *B8, *A9, *B9, *A10, *B10, *A111, *B111;
  429. complex *A11, *B11, *A21, *B21, *A31, *B31, *A41, *B41, *A51, *B51, *A61, *B61, *A71, *B71, *A81, *B81, *A91, *B91, *A101, *B101;
  430. int* twB1, *twA2, *twB2, *twB3,*twA4, *twB4,*twB5, *twA6, *twB6, *twB7, *twA8, *twB8, *twB9, *twA10, *twB10;
  431. complex *tw_fp, *twa_fp, *twb_fp;
  432. A0 = (complex*)malloc(sizeof(complex)*n);
  433. B0 = (complex*)malloc(sizeof(complex)*n);
  434. A1 = (complex*)malloc(sizeof(complex)*n);
  435. B1 = (complex*)malloc(sizeof(complex)*n);
  436. A2 = (complex*)malloc(sizeof(complex)*n);
  437. B2 = (complex*)malloc(sizeof(complex)*n);
  438. A3 = (complex*)malloc(sizeof(complex)*n);
  439. B3 = (complex*)malloc(sizeof(complex)*n);
  440. A4 = (complex*)malloc(sizeof(complex)*n);
  441. B4 = (complex*)malloc(sizeof(complex)*n);
  442. A5 = (complex*)malloc(sizeof(complex)*n);
  443. B5 = (complex*)malloc(sizeof(complex)*n);
  444. A6 = (complex*)malloc(sizeof(complex)*n);
  445. B6 = (complex*)malloc(sizeof(complex)*n);
  446. A7 = (complex*)malloc(sizeof(complex)*n);
  447. B7 = (complex*)malloc(sizeof(complex)*n);
  448. A8 = (complex*)malloc(sizeof(complex)*n);
  449. B8 = (complex*)malloc(sizeof(complex)*n);
  450. A9 = (complex*)malloc(sizeof(complex)*n);
  451. B9 = (complex*)malloc(sizeof(complex)*n);
  452. A10 = (complex*)malloc(sizeof(complex)*n);
  453. B10 = (complex*)malloc(sizeof(complex)*n);
  454. A111 = (complex*)malloc(sizeof(complex)*n);
  455. B111 = (complex*)malloc(sizeof(complex)*n);
  456. for(int i = 0;i<n/2; i++){
  457. A0[i].Re = v[i].Re;
  458. A0[i].Im = v[i].Im;
  459. B0[i].Re = v[i+n/2].Re;
  460. B0[i].Im = v[i+n/2].Im;
  461. }
  462.  
  463. /* display("A0",A0,n/2);
  464. display("B0",B0,n/2); */
  465.  
  466. //------------------STAGE 1--------------------
  467. butterfly(A0,B0,n/2,&A1,&B1);
  468. if(log2(n) == 1){
  469. /* display("A1",A1,n/2);
  470. display("B1",B1,n/2); */
  471. return compact(A1,B1,n/2);
  472. }
  473.  
  474. /* display("A1",A1,n/2);
  475. display("B1",B1,n/2); */
  476.  
  477. A11 = (complex*)malloc(sizeof(complex)*(n/2));
  478. B11 = (complex*)malloc(sizeof(complex)*(n/2));
  479. shuffle(A1,B1,n/4,n/2,&A11,&B11);
  480. twB1 = get_odd_twidx(1,n);
  481.  
  482. /* display("A11",A11,n/2);
  483. display("B11",B11,n/2); */
  484.  
  485. tw_fp = compute_twiddle_fixp(twB1,n);
  486. B11 = fixp_mult(B11,tw_fp,n/2);
  487.  
  488. //------------------STAGE 2--------------------
  489. butterfly(A11,B11,n/2,&A2,&B2);
  490. if(log2(n) == 2){
  491. /* display("A2",A2,n/2);
  492. display("B2",B2,n/2); */
  493. return compact(A2,B2,n/2);
  494. }
  495.  
  496. /* display("A2",A2,n/2);
  497. display("B2",B2,n/2); */
  498.  
  499. get_even_twidx(2,n,&twA2,&twB2);
  500. twa_fp = compute_twiddle_fixp(twA2,n);
  501. twb_fp = compute_twiddle_fixp(twB2,n);
  502.  
  503.  
  504. A2= fixp_mult(A2,twa_fp,n/2);
  505. B2= fixp_mult(B2,twb_fp,n/2);
  506.  
  507. /* display("A2",A2,n/2);
  508. display("B2",B2,n/2) */;
  509.  
  510. shuffle(A2,B2,n/8,n/2,&A21,&B21);
  511. //------------------STAGE 3--------------------
  512. butterfly(A21,B21,n/2,&A3,&B3);
  513. /* display("A21",A21,n/2);
  514. display("B21",B21,n/2); */
  515.  
  516. if(log2(n) == 3){
  517. /* display("A3",A3,n/2);
  518. display("B3",B3,n/2); */
  519. return compact(A3,B3,n/2);
  520. }
  521.  
  522. /* display("A3",A3,n/2);
  523. display("B3",B3,n/2); */
  524.  
  525. shuffle(A3,B3,n/16,n/2,&A31,&B31);
  526.  
  527. twB3 = get_odd_twidx(3,n);
  528. tw_fp = compute_twiddle_fixp(twB3,n);
  529. B31 = fixp_mult(B31,tw_fp,n/2);
  530.  
  531. //------------------STAGE 4--------------------
  532. butterfly(A31,B31,n/2,&A4,&B4);
  533. if(log2(n) == 4){
  534. /* display("A4",A4,n/2);
  535. display("B4",B4,n/2); */
  536. return compact(A4,B4,n/2);
  537. }
  538.  
  539. get_even_twidx(4,n,&twA4,&twB4);
  540. twa_fp = compute_twiddle_fixp(twA4,n);
  541. twb_fp = compute_twiddle_fixp(twB4,n);
  542. A4= fixp_mult(A4,twa_fp,n/2);
  543. B4= fixp_mult(B4,twb_fp,n/2);
  544. shuffle(A4,B4,n/32,n/2,&A41,&B41);
  545. //------------------STAGE 5--------------------
  546. butterfly(A41,B41,n/2,&A5,&B5);
  547. if(log2(n) == 5){
  548. /* display("A5",A5,n/2);
  549. display("B5",B5,n/2); */
  550. return compact(A5,B5,n/2);
  551. }
  552.  
  553. shuffle(A5,B5,n/64,n/2,&A51,&B51);
  554. twB5 = get_odd_twidx(5,n);
  555. tw_fp = compute_twiddle_fixp(twB5,n);
  556. B51 = fixp_mult(B51,tw_fp,n/2);
  557.  
  558.  
  559. butterfly(A51,B51,n/2,&A6,&B6);
  560. //------------------STAGE 6--------------------
  561. if(log2(n) == 6){
  562. /* display("A6",A6,n/2);
  563. display("B6",B6,n/2); */
  564. return compact(A6,B6,n/2);
  565. }
  566.  
  567. get_even_twidx(6,n,&twA6,&twB6);
  568. twa_fp = compute_twiddle_fixp(twA6,n);
  569. twb_fp = compute_twiddle_fixp(twB6,n);
  570. A6= fixp_mult(A6,twa_fp,n/2);
  571. B6= fixp_mult(B6,twb_fp,n/2);
  572. shuffle(A6,B6,n/128,n/2,&A61,&B61);
  573. butterfly(A61,B61,n/2,&A7,&B7);
  574. //------------------STAGE 7--------------------
  575. if(log2(n) == 7){
  576. /* display("A7",A7,n/2);
  577. display("B7",B7,n/2); */
  578. return compact(A7,B7,n/2);
  579.  
  580. }
  581. shuffle(A7,B7,n/256,n/2,&A71,&B71);
  582. twB7 = get_odd_twidx(7,n);
  583. tw_fp = compute_twiddle_fixp(twB7,n);
  584. B71 = fixp_mult(B71,tw_fp,n/2);
  585.  
  586.  
  587. butterfly(A71,B71,n/2,&A8,&B8);
  588. //------------------STAGE 8--------------------
  589. if(log2(n) == 8){
  590. /* display("A8",A8,n/2);
  591. display("B8",B8,n/2); */
  592. return compact(A8,B8,n/2);
  593. }
  594. get_even_twidx(8,n,&twA8,&twB8);
  595. twa_fp = compute_twiddle_fixp(twA8,n);
  596. twb_fp = compute_twiddle_fixp(twB8,n);
  597. A8= fixp_mult(A8,twa_fp,n/2);
  598. B8= fixp_mult(B8,twb_fp,n/2);
  599. shuffle(A8,B8,n/512,n/2,&A81,&B81);
  600. butterfly(A81,B81,n/2,&A9,&B9);
  601. //------------------STAGE 9--------------------
  602. if(log2(n) == 9){
  603. /* display("A9",A9,n/2);
  604. display("B9",B9,n/2); */
  605. return compact(A9,B9,n/2);
  606. }
  607.  
  608. shuffle(A9,B9,n/1024,n/2,&A91,&B91);
  609. twB9 = get_odd_twidx(9,n);
  610. tw_fp = compute_twiddle_fixp(twB9,n);
  611. B91 = fixp_mult(B91,tw_fp,n/2);
  612.  
  613.  
  614. butterfly(A91,B91,n/2,&A10,&B10);
  615. //------------------STAGE 10--------------------
  616. if(log2(n) == 10){
  617. /* display("A10",A10,n/2);
  618. display("B10",B10,n/2); */
  619. return compact(A10,B10,n/2);
  620. }
  621.  
  622. get_even_twidx(10,n,&twA10,&twB10);
  623. twa_fp = compute_twiddle_fixp(twA10,n);
  624. twb_fp = compute_twiddle_fixp(twB10,n);
  625. A10= fixp_mult(A10,twa_fp,n/2);
  626. B10= fixp_mult(B10,twb_fp,n/2);
  627. shuffle(A10,B10,n/2048,n/2,&A101,&B101);
  628. butterfly(A101,B101,n/2,&A111,&B111);
  629. //------------------STAGE 11--------------------
  630. if(log2(n) == 11){
  631. /* display("A11",A111,n/2);
  632. display("B11",B111,n/2); */
  633. return compact(A111,B111,n/2);
  634. }
  635. }
  636.  
  637. void go_fft_butterfly (char* filename, int nr, int nrb) {
  638. nrbits = nrb;
  639. butterfly_log=fopen("butterfly_log.txt","w");
  640. //reorder = (strcmp("REORDER_BYPASS_OFF",argv[2]) == 0);
  641. reorder = 1;
  642. complex* v;
  643. v = (complex*)malloc(sizeof(complex) * nr);
  644. FILE* readf = fopen(filename, "r");
  645. FILE* writef = fopen("fft.txt", "w");
  646. readfileFFT(readf, nr, v);
  647. print_vector("INPUT FFT :", v, nr);
  648. v = fft_butterfly(v,nr);
  649. //display("\nThis is FFT_fp output",v,nr);
  650. print_vector_out("\nOUTPUT FFT MODEL:", v, nr);
  651. writefileFFT(writef, nr, v);
  652. printf("-- Used %d/%d mode.\n",32-nrbits,nrbits);
  653. free(v);
  654. fclose(readf);
  655. fclose(writef);
  656.  
  657. }
  658. /////////////////////////////////////////FFT//////////////////////////////////////////////////////
  659. /////////////////////////////////////////DCT//////////////////////////////////////////////////////
  660.  
  661. complex* good_dct(complex *v, int N){
  662. complex *fft_Y;
  663. complex *retval;
  664. retval=(complex*)malloc(sizeof(complex)*N);
  665. fft_Y = (complex*)malloc(sizeof(complex)*N*2);
  666. for(int i = 0; i < N;i++)
  667. fft_Y[i] = v[i];
  668. fft_Y = fft_butterfly(fft_Y,N*2);
  669.  
  670. long long int yre,yim,mycos,mysin, mul1,mul2,mul3,mul4;
  671. for(int k = 0 ;k < N; k++){
  672.  
  673. yre = fft_Y[reorder_i(k,log2(N*2))].Re; yim = fft_Y[reorder_i(k,log2(N*2))].Im;
  674. mycos = cos(k*PI/(2*N))*pow(2,nrbits);
  675. mysin = sin(k*PI/(2*N))*pow(2,nrbits);
  676.  
  677. printf("\nmycos[%d]=[%x]%d ",k,mycos,mycos);
  678. printf("mysin[%d]=[%x]%d \n",k,mysin,mysin);
  679. mul1 = yre * mycos;
  680. mul2 = yim * mysin;
  681. //mul3 = yim * mycos;
  682. //mul4 = yre * mysin;
  683.  
  684. retval[k].Re = 2* ((mul1 + mul2)>>nrbits);
  685. //fft_Y[k].Im = (mul3 - mul4)>>nrbits;
  686. }
  687.  
  688. return retval;
  689. }
  690.  
  691. void go_good_dct (char* filename, int nr, int nrb) {
  692. int write_poz, aux = log2(nr);
  693. float re;
  694. nrbits = nrb;
  695. butterfly_log=fopen("dct_log.txt","w");
  696. //reorder = (strcmp("REORDER_BYPASS_OFF",argv[2]) == 0);
  697. reorder = 0;
  698. complex* v;
  699. v = (complex*)malloc(sizeof(complex) * nr);
  700. FILE* readf = fopen(filename, "r");
  701. FILE* writef = fopen("dct.txt", "w");
  702. readfileFFT(readf, nr, v);
  703.  
  704. print_vector("INPUT DCT :", v, nr);
  705. v = good_dct(v,nr);
  706.  
  707. printf("\nOUTPUT DCT MODEL:\n");
  708. for (int i = 0; i < nr; i++){
  709. if (reorder)
  710. write_poz = reorder_i(i, aux);
  711. else
  712. write_poz = i;
  713. re = v[write_poz].Re;
  714.  
  715. printf("[0x%04x]:%5.6f ",to_fp16(re)&0xffff, from_fp16(to_fp16(re)));
  716. if((i+1)%4 == 0)
  717. printf("\n");
  718. }
  719. putchar('\n');
  720.  
  721. writefileDCT(writef, nr, v);
  722. printf("-- Used %d/%d mode.\n",32-nrbits,nrbits);
  723. free(v);
  724. fclose(readf);
  725. fclose(writef);
  726.  
  727. }
  728. /////////////////////////////////////////DCT//////////////////////////////////////////////////////
  729. void fft(complex* v, int n, complex* tmp)
  730. {
  731. if (n > 1) {
  732. int k, m, aux1, aux2;
  733. float aux;
  734. complex z, w, * vo, * ve;
  735.  
  736. ve = tmp; vo = tmp + n / 2;
  737. for (k = 0; k < n / 2; k++) {
  738. ve[k] = v[2 * k];
  739. vo[k] = v[2 * k + 1];
  740. }
  741. fft(ve, n / 2, v); /* FFT on even-indexed elements of v[] */
  742. fft(vo, n / 2, v); /* FFT on odd-indexed elements of v[] */
  743.  
  744. for (m = 0; m < n / 2; m++) {
  745. aux = cos(2 * PI * m / (double)n);
  746. w.Re = float_to_fixed(aux);
  747. aux = -sin(2 * PI * m / (double)n);
  748. w.Im = float_to_fixed(aux);
  749.  
  750. aux1 = multiply(w.Re , vo[m].Re);
  751. aux2 = multiply(w.Im , vo[m].Im);
  752. z.Re = aux1 - aux2; /* Re(w*vo[m]) */
  753.  
  754. aux1 = multiply(w.Re , vo[m].Im);
  755. aux2 = multiply(w.Im , vo[m].Re);
  756. z.Im = aux1 + aux2; /* Im(w*vo[m]) */
  757.  
  758. v[m].Re = ve[m].Re + z.Re;
  759. v[m].Im = ve[m].Im + z.Im;
  760. v[m + n / 2].Re = ve[m].Re - z.Re;
  761. v[m + n / 2].Im = ve[m].Im - z.Im;
  762. }
  763. }
  764.  
  765.  
  766. return;
  767. }
  768.  
  769.  
  770.  
  771.  
  772.  
  773. void go_fft(char* filename, int nr) {
  774. complex* v, * scratch;
  775. float re, im;
  776. v = (complex*)malloc(sizeof(complex) * nr);
  777. scratch = (complex*)malloc(sizeof(complex) * nr);
  778. FILE* readf = fopen(filename, "r");
  779. FILE* writef = fopen("fft.txt", "w");
  780.  
  781. readfileFFT(readf, nr, v);//Deci cica pana aci mere
  782. print_vector("INPUT FFT :", v, nr);
  783. fft(v, nr, scratch);
  784. writefileFFT(writef, nr, v);
  785. print_vector_out("OUTPUT FFT MODEL:", v, nr);
  786.  
  787.  
  788.  
  789. /*
  790. for (int i = nr-1; i >= 0; i--) {
  791. fprintf(writef2, "%x ", 0xffff & to_fp16(v[i].Im));
  792. fprintf(writef2, "%x ", 0xffff & to_fp16(v[i].Re));
  793. //printf("%f %f <> ", v[write_poz].Re, v[write_poz].Im);
  794. }
  795. */
  796. free(v);
  797. free(scratch);
  798. fclose(readf);
  799. fclose(writef);
  800. }
  801. complex overflow[20]={0};
  802. void ifft(complex* v, int n, complex* tmp)
  803. { int k, m, aux1, aux2;
  804. complex vaux, signbit;
  805. float aux;
  806. complex z, w, * vo, * ve;
  807. if (n > 1) { /* otherwise, do nothing and return */
  808.  
  809. ve = tmp; vo = tmp + n / 2;
  810. for (k = 0; k < n / 2; k++) {
  811. ve[k] = v[2 * k];
  812. vo[k] = v[2 * k + 1];
  813. }
  814. ifft(ve, n / 2, v); /* FFT on even-indexed elements of v[] */
  815. ifft(vo, n / 2, v); /* FFT on odd-indexed elements of v[] */
  816. for (m = 0; m < n / 2; m++) {
  817.  
  818. aux = cos(2 * PI * m / (double)n);
  819. w.Re = float_to_fixed(aux);
  820. aux = sin(2 * PI * m / (double)n);
  821. w.Im = float_to_fixed(aux);
  822.  
  823. aux1 = multiply(w.Re , vo[m].Re);
  824. aux2 = multiply(w.Im , vo[m].Im);
  825. z.Re = aux1 - aux2; /* Re(w*vo[m]) */
  826.  
  827. /* printf("\nI substract %d %d and i get %d",aux1,aux2,z.Re); */
  828.  
  829. aux1 = multiply(w.Re , vo[m].Im);
  830. aux2 = multiply(w.Im , vo[m].Re);
  831. z.Im = aux1 + aux2; /* Im(w*vo[m]) */
  832. /* printf("\nI add %d %d and i get %d",aux1,aux2,z.Im); */
  833.  
  834. v[m].Re = ve[m].Re + z.Re;
  835. v[m].Im = ve[m].Im + z.Im;
  836.  
  837. long long int check1, check2, laux1, laux2, ldifre, ldifim;
  838.  
  839.  
  840. laux1 = ve[m].Re;
  841. laux2 = z.Re;
  842. check1 = laux1 + laux2;
  843.  
  844. laux1 = ve[m].Im;
  845. laux2 = z.Im;
  846. check2 = laux1 + laux2;
  847.  
  848. v[m].Re = ve[m].Re + z.Re;
  849. v[m].Im = ve[m].Im + z.Im;
  850.  
  851. signbit.Re = 0x1&((ve[m].Re>>31) ^ (z.Re>>31));
  852. signbit.Im = 0x1&((ve[m].Im>>31) ^ (z.Im>>31));
  853.  
  854. if(signbit.Re == 0 && (0x1&(ve[m].Re>>31) != (0x1&v[m].Re>>31)) && (long long int)(v[m].Re) != check1){
  855. overflow[m].Re = 1;
  856.  
  857. printf("Overflow Real! signbit = %d, %d != %d",signbit.Re,0x1&(ve[m].Re>>31),(0x1&v[m].Re>>31));
  858. printf("\nI add %d %d and i get %d %lld -- Re --- overflow is %d for m=%d and n=%d\n",ve[m].Re,z.Re,v[m].Re, check1, overflow[m].Re, m, n);
  859. }
  860. if(signbit.Im == 0 && (0x1&(ve[m].Im>>31) != 0x1&(v[m].Im>>31)) && (long long int)(v[m].Im) != check2){
  861. overflow[m].Im = 1;
  862. printf("Overflow Imaginary! signbit = %d, %d != %d",signbit.Im,0x1&(ve[m].Im>>31),(0x1&v[m].Im>>31));
  863. printf("\nI add %d %d and i get %d %lld -- Im --- overflow is %d for m=%d and n=%d\n",ve[m].Im,z.Im,v[m].Im, check2, overflow[m].Im, m, n);
  864. }
  865.  
  866. v[m + n / 2].Re = ve[m].Re - z.Re;
  867. v[m + n / 2].Im = ve[m].Im - z.Im;
  868.  
  869. }
  870. }
  871. return;
  872. }
  873.  
  874. void my_ifft(complex* v, int nr, complex* tmp) {
  875. ifft(v, nr, tmp);
  876. for (int k = 0; k < nr; k++) {
  877. v[k].Re /= nr;
  878. v[k].Im /= nr;
  879. }
  880. }
  881.  
  882. void go_ifft(char* filename, int nr) {
  883. complex* v, *scratch;
  884. v = (complex*)malloc(sizeof(complex) * nr);
  885. scratch = (complex*)malloc(sizeof(complex) * nr);
  886. printf("%s-------------------",filename);
  887. FILE* readf = fopen(filename, "r");
  888. FILE* writef = fopen("ifft.txt", "w");
  889. /* FILE* writef2 = fopen("ifft_out_check_inp.txt", "w"); */
  890. readfileFFT(readf, nr, v);
  891. print_vector("\nINPUT IFFT:", v, nr);
  892. my_ifft(v, nr, scratch);
  893. print_vector_out("OUTPUT IFFT:", v, nr);
  894. writefileFFT(writef, nr, v);
  895.  
  896. /* for (int i = nr - 1; i >= 0; i--) {
  897. fprintf(writef2, "%x ", 0xffff & to_fp16(v[i].Im));
  898. fprintf(writef2, "%x ", 0xffff & to_fp16(v[i].Re));
  899. //printf("%f %f <> ", v[write_poz].Re, v[write_poz].Im);
  900. } */
  901.  
  902. free(v);
  903. free(scratch);
  904. fclose(readf);
  905. fclose(writef);
  906. /* fclose(writef2); */
  907. //printf("\n");
  908. }
  909.  
  910. /* int main(int argc, char** argv)
  911. {
  912.  
  913. reorder = (strcmp("REORDER_BYPASS_OFF",argv[2]) == 0);
  914. //go_fft(argv[3], atoi(argv[1]));
  915. //go_ifft(argv[4], atoi(argv[1]));
  916. go_good_dct(argv[3],atoi(argv[1]),24);
  917. //go_dct(argv[5], atoi(argv[1]));
  918. //go_idct(argv[6], atoi(argv[1]));
  919. } */
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!