File: earplug~.c

package info (click to toggle)
pd-earplug 0.2-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 2,080 kB
  • sloc: ansic: 2,096; makefile: 487; perl: 27
file content (259 lines) | stat: -rw-r--r-- 9,327 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
/* RT binaural filter: earplug~        */
/* based on KEMAR impulse measurement  */
/* Pei Xiang, summer 2004              */
/* Revised in Fall 2006 by Jorge Castellanos */
/* Revised in Spring 2009 by Hans-Christoph Steiner to compile in the data file */

#include <stdio.h>
#include <math.h>
#include "m_pd.h"
#include "earplug~.h"
#ifdef _MSC_VER /* Thes pragmas only apply to Microsoft's compiler */
#pragma warning( disable : 4244 )
#pragma warning( disable : 4305 )
#endif

/* elevation degree:      -40  -30  -20  -10   0   10  20  30  40  50  60  70  80  90   */
/* index array:             0    1    2    3   4    5   6   7   8   9  10  11  12  13   */
/* impulse response number: 29   31   37   37  37   37  37  31  29  23  19  13   7   1   */ 
/* 0 degree response index:  0   29   60   97 134  171 208 245 276 305 328 347 360 367   */

static t_class *earplug_class;

typedef struct _earplug
{
    t_object x_obj; 
    t_outlet *left_channel ; 
    t_outlet *right_channel ; 

    t_float azimuth ;             /* from 0 to 360 degrees */
    t_float elevation ;           /* from -40 to 90 (degrees) */
    t_float azi ;
    t_float ele ;
     
    t_float crossCoef[8192] ; 
    t_float azimScale[13] ;
    t_int azimOffset[13] ;

    t_float previousImpulse[2][128] ;
    t_float currentImpulse[2][128] ;
    t_float convBuffer[128] ;
    t_float (*impulses)[2][128];          /* a 3D array of 368x2x128 */
    t_float f ;                   /* dummy float for dsp */
    int bufferPin;
} t_earplug;

