1 #include "ltfat.h" 2 #include "ltfat/types.h" 3 #include "ltfat/macros.h" 4 5 6 LTFAT_API void LTFAT_NAME(pfilt_fir_rr)7LTFAT_NAME(pfilt_fir_rr)(const LTFAT_REAL *f, const LTFAT_REAL *g, 8 ltfat_int L, ltfat_int gl, 9 ltfat_int W, ltfat_int a, 10 LTFAT_REAL *cout) 11 { 12 /* --------- initial declarations -------------- */ 13 14 ltfat_int l, n, w; 15 16 LTFAT_REAL *gw; 17 18 LTFAT_REAL *gb; 19 LTFAT_REAL fw; 20 21 const LTFAT_REAL *fbd; 22 23 /* ----------- calculation of parameters and plans -------- */ 24 25 ltfat_int N=L/a; 26 27 /* These are floor operations. */ 28 ltfat_int glh=gl/2; 29 30 /* This is a ceil operation. */ 31 ltfat_int glh_d_a=(ltfat_int)ceil((glh*1.0)/(a)); 32 33 gw = (LTFAT_REAL*)ltfat_malloc(gl*sizeof(LTFAT_REAL)); 34 35 /* ---------- main body ----------- */ 36 37 /* Do the fftshift of g to place the center in the middle and 38 * conjugate it. 39 */ 40 41 for (l=0; l<glh; l++) 42 { 43 gw[l]=g[l+(gl-glh)]; 44 } 45 for (l=glh; l<gl; l++) 46 { 47 gw[l]=g[l-glh]; 48 } 49 50 for (w=0; w<W; w++) 51 { 52 /*----- Handle the first boundary using periodic boundary conditions.*/ 53 for (n=0; n<glh_d_a; n++) 54 { 55 gb=gw; 56 fbd=f+L-(glh-n*a)+L*w; 57 fw=0.0; 58 for (l=0; l<glh-n*a; l++) 59 { 60 fw +=fbd[l]*gb[l]; 61 } 62 fbd=f-(glh-n*a)+L*w; 63 for (l=glh-n*a; l<gl; l++) 64 { 65 fw +=fbd[l]*gb[l]; 66 } 67 cout[n+w*N]=fw; 68 } 69 70 /* ----- Handle the middle case. --------------------- */ 71 for (n=glh_d_a; n<(L-(gl+1)/2)/a+1; n++) 72 { 73 gb=gw; 74 fbd=f+(n*a-glh+L*w); 75 fw=0; 76 for (l=0; l<gl; l++) 77 { 78 fw +=fbd[l]*gb[l]; 79 } 80 cout[n+w*N]=fw; 81 } 82 83 /* Handle the last boundary using periodic boundary conditions. */ 84 for (n=(L-(gl+1)/2)/a+1; n<N; n++) 85 { 86 gb=gw; 87 fbd=f+(n*a-glh+L*w); 88 fw=0; 89 for (l=0; l<L-n*a+glh; l++) 90 { 91 fw +=fbd[l]*gb[l]; 92 } 93 fbd=f-(L-n*a+glh)+L*w; 94 for (l=L-n*a+glh; l<gl; l++) 95 { 96 fw +=fbd[l]*gb[l]; 97 } 98 cout[n+w*N]=fw; 99 } 100 } 101 102 /* ----------- Clean up ----------------- */ 103 ltfat_free(gw); 104 105 } 106