~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/vector/vec_combine.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)Vector maths (add, subtract, etc.)
  2  */
  3 
  4 #include <tina/sys.h>
  5 #include <tina/sysfuncs.h>
  6 #include <tina/math.h>
  7 #include <tina/mathfuncs.h>
  8 
  9 void    vector_sum_inplace(Vector * v1, Vector * v2)
 10 {
 11     int     i, n;
 12 
 13     if (v1 == NULL || v2 == NULL)
 14         return;
 15 
 16     n = MIN(v1->n, v2->n);
 17     switch (v1->vtype)
 18     {
 19     case char_v:
 20     case uchar_v:
 21     case short_v:
 22     case ushort_v:
 23     case int_v:
 24     case uint_v:
 25         for (i = 0; i < n; i++)
 26         {
 27             double     v1i, v2i;
 28 
 29             VECTOR_GET(v1, i, v1i);
 30             VECTOR_GET(v2, i, v2i);
 31             v1i += v2i;
 32             VECTOR_SET(v1, i, v1i);
 33         }
 34         break;
 35     case float_v:
 36     case double_v:
 37         for (i = 0; i < n; i++)
 38         {
 39             double  v1i, v2i;
 40 
 41             VECTOR_GET(v1, i, v1i);
 42             VECTOR_GET(v2, i, v2i);
 43             v1i += v2i;
 44             VECTOR_SET(v1, i, v1i);
 45         }
 46         break;
 47     case complex_v:
 48         for (i = 0; i < n; i++)
 49         {
 50             Complex v1i = {Complex_id};
 51             Complex v2i = {Complex_id};
 52 
 53             VECTOR_GETZ(v1, i, v1i);
 54             VECTOR_GETZ(v2, i, v2i);
 55             v1i = cmplx_sum(v1i, v2i);
 56             VECTOR_SETZ(v1, i, v1i);
 57         }
 58         break;
 59     default:
 60         break;
 61     }
 62 }
 63 
 64 Vector *vector_sum(Vector * v1, Vector * v2)
 65 {
 66     Vector *w;
 67     Vartype vtype;
 68     int     n;
 69 
 70     if (v1 == NULL && v2 == NULL)
 71         return (NULL);
 72     if (v1 == NULL)
 73         return (vector_copy(v2));
 74     if (v2 == NULL)
 75         return (vector_copy(v1));
 76 
 77     vtype = vector_sup_vtype(v1->vtype, v2->vtype);
 78     n = MAX(v1->n, v2->n);
 79     w = vector_extend(v1, n, vtype);
 80     vector_sum_inplace(w, v2);
 81 
 82     return (w);
 83 }
 84 
 85 void    vector_diff_inplace(Vector * v1, Vector * v2)
 86 {
 87     int     i, n;
 88 
 89     if (v1 == NULL || v2 == NULL)
 90         return;
 91 
 92     n = MIN(v1->n, v2->n);
 93     switch (v1->vtype)
 94     {
 95     case char_v:
 96     case uchar_v:
 97     case short_v:
 98     case ushort_v:
 99     case int_v:
100     case uint_v:
101         for (i = 0; i < n; i++)
102         {
103             double     v1i, v2i;
104 
105             VECTOR_GET(v1, i, v1i);
106             VECTOR_GET(v2, i, v2i);
107             v1i -= v2i;
108             VECTOR_SET(v1, i, v1i);
109         }
110         break;
111     case float_v:
112     case double_v:
113         for (i = 0; i < n; i++)
114         {
115             double  v1i, v2i;
116 
117             VECTOR_GET(v1, i, v1i);
118             VECTOR_GET(v2, i, v2i);
119             v1i -= v2i;
120             VECTOR_SET(v1, i, v1i);
121         }
122         break;
123     case complex_v:
124         for (i = 0; i < n; i++)
125         {
126             Complex v1i = {Complex_id};
127             Complex v2i = {Complex_id};
128 
129             VECTOR_GETZ(v1, i, v1i);
130             VECTOR_GETZ(v2, i, v2i);
131             v1i = cmplx_diff(v1i, v2i);
132             VECTOR_SETZ(v1, i, v1i);
133         }
134         break;
135     default:
136         break;
137     }
138 }
139 
140 Vector *vector_diff(Vector * v1, Vector * v2)
141 {
142     Vector *w;
143     Vartype vtype;
144     int     n;
145 
146     if (v1 == NULL && v2 == NULL)
147         return (NULL);
148     if (v1 == NULL)
149         return (vector_copy(v2));
150     if (v2 == NULL)
151         return (vector_minus(v1));
152 
153     vtype = vector_sup_vtype(v1->vtype, v2->vtype);
154     n = MAX(v1->n, v2->n);
155     w = vector_extend(v1, n, vtype);
156     vector_diff_inplace(w, v2);
157 
158     return (w);
159 }
160 
161 void    vector_prod_inplace(Vector * v1, Vector * v2)
162 {
163     int     i, n;
164 
165     if (v1 == NULL || v2 == NULL)
166         return;
167 
168     n = MIN(v1->n, v2->n);
169     switch (v1->vtype)
170     {
171     case char_v:
172     case uchar_v:
173     case short_v:
174     case ushort_v:
175     case int_v:
176     case uint_v:
177         for (i = 0; i < n; i++)
178         {
179             double  v1i, v2i;
180 
181             VECTOR_GET(v1, i, v1i);
182             VECTOR_GET(v2, i, v2i);
183             v1i *= v2i;
184             VECTOR_SET(v1, i, v1i);
185         }
186         break;
187     case float_v:
188     case double_v:
189         for (i = 0; i < n; i++)
190         {
191             double  v1i, v2i;
192 
193             VECTOR_GET(v1, i, v1i);
194             VECTOR_GET(v2, i, v2i);
195             v1i *= v2i;
196             VECTOR_SET(v1, i, v1i);
197         }
198         break;
199     case complex_v:
200         for (i = 0; i < n; i++)
201         {
202             Complex v1i = {Complex_id};
203             Complex v2i = {Complex_id};
204 
205             VECTOR_GETZ(v1, i, v1i);
206             VECTOR_GETZ(v2, i, v2i);
207             v1i = cmplx_prod(v1i, v2i);
208             VECTOR_SETZ(v1, i, v1i);
209         }
210         break;
211     default:
212         break;
213     }
214 }
215 
216 Vector *vector_prod(Vector * v1, Vector * v2)
217 {
218     Vector *w;
219     Vartype vtype;
220     int     n;
221 
222     if (v1 == NULL || v2 == NULL)
223         return (NULL);
224 
225     vtype = vector_sup_vtype(v1->vtype, v2->vtype);
226     n = MAX(v1->n, v2->n);
227     w = vector_extend(v1, n, vtype);
228     vector_prod_inplace(w, v2);
229     return (w);
230 }
231 
232 /* Product of v1 and complex conjugate of v2 */
233 void    vector_cprod_inplace(Vector * v1, Vector * v2)
234 {
235     int     i, n;
236 
237     if (v1 == NULL || v2 == NULL)
238         return;
239 
240     n = MIN(v1->n, v2->n);
241     switch (v1->vtype)
242     {
243     case char_v:
244     case uchar_v:
245     case short_v:
246     case ushort_v:
247     case int_v:
248     case uint_v:
249         for (i = 0; i < n; i++)
250         {
251             double     v1i, v2i;
252 
253             VECTOR_GET(v1, i, v1i);
254             VECTOR_GET(v2, i, v2i);
255             v1i *= v2i;
256             VECTOR_SET(v1, i, v1i);
257         }
258         break;
259     case float_v:
260     case double_v:
261         for (i = 0; i < n; i++)
262         {
263             double  v1i, v2i;
264 
265             VECTOR_GET(v1, i, v1i);
266             VECTOR_GET(v2, i, v2i);
267             v1i *= v2i;
268             VECTOR_SET(v1, i, v1i);
269         }
270         break;
271     case complex_v:
272         for (i = 0; i < n; i++)
273         {
274             Complex v1i = {Complex_id};
275             Complex v2i = {Complex_id};
276 
277             VECTOR_GETZ(v1, i, v1i);
278             VECTOR_GETZ(v2, i, v2i);
279             v1i = cmplx_cprod(v1i, v2i);
280             VECTOR_SETZ(v1, i, v1i);
281         }
282         break;
283     default:
284         break;
285     }
286 }
287 
288 Vector *vector_cprod(Vector * v1, Vector * v2)
289 {
290     Vector *w;
291     Vartype vtype;
292     int     n;
293 
294     if (v1 == NULL || v2 == NULL)
295         return (NULL);
296 
297     vtype = vector_sup_vtype(v1->vtype, v2->vtype);
298     n = MAX(v1->n, v2->n);
299     w = vector_extend(v1, n, vtype);
300     vector_cprod_inplace(w, v2);
301     return (w);
302 }
303 
304 void    vector_div_inplace(Vector * v1, Vector * v2)
305 {
306     int     i, n;
307 
308     if (v1 == NULL || v2 == NULL)
309         return;
310 
311     n = MIN(v1->n, v2->n);
312     switch (v1->vtype)
313     {
314     case char_v:
315     case uchar_v:
316     case short_v:
317     case ushort_v:
318     case int_v:
319     case uint_v:
320         for (i = 0; i < n; i++)
321         {
322             double     v1i, v2i;
323 
324             VECTOR_GET(v1, i, v1i);
325             VECTOR_GET(v2, i, v2i);
326             v1i /= v2i;
327             VECTOR_SET(v1, i, v1i);
328         }
329         break;
330     case float_v:
331     case double_v:
332         for (i = 0; i < n; i++)
333         {
334             double  v1i, v2i;
335 
336             VECTOR_GET(v1, i, v1i);
337             VECTOR_GET(v2, i, v2i);
338             v1i /= v2i;
339             VECTOR_SET(v1, i, v1i);
340         }
341         break;
342     case complex_v:
343         for (i = 0; i < n; i++)
344         {
345             Complex v1i = {Complex_id};
346             Complex v2i = {Complex_id};
347 
348             VECTOR_GETZ(v1, i, v1i);
349             VECTOR_GETZ(v2, i, v2i);
350             v1i = cmplx_div(v1i, v2i);
351             VECTOR_SETZ(v1, i, v1i);
352         }
353         break;
354     default:
355         break;
356     }
357 }
358 
359 Vector *vector_div(Vector * v1, Vector * v2)
360 {
361     Vector *w;
362     Vartype vtype;
363     int     n;
364 
365     if (v1 == NULL || v2 == NULL)
366         return (NULL);
367 
368     vtype = vector_sup_vtype(v1->vtype, v2->vtype);
369     n = MAX(v1->n, v2->n);
370     w = vector_extend(v1, n, vtype);
371     vector_div_inplace(w, v2);
372     return (w);
373 }
374 
375 Vector *vector_conv(Vector * v1, Vector * v2)
376 {
377     Vector *f1;
378     Vector *f2;
379     Vector *fc;
380     Vector *c;
381 
382     if (v1 == NULL || v2 == NULL)
383         return (NULL);
384 
385     f1 = vector_fft(v1);
386     f2 = vector_fft(v2);
387     fc = vector_prod(v1, v2);
388     vector_free(f1);
389     vector_free(f2);
390     c = vector_fft_inverse(fc);
391     vector_free(fc);
392     return (c);
393 }
394 
395 /* un-normalised correlation of two vectors calculated by fft */
396 Vector *vector_corr(Vector * v1, Vector * v2)
397 {
398     Vector *f1;
399     Vector *f2;
400     Vector *fc;
401     Vector *c;
402     Vector *cr;
403     Vector *e1 = NULL;
404     Vector *e2 = NULL;
405     int     n = 1;
406 
407     if (v1 == NULL || v2 == NULL)
408         return (NULL);
409     while (n < v1->n || n < v2->n)
410         n *= 2;;
411 
412     if (n != v1->n)
413         v1 = e1 = vector_extend(v1, n, v1->vtype);
414     if (n != v2->n)
415         v2 = e2 = vector_extend(v2, n, v2->vtype);
416 
417     f1 = vector_fft(v1);
418     f2 = vector_fft(v2);
419     fc = vector_cprod(f1, f2);
420     vector_free(f1);
421     vector_free(f2);
422     c = vector_fft_inverse(fc);
423     vector_free(fc);
424     cr = vector_cast(c, double_v);
425     vector_free(c);
426     vector_free(e1);
427     vector_free(e2);
428     return (cr);
429 }
430 
431 /* un-normalised correlation of two vectors calculated by fft, with cut
432  * lowest frequencies cut off */
433 Vector *vector_corr_cutoff(Vector * v1, Vector * v2, int cut)
434 {
435     Vector *f1;
436     Vector *f2;
437     Vector *fc;
438     Vector *c;
439     Vector *cr;
440     Vector *e1 = NULL;
441     Vector *e2 = NULL;
442     int     n = 1, i;
443 
444     if (v1 == NULL || v2 == NULL)
445         return (NULL);
446     while (n < v1->n || n < v2->n)
447         n *= 2;;
448 
449     if (n != v1->n)
450         v1 = e1 = vector_extend(v1, v1->vtype, n);
451     if (n != v2->n)
452         v2 = e2 = vector_extend(v2, v2->vtype, n);
453 
454     f1 = vector_fft(v1);
455     f2 = vector_fft(v2);
456     for (i = 0; i < cut; i++)
457     {
458         VECTOR_SETZ(f1, i, cmplx_zero());
459         if (i != 0)
460             VECTOR_SETZ(f1, n - i, cmplx_zero());
461 
462         VECTOR_SETZ(f2, i, cmplx_zero());
463         if (i != 0)
464             VECTOR_SETZ(f2, n - i, cmplx_zero());
465     }
466     fc = vector_cprod(f1, f2);
467     vector_free(f1);
468     vector_free(f2);
469     c = vector_fft_inverse(fc);
470     vector_free(fc);
471     cr = vector_cast(c, double_v);
472     vector_free(c);
473     vector_free(e1);
474     vector_free(e2);
475     return (cr);
476 }
477 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.