static t_int *earplug_perform(t_int *w)
{
    t_earplug *x = (t_earplug *)(w[1]) ;
    t_float *in = (t_float *)(w[2]);
    t_float *right_out = (t_float *)(w[3]);
    t_float *left_out = (t_float *)(w[4]); 
    int blocksize = (int)(w[5]);

    unsigned ch_L, ch_R, i;

    x->azi = x->azimuth ; 
    x->ele = x->elevation ;
	
    if (x->ele < -40) 
		x->ele = -40;
    if (x->ele > 90)
		x->ele = 90;
    if (x->azi < 0 || x->azi > 360) 
		x->azi = 0;
    if (x->azi <= 180)
    {
		ch_L = 0;
		ch_R = 1;
	}
    else
    { 
		ch_L = 1;
		ch_R = 0;
		x->azi = 360.0 - x->azi;
	}
	
    x->ele *= 0.1;   /* divided by 10 since each elevation is 10 degrees apart */

    if (x->ele < 8.) // if elevation is less than 80 degrees...
    { 
		int elevInt = (int)floor(x->ele); // A quantized version of the elevation    
		unsigned elevGridIndex = elevInt + 4;  // Used as the index to the array of scaling factors for the azimuth (adding 4 because the lowest elevation is -4, so it starts at 0)
		unsigned azimIntUp = (unsigned)(x->azi * x->azimScale[elevGridIndex+1]); // 
        float azimFracUp = azimIntUp + 1.0 - x->azi * x->azimScale[elevGridIndex+1];
		float azimFracUpInv = 1.0 - azimFracUp;
		float elevFracUp = x->ele - elevInt * 1.0;
    	unsigned azimIntDown = (unsigned)(x->azi * x->azimScale[elevGridIndex]);
        float azimFracDown = azimIntDown + 1.0 - x->azi * x->azimScale[elevGridIndex];
		float azimFracDownInv = 1.0 - azimFracDown;
		float elevFracDown = 1.0 - elevFracUp;
		unsigned lowerIdx = x->azimOffset[elevGridIndex] + azimIntDown;
		unsigned upperIdx = x->azimOffset[elevGridIndex + 1] + azimIntUp;
		
		for (i = 0; i < 128; i++)
		{
			x->currentImpulse[ch_L][i] = elevFracDown *		// Interpolate the lower two HRIRs and multiply them by their "fraction"
										(azimFracDown * x->impulses[lowerIdx][0][i] + 
										azimFracDownInv * x->impulses[lowerIdx + 1][0][i]) + 
										elevFracUp *		// Interpolate the upper two HRIRs and multiply them by their "fraction"
										(azimFracUp * x->impulses[upperIdx][0][i] + 
										azimFracUpInv * x->impulses[upperIdx + 1][0][i]);   

			x->currentImpulse[ch_R][i] = elevFracDown *
										(azimFracDown * x->impulses[lowerIdx][1][i] + 
										azimFracDownInv * x->impulses[lowerIdx + 1][1][i]) + 
										elevFracUp * 
										(azimFracUp * x->impulses[upperIdx][1][i] + 
										azimFracUpInv * x->impulses[upperIdx + 1][1][i]);  
		}   

    }
    else	// if elevation is 80 degrees or more the interpolation requires only three points (because there's only one HRIR at 90 deg)
    {
    	unsigned azimIntDown = (unsigned)(x->azi * 0.033333); // Scale the azimuth to 12 (the number of HRIRs at 80 deg) discreet points
        float azimFracDown = azimIntDown + 1.0 - x->azi * 0.033333;
        float elevFracUp = x->ele - 8.0;
		float elevFracDown = 9.0 - x->ele;
		for (i = 0; i < 128; i++) {
			
			x->currentImpulse[ch_L][i] = elevFracDown * 
										( azimFracDown * x->impulses[360+azimIntDown][0][i] +	// These two lines interpolate the lower two HRIRs
										(1.0 - azimFracDown) * x->impulses[361+azimIntDown][0][i] )
										+ elevFracUp * x->impulses[367][0][i];	// multiply the 90 degree HRIR with its corresponding fraction
			x->currentImpulse[ch_R][i] = elevFracDown * 
										(azimFracDown * x->impulses[360+azimIntDown][1][i]  +
										(1.0 - azimFracDown) * x->impulses[361+azimIntDown][1][i])
										+ elevFracUp * x->impulses[367][1][i]; 
		}

	}

    float inSample;
	float convSum[2]; // to accumulate the sum during convolution.
    int blockScale = 8192 / blocksize;

	// Convolve the - interpolated - HRIRs (Left and Right) with the input signal.
    while (blocksize--)
    {
		convSum[0] = 0; 
		convSum[1] = 0;	

		inSample = *(in++);

		x->convBuffer[x->bufferPin] = inSample;
		unsigned scaledBlocksize = blocksize * blockScale;
		unsigned blocksizeDelta = 8191 - scaledBlocksize;
		for (i = 0; i < 128; i++)
		{ 
			convSum[0] += (x->previousImpulse[0][i] * x->crossCoef[blocksizeDelta] + 
							x->currentImpulse[0][i] * x->crossCoef[scaledBlocksize]) * 
							x->convBuffer[(x->bufferPin - i) &127];
			convSum[1] += (x->previousImpulse[1][i] * x->crossCoef[blocksizeDelta] + 
							x->currentImpulse[1][i] * x->crossCoef[scaledBlocksize]) * 
							x->convBuffer[(x->bufferPin - i) &127];

			x->previousImpulse[0][i] = x->currentImpulse[0][i];
			x->previousImpulse[1][i] = x->currentImpulse[1][i];
		}	
		x->bufferPin = (x->bufferPin + 1) & 127;
  
        *left_out++ = convSum[0];
     	*right_out++ = convSum[1];
    }
    return (w+6);
}

static void earplug_dsp(t_earplug *x, t_signal **sp)
{
		// callback, params, userdata, in_samples, out_L,		out_R,		blocksize.
    dsp_add(earplug_perform, 5, x,  sp[0]->s_vec, sp[1]->s_vec, sp[2]->s_vec, sp[0]->s_n);
}

