API tools faq
paste
Advertisement
Guest User

Untitled

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