~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/sandbox/cdavid/src/autocorr_nofft.tpl

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
[+ AutoGen5 template c +]
 
2
/*
 
3
 * Last Change: Tue Nov 28 03:00 PM 2006 J
 
4
 * vim:syntax=c
 
5
 *
 
6
 * TODO: is size_t 64 bits long on 64 bits machines ?
 
7
 */
 
8
#include <stddef.h> /* for size_t */
 
9
#include <stdio.h> /* for size_t */
 
10
 
 
11
#include "autocorr_nofft.h"
 
12
 
 
13
/*
 
14
 * NOFFT auto correlation
 
15
 */
 
16
 
 
17
[+ For float_type +]
 
18
/* 
 
19
 * [+ (get "type_name") +] version; out must have a size of lag+1, already pre allocated 
 
20
 *
 
21
 * lag should be < size
 
22
 *
 
23
 * returns 0 is succesfull, other value otherwise.
 
24
 */
 
25
int [+ (get "short_name") +]_xcorr_nofft_1d(const [+ (get "type_name") +] *in, 
 
26
    const size_t size, [+ (get "type_name") +] *out, const size_t lag)
 
27
{
 
28
        size_t  i, j;
 
29
        [+ (get "type_name") +] acc;
 
30
        
 
31
        /* lag 0 */
 
32
        acc     = 0;
 
33
        for (i = 0; i <  size; ++i) {
 
34
                acc     += in[i]*in[i];
 
35
        }
 
36
        out[0] = acc;
 
37
 
 
38
        /* lag : 1 -> lag */
 
39
        for (i = 1; i <= lag; ++i) {
 
40
                acc     = 0;
 
41
                for (j = i; j < size; ++j) {
 
42
                        acc     +=      in[j-i]*in[j];
 
43
                }
 
44
                out[i]  = acc;
 
45
        }
 
46
 
 
47
        return 0;
 
48
}
 
49
[+ ENDFOR float_type +]
 
50
 
 
51
[+ For float_type +]
 
52
/* 
 
53
 * [+ (get "type_name") +] version for non contiguous arrays; the corresponding 
 
54
 * array should have at least in_size elements. 
 
55
 *
 
56
 * Constraints:
 
57
 *  - lag should be < in_size
 
58
 *  - strides in bytes
 
59
 *  - TODO: check if should be aligned ?
 
60
 * 
 
61
 * returns 0 is succesfull, other value otherwise.
 
62
 */
 
63
int [+ (get "short_name") +]_xcorr_nofft_1d_noncontiguous(const [+ (get "type_name") +] *in, size_t in_size, 
 
64
    size_t in_stride, [+ (get "type_name") +] *out, size_t out_stride, size_t lag)
 
65
{
 
66
        size_t  i, j, clag;
 
67
    size_t  istride = in_stride / sizeof([+ (get "type_name") +]);
 
68
    size_t  ostride = out_stride / sizeof([+ (get "type_name") +]);
 
69
        [+ (get "type_name") +] acc;
 
70
        
 
71
        /* lag 0 */
 
72
        acc     = 0;
 
73
        for (i = 0; i < in_size * istride; i+= istride) {
 
74
                acc     += in[i]*in[i];
 
75
        }
 
76
        out[0] = acc;
 
77
 
 
78
        /* lag : 1 -> lag */
 
79
        for (i = 1; i <= lag ; ++i) {
 
80
                acc         = 0;
 
81
        clag    = i * istride;
 
82
        for (j = clag; j < in_size * istride; j += istride) {
 
83
                        acc     +=      in[j-clag]*in[j];
 
84
                }
 
85
                out[i * ostride]        = acc;
 
86
        }
 
87
 
 
88
    return 0;
 
89
}
 
90
[+ ENDFOR float_type +]
 
91
 
 
92
[+ For float_type +]
 
93
/* 
 
94
 * For rank 2 arrays, contiguous cases
 
95
 * [+ (get "type_name") +] version; out must have a size of lag+1, already pre allocated 
 
96
 *
 
97
 * lag should be < size
 
98
 *
 
99
 * returns 0 is succesfull, other value otherwise.
 
100
 */
 
101
int [+ (get "short_name") +]_xcorr_nofft_2d(const [+ (get "type_name") +] *in, 
 
102
    size_t dim0, size_t dim1, [+ (get "type_name") +] *out, const size_t lag)
 
103
{
 
104
    size_t  i;
 
105
    [+ (get "type_name") +]  *coaxis;
 
106
 
 
107
#if 0
 
108
    for(i = 0; i < dim0; ++i) {
 
109
        fprintf(stdout, "%d 1d autocorr, first element is %f\n", i, in[i * dim1]);
 
110
    }
 
111
#endif
 
112
    for(i = 0; i < dim0; ++i) {
 
113
        coaxis  = out + i * (lag + 1);
 
114
        [+ (get "short_name") +]_xcorr_nofft_1d(in + i * dim1, dim1, coaxis, lag);
 
115
    }
 
116
 
 
117
        return 0;
 
118
}
 
119
[+ ENDFOR float_type +]
 
120
 
 
121
[+ For float_type +]
 
122
/* 
 
123
 * For rank 2 arrays, non contiguous cases
 
124
 * [+ (get "type_name") +] version; out must have a size of lag+1, already pre allocated 
 
125
 *
 
126
 * lag should be < size
 
127
 *
 
128
 * returns 0 is succesfull, other value otherwise.
 
129
 */
 
130
int [+ (get "short_name") +]_xcorr_nofft_2d_noncontiguous(const [+ (get "type_name") +] *in, 
 
131
    size_t dim0, size_t dim1, size_t in_stride0, size_t in_stride1, 
 
132
    [+ (get "type_name") +] *out, size_t out_stride0, size_t out_stride1,
 
133
    const size_t lag)
 
134
{
 
135
    size_t  i;
 
136
 
 
137
    size_t  istride0    = in_stride0 / sizeof([+ (get "type_name") +]);
 
138
    size_t  ostride0    = out_stride0 / sizeof([+ (get "type_name") +]);
 
139
 
 
140
    [+ (get "type_name") +]  *coaxis;
 
141
#if 0
 
142
    fprintf(stdout, "%s: shape is (%d, %d)\n", __func__, dim0, dim1);
 
143
    fprintf(stdout, "%s: istrides are (%d, %d)\n", __func__, istride0, istride1);
 
144
    
 
145
    fprintf(stdout, "%s: ostrides are (%d, %d)\n", __func__, ostride0, ostride1);
 
146
    for(i = 0; i < dim0; ++i) {
 
147
        ciaxis  = in + i * istride0;
 
148
        coaxis  = out + i * istride0;
 
149
        fprintf(stdout, "%d 1d autocorr, first element is %f, last is %f (%d el)\n", 
 
150
            i, ciaxis[0], ciaxis[(dim1-1) * istride1], dim1);
 
151
    }
 
152
#endif
 
153
 
 
154
    for(i = 0; i < dim0; ++i) {
 
155
        coaxis  = out + i * ostride0;
 
156
        [+ (get "short_name") +]_xcorr_nofft_1d_noncontiguous(in + i * istride0, dim1, in_stride1,
 
157
            coaxis, out_stride1, lag);
 
158
    }
 
159
        return 0;
 
160
}
 
161
[+ ENDFOR float_type +]
 
162