1 /*
2     dcblockr.c:
3 
4     Copyright (C) 1998 John ffitch, 2008 V Lazzarini
5 
6     This file is part of Csound.
7 
8     The Csound Library is free software; you can redistribute it
9     and/or modify it under the terms of the GNU Lesser General Public
10     License as published by the Free Software Foundation; either
11     version 2.1 of the License, or (at your option) any later version.
12 
13     Csound is distributed in the hope that it will be useful,
14     but WITHOUT ANY WARRANTY; without even the implied warranty of
15     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16     GNU Lesser General Public License for more details.
17 
18     You should have received a copy of the GNU Lesser General Public
19     License along with Csound; if not, write to the Free Software
20     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
21     02110-1301 USA
22 */
23 
24 /*******************************************/
25 /*  DC Blocking Filter                     */
26 /*  by Perry R. Cook, 1995-96              */
27 /*  This guy is very helpful in, uh,       */
28 /*  blocking DC.  Needed because a simple  */
29 /*  low-pass reflection filter allows DC   */
30 /*  to build up inside recursive           */
31 /*  structures.                            */
32 /*******************************************/
33 
34 #include "stdopcod.h"
35 #include "dcblockr.h"
36 
dcblockrset(CSOUND * csound,DCBlocker * p)37 static int32_t dcblockrset(CSOUND *csound, DCBlocker* p)
38 {
39     IGN(csound);
40     p->outputs = 0.0;
41     p->inputs = 0.0;
42     p->gain = (double)*p->gg;
43     if (p->gain == 0.0 || p->gain>=1.0 || p->gain<=-1.0)
44       p->gain = 0.99;
45     return OK;
46 }
47 
dcblockr(CSOUND * csound,DCBlocker * p)48 static int32_t dcblockr(CSOUND *csound, DCBlocker* p)
49 {
50     IGN(csound);
51     MYFLT       *ar = p->ar;
52     uint32_t    offset = p->h.insdshead->ksmps_offset;
53     uint32_t    early  = p->h.insdshead->ksmps_no_end;
54     uint32_t    n, nsmps = CS_KSMPS;
55     double      gain = p->gain;
56     double      outputs = p->outputs;
57     double      inputs = p->inputs;
58     MYFLT       *samp = p->in;
59 
60     if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
61     if (UNLIKELY(early)) {
62       nsmps -= early;
63       memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
64     }
65     for (n=offset; n<nsmps; n++) {
66       double sample = (double)samp[n];
67       outputs = sample - inputs + (gain * outputs);
68       inputs = sample;
69       ar[n] = (MYFLT)outputs;
70     }
71     p->outputs = outputs;
72     p->inputs = inputs;
73     return OK;
74 }
75 
76 /*******************************************/
77 /*  DC Blocking Filter                     */
78 /*  Improved DC attenuation                */
79 /*  V Lazzarini                            */
80 /*******************************************/
81 
82 typedef struct _dcblk2 {
83   OPDS    h;
84   MYFLT   *output;
85   MYFLT   *input, *order, *iskip;
86   AUXCH   delay1;
87   AUXCH   iirdelay1, iirdelay2, iirdelay3, iirdelay4;
88   double  ydels[4];
89   int32_t dp1,dp2;
90   double  scaler;
91 } DCBlock2;
92 
93 
dcblock2set(CSOUND * csound,DCBlock2 * p)94 static int32_t dcblock2set(CSOUND *csound, DCBlock2* p)
95 {
96     int32_t order = (int32_t) *p->order;
97     if (order == 0) order = 128;
98     else if (order < 4) order = 4;
99 
100     if (p->delay1.auxp == NULL ||
101         p->delay1.size < (order-1)*2*sizeof(double))
102       csound->AuxAlloc(csound, (order-1)*2*sizeof(double),
103                        &p->delay1);
104 
105     if (p->iirdelay1.auxp == NULL ||
106         p->iirdelay1.size < (order)*sizeof(double))
107       csound->AuxAlloc(csound,
108                        (order)*sizeof(double), &p->iirdelay1);
109 
110     if (p->iirdelay2.auxp == NULL ||
111         p->iirdelay2.size < (order)*sizeof(double))
112       csound->AuxAlloc(csound,
113                        (order)*sizeof(double), &p->iirdelay2);
114 
115     if (p->iirdelay3.auxp == NULL ||
116         p->iirdelay3.size < (order)*sizeof(double))
117       csound->AuxAlloc(csound,
118                        (order)*sizeof(double), &p->iirdelay3);
119 
120     if (p->iirdelay4.auxp == NULL ||
121         p->iirdelay4.size < (order)*sizeof(double))
122       csound->AuxAlloc(csound,
123                        (order)*sizeof(double), &p->iirdelay4);
124 
125     p->scaler = 1.0/order;
126     if (!*p->iskip) {
127       memset(p->ydels, 0, 4*sizeof(double));
128       /* p->ydels[0] = 0.0;   p->ydels[1] = 0.0; */
129       /* p->ydels[2] = 0.0;   p->ydels[3] = 0.0; */
130       memset(p->delay1.auxp, 0, sizeof(double)*(order-1)*2);
131       memset(p->iirdelay1.auxp, 0, sizeof(double)*(order));
132       memset(p->iirdelay2.auxp, 0, sizeof(double)*(order));
133       memset(p->iirdelay3.auxp, 0, sizeof(double)*(order));
134       memset(p->iirdelay4.auxp, 0, sizeof(double)*(order));
135       p->dp1 = 0; p->dp2 = 0;
136     }
137     return OK;
138 }
139 
dcblock2(CSOUND * csound,DCBlock2 * p)140 static int32_t dcblock2(CSOUND *csound, DCBlock2* p)
141 {
142     IGN(csound);
143     uint32_t offset = p->h.insdshead->ksmps_offset;
144     uint32_t early  = p->h.insdshead->ksmps_no_end;
145     uint32_t i, nsmps = CS_KSMPS;
146     MYFLT    *in = p->input;
147     MYFLT    *out = p->output;
148     double   *del1 = (double *)p->delay1.auxp;
149     double   *iirdel[4],x1,x2,y,del;
150     double   *ydels = p->ydels;
151     double   scale = p->scaler;
152     int32_t      p1 = p->dp1;
153     int32_t      p2 = p->dp2;
154     int32_t      j,del1size, iirdelsize;
155 
156     iirdel[0] = (double *) p->iirdelay1.auxp;
157     iirdel[1] = (double *) p->iirdelay2.auxp;
158     iirdel[2] = (double *) p->iirdelay3.auxp;
159     iirdel[3] = (double *) p->iirdelay4.auxp;
160 
161     del1size = p->delay1.size/sizeof(double);
162     iirdelsize = p->iirdelay1.size/sizeof(double);
163 
164     if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
165     if (UNLIKELY(early)) {
166       nsmps -= early;
167       memset(&out[nsmps], '\0', early*sizeof(MYFLT));
168     }
169     for (i=offset; i < nsmps; i++) {
170 
171       /* long delay */
172       del = del1[p1];
173       del1[p1] = x1 = (double)in[i];
174 
175       /* IIR cascade */
176       for (j=0; j < 4; j++) {
177         x2 = iirdel[j][p2];
178         iirdel[j][p2] = x1;
179         y = x1 - x2 + ydels[j];
180         ydels[j] = y;
181         x1 = y*scale;
182       }
183       out[i] = (MYFLT)(del - x1);
184 
185       p1 = (p1 == del1size-1 ? 0 : p1 + 1);
186       p2 = (p2 == iirdelsize-1 ? 0 : p2 + 1);
187     }
188 
189     p->dp1 = p1; p->dp2 = p2;
190     return OK;
191 }
192 
193 
194 
195 
196 #define S(x)    sizeof(x)
197 
198 static OENTRY localops[] = {
199   { "dcblock", S(DCBlocker), 0, 3, "a", "ao",
200                                    (SUBR)dcblockrset, (SUBR)dcblockr},
201   { "dcblock2", S(DCBlock2), 0, 3, "a", "aoo",
202                                    (SUBR)dcblock2set, (SUBR)dcblock2}
203 };
204 
dcblockr_init_(CSOUND * csound)205 int32_t dcblockr_init_(CSOUND *csound)
206 {
207     return csound->AppendOpcodes(csound, &(localops[0]),
208                                  (int32_t
209                                   ) (sizeof(localops) / sizeof(OENTRY)));
210 }
211