static void *earplug_new(t_floatarg azimArg, t_floatarg elevArg)
{
    t_earplug *x = (t_earplug *)pd_new(earplug_class);
    x->left_channel = outlet_new(&x->x_obj, gensym("signal"));
    x->right_channel = outlet_new(&x->x_obj, gensym("signal"));
    floatinlet_new(&x->x_obj, &x->azimuth) ;     /* 0 to 360 degrees */
    floatinlet_new(&x->x_obj, &x->elevation) ;   /* -40 to 90 degrees */

    x->azimuth = azimArg ;
    x->elevation = elevArg ;

    int i, j;
    FILE *fp;
    char buff[1024], *bufptr;
    int filedesc;

	filedesc = open_via_path(canvas_getdir(canvas_getcurrent())->s_name, "earplug_data.txt", "", buff, &bufptr, 1024, 0 );

    if (filedesc >= 0) // If there was no error opening the text file...
    {	
		post("[earplug~] found impulse response file, overriding defaults:") ;
        post("let's try loading %s/earplug_data.txt", buff);
        fp = fdopen(filedesc, "r") ;  
        for (i = 0; i < 368; i++) 
        {
            while(fgetc(fp) != 10) ;
            for (j = 0 ; j < 128 ; j++) {
				fscanf(fp, "%f %f ", &earplug_impulses[i][0][j], &earplug_impulses[i][1][j]);
            }
        }
        fclose(fp) ;
    }
    x->impulses = earplug_impulses;
    
    post("        earplug~: binaural filter with measured responses\n") ;
    post("        elevation: -40 to 90 degrees. azimuth: 360") ;
    post("        don't let blocksize > 8192\n"); 
       
    for (i = 0; i < 128 ; i++)
    {
         x->convBuffer[i] = 0; 
		 x->previousImpulse[0][i] = 0; 
		 x->previousImpulse[1][i] = 0;
    }
	
    x->bufferPin = 0;

    for (i = 0; i < 8192 ; i++)
    {	
		x->crossCoef[i] = 1.0 * i / 8192;
	}

	// This is the scaling factor for the azimuth so that it corresponds to an HRTF in the KEMAR database
    x->azimScale[0] = 0.153846153; x->azimScale[8] = 0.153846153;   /* -40 and 40 degree */
    x->azimScale[1] = 0.166666666; x->azimScale[7] = 0.166666666;   /* -30 and 30 degree */
    x->azimScale[2] = 0.2; x->azimScale[3]=0.2; x->azimScale[4]=0.2; x->azimScale[5]=0.2; x->azimScale[6]=0.2;  /* -20 to 20 degree */
    x->azimScale[9] = 0.125;  		/* 50 degree */
    x->azimScale[10] = 0.1;		/* 60 degree */
    x->azimScale[11] = 0.066666666;      /* 70 degree */
    x->azimScale[12] = 0.033333333;	/* 80 degree */

    x->azimOffset[0] = 0 ; 
    x->azimOffset[1] = 29 ;
    x->azimOffset[2] = 60 ;
    x->azimOffset[3] = 97 ;
    x->azimOffset[4] = 134 ;
    x->azimOffset[5] = 171 ;
    x->azimOffset[6] = 208 ;
    x->azimOffset[7] = 245 ;
    x->azimOffset[8] = 276 ;
    x->azimOffset[9] = 305 ;
    x->azimOffset[10] = 328 ;
    x->azimOffset[11] = 347 ;
    x->azimOffset[12] = 360 ;

    return (x);
}

void earplug_tilde_setup(void)
{
    earplug_class = class_new(gensym("earplug~"), (t_newmethod)earplug_new, 0,
    	sizeof(t_earplug), CLASS_DEFAULT, A_DEFFLOAT, A_DEFFLOAT, 0);
   
    CLASS_MAINSIGNALIN(earplug_class, t_earplug, f);
   
    class_addmethod(earplug_class, (t_method)earplug_dsp, gensym("dsp"), A_CANT, 0);
}