Mercurial > hg > libxtract

comparison src/vector.c @ 70:adcecb0b5d99

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Further updated xtract_spectrum() to hopefully fix fft iteration bug and nyquist/DC inclusion. Added new boolean argument 'withDC' to select whether the DC component is required in the output
author Jamie Bullock <jamie@postlude.co.uk>
date Mon, 19 Mar 2007 18:58:21 +0000
parents 99ea1aae68ec
children 06ee24d94059
comparison
equal deleted inserted replaced
69:99ea1aae68ec 70:adcecb0b5d99
33 33
34 int xtract_spectrum(const float *data, const int N, const void *argv, float *result){ 34 int xtract_spectrum(const float *data, const int N, const void *argv, float *result){
35 35
36 float *input, *rfft, q, temp; 36 float *input, *rfft, q, temp;
37 size_t bytes; 37 size_t bytes;
38 int n , NxN, M, vector; 38 int n , NxN, M, vector, withDC;
39 fftwf_plan plan; 39 fftwf_plan plan;
40 40
41 M = N >> 1; 41 M = N >> 1;
42 NxN = XTRACT_SQ(N); 42 NxN = XTRACT_SQ(N);
43 withDC = 0;
43 44
44 rfft = (float *)fftwf_malloc(N * sizeof(float)); 45 rfft = (float *)fftwf_malloc(N * sizeof(float));
45 input = (float *)malloc(bytes = N * sizeof(float)); 46 input = (float *)malloc(bytes = N * sizeof(float));
46 input = memcpy(input, data, bytes); 47 input = memcpy(input, data, bytes);
47 48
48 q = *(float *)argv; 49 q = *(float *)argv;
49 vector = (int)*((float *)argv+1); 50 vector = (int)*((float *)argv+1);
51 withDC = (int)*((float *)argv+2);
50 52
51 XTRACT_CHECK_q; 53 XTRACT_CHECK_q;
52 54
53 plan = fftwf_plan_r2r_1d(N, input, rfft, FFTW_R2HC, FFTW_ESTIMATE); 55 plan = fftwf_plan_r2r_1d(N, input, rfft, FFTW_R2HC, FFTW_ESTIMATE);
54 56
55 fftwf_execute(plan); 57 fftwf_execute(plan);
56 58
57 switch(vector){ 59 switch(vector){
58 60
59 /* case XTRACT_MAGNITUDE_SPECTRUM:
60 for(n = 1; n < M; n++){
61 result[n] = sqrt(XTRACT_SQ(rfft[n]) +
62 XTRACT_SQ(rfft[N - n - 1])) / N;
63 result[M + n] = n * q;
64 }
65 break;
66 */
67 case XTRACT_LOG_MAGNITUDE_SPECTRUM: 61 case XTRACT_LOG_MAGNITUDE_SPECTRUM:
68 for(n = 1; n < M; n++){ 62 for(n = 1; n < M; n++){
69 if ((temp = XTRACT_SQ(rfft[n]) + 63 if ((temp = XTRACT_SQ(rfft[n]) +
70 XTRACT_SQ(rfft[N - n - 1])) > XTRACT_LOG_LIMIT) 64 XTRACT_SQ(rfft[N - n])) > XTRACT_LOG_LIMIT)
71 temp = log(sqrt(temp) / N); 65 temp = log(sqrt(temp) / N);
72 else 66 else
73 temp = XTRACT_LOG_LIMIT_DB; 67 temp = XTRACT_LOG_LIMIT_DB;
74 /*Normalise*/ 68 if(withDC) {
75 result[n] = 69 result[n] =
76 (temp + XTRACT_DB_SCALE_OFFSET) / XTRACT_DB_SCALE_OFFSET; 70 /*Normalise*/
77 result[M + n] = n * q; 71 (temp + XTRACT_DB_SCALE_OFFSET) /
72 XTRACT_DB_SCALE_OFFSET;
73 result[M + n + 1] = n * q;
74 }
75 else {
76 result[n - 1] =
77 (temp + XTRACT_DB_SCALE_OFFSET) /
78 XTRACT_DB_SCALE_OFFSET;
79 result[M + n - 1] = n * q;
80 }
78 } 81 }
79 break; 82 break;
80 83
81 case XTRACT_POWER_SPECTRUM: 84 case XTRACT_POWER_SPECTRUM:
82 for(n = 1; n < M; n++){ 85 for(n = 1; n < M; n++){
83 result[n] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n - 1])) 86 if(withDC){
84 / NxN; 87 result[n] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n]))
85 result[M + n] = n * q; 88 / NxN;
89 result[M + n + 1] = n * q;
90 }
91 else {
92 result[n - 1] =
93 (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) / NxN;
94 result[M + n - 1] = n * q;
95 }
86 } 96 }
87 break; 97 break;
88 98
89 case XTRACT_LOG_POWER_SPECTRUM: 99 case XTRACT_LOG_POWER_SPECTRUM:
90 for(n = 1; n < M; n++){ 100 for(n = 1; n < M; n++){
91 if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n - 1])) > 101 if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) >
92 XTRACT_LOG_LIMIT) 102 XTRACT_LOG_LIMIT)
93 temp = log(temp / NxN); 103 temp = log(temp / NxN);
94 else 104 else
95 temp = XTRACT_LOG_LIMIT_DB; 105 temp = XTRACT_LOG_LIMIT_DB;
96 result[n] = (temp + XTRACT_DB_SCALE_OFFSET) / 106 if(withDC){
97 XTRACT_DB_SCALE_OFFSET; 107 result[n] = (temp + XTRACT_DB_SCALE_OFFSET) /
98 result[M + n] = n * q; 108 XTRACT_DB_SCALE_OFFSET;
109 result[M + n + 1] = n * q;
110 }
111 else {
112 result[n - 1] = (temp + XTRACT_DB_SCALE_OFFSET) /
113 XTRACT_DB_SCALE_OFFSET;
114 result[M + n - 1] = n * q;
115 }
99 } 116 }
100 break; 117 break;
101 118
102 default: 119 default:
103 /* MAGNITUDE_SPECTRUM */ 120 /* MAGNITUDE_SPECTRUM */
104 for(n = 1; n < M; n++){ 121 for(n = 1; n < M; n++){
105 result[n] = sqrt(XTRACT_SQ(rfft[n]) + 122 if(withDC){
106 XTRACT_SQ(rfft[N - n - 1])) / N; 123 result[n] = sqrt(XTRACT_SQ(rfft[n]) +
107 result[M + n] = n * q; 124 XTRACT_SQ(rfft[N - n])) / N;
125 result[M + n + 1] = n * q;
126 }
127 else {
128 result[n - 1] = sqrt(XTRACT_SQ(rfft[n]) +
129 XTRACT_SQ(rfft[N - n])) / N;
130 result[M + n - 1] = n * q;
131 }
108 } 132 }
109 break; 133 break;
110 } 134 }
111 135
112 /* Set the DC component to 0 */ 136 if(withDC){
113 result[0] = result[M] = 0.f; 137 /* The DC component */
114 /* Set the Nyquist */ 138 result[0] = XTRACT_SQ(rfft[0]);
115 result[N] = q * M; 139 result[M + 1] = 0.f;
140 /* The Nyquist */
141 result[M] = XTRACT_SQ(rfft[M]);
142 result[N + 1] = q * M;
143 }
144 else {
145 /* The Nyquist */
146 result[M - 1] = XTRACT_SQ(rfft[M]);
147 result[N - 1] = q * M;
148 }
116 149
117 fftwf_destroy_plan(plan); 150 fftwf_destroy_plan(plan);
118 fftwf_free(rfft); 151 fftwf_free(rfft);
119 free(input); 152 free(input);
120 153