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)7 LTFAT_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