Adding sub-band CELP (SB-CELP) -like encoding. Still incomplete.
[speexdsp.git] / libspeex / filters.c
1 /* Copyright (C) 2002 Jean-Marc Valin 
2    File: filters.c
3    Various analysis/synthesis filters
4
5    This library is free software; you can redistribute it and/or
6    modify it under the terms of the GNU Lesser General Public
7    License as published by the Free Software Foundation; either
8    version 2.1 of the License, or (at your option) any later version.
9    
10    This library is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    Lesser General Public License for more details.
14    
15    You should have received a copy of the GNU Lesser General Public
16    License along with this library; if not, write to the Free Software
17    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18 */
19
20 #include "filters.h"
21
22 #define min(a,b) ((a) < (b) ? (a) : (b))
23
24 void bw_lpc(float gamma, float *lpc_in, float *lpc_out, int order)
25 {
26    int i;
27    float tmp=1;
28    for (i=0;i<order+1;i++)
29    {
30       lpc_out[i] = tmp * lpc_in[i];
31       tmp *= gamma;
32    }
33 }
34
35 void syn_filt(float *x, float *a, float *y, int N, int ord)
36 {
37    int i,j;
38    for (i=0;i<N;i++)
39    {
40       y[i]=x[i];
41       for (j=1;j<=ord;j++)
42          y[i] -= a[j]*y[i-j];
43    }
44 }
45
46 void syn_filt_zero(float *x, float *a, float *y, int N, int ord)
47 {
48    int i,j;
49    for (i=0;i<N;i++)
50    {
51       y[i]=x[i];
52       for (j=1;j<=min(ord,i);j++)
53          y[i] -= a[j]*y[i-j];
54    }
55 }
56
57 void syn_filt_mem(float *x, float *a, float *y, int N, int ord, float *mem)
58 {
59    int i,j;
60    for (i=0;i<N;i++)
61    {
62       y[i]=x[i];
63       for (j=1;j<=min(ord,i);j++)
64          y[i] -= a[j]*y[i-j];
65       for (j=i+1;j<=ord;j++)
66          y[i] -= a[j]*mem[j-i-1];
67    }
68    for (i=0;i<ord;i++)
69       mem[i]=y[N-i-1];
70 }
71
72
73 void residue(float *x, float *a, float *y, int N, int ord)
74 {
75    int i,j;
76    for (i=N-1;i>=0;i--)
77    {
78       y[i]=x[i];
79       for (j=1;j<=ord;j++)
80          y[i] += a[j]*x[i-j];
81    }
82 }
83
84 void residue_zero(float *x, float *a, float *y, int N, int ord)
85 {
86    int i,j;
87    for (i=N-1;i>=0;i--)
88    {
89       y[i]=x[i];
90       for (j=1;j<=min(ord,i);j++)
91          y[i] += a[j]*x[i-j];
92    }
93 }
94
95 void residue_mem(float *x, float *a, float *y, int N, int ord, float *mem)
96 {
97    int i,j;
98    for (i=N-1;i>=0;i--)
99    {
100       y[i]=x[i];
101       for (j=1;j<=min(ord,i);j++)
102          y[i] += a[j]*x[i-j];
103       for (j=i+1;j<=ord;j++)
104          y[i] += a[j]*mem[j-i-1];
105    }
106    for (i=0;i<ord;i++)
107       mem[i]=x[N-i-1];
108 }
109
110 float xcorr(float *x, float *y, int len)
111 {
112    int i;
113    float sum=0;
114    for (i=0;i<len;i++)
115       sum += x[i]*y[i];
116    return sum;
117 }
118
119 void fir_mem(float *x, float *a, float *y, int N, int M, float *mem)
120 {
121    int i,j;
122    for (i=0;i<N;i++)
123    {
124       y[i]=0;
125       for (j=0;j<=min(M-1,i);j++)
126          y[i] += a[j]*x[i-j];
127       for (j=i+1;j<=M-1;j++)
128          y[i] += a[j]*mem[j-i-1];
129    }
130    for (i=0;i<M-1;i++)
131       mem[i]=x[N-i-1];
132 }