Mercurial > hg > libxtract

comparison src/vector.c @ 67:2c2ea9326c30

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed bug in xtract_spectrum() FFTW_R2HC now interpreted correctly. Updated ChangeLog and AUTHORS files.
author Jamie Bullock <jamie@postlude.co.uk>
date Wed, 14 Mar 2007 17:20:14 +0000
parents af594496da53
children 9de5628b69a8
comparison
equal deleted inserted replaced
66:af594496da53 67:2c2ea9326c30
53 plan = fftwf_plan_r2r_1d(N, input, rfft, FFTW_R2HC, FFTW_ESTIMATE); 53 plan = fftwf_plan_r2r_1d(N, input, rfft, FFTW_R2HC, FFTW_ESTIMATE);
54 54
55 fftwf_execute(plan); 55 fftwf_execute(plan);
56 56
57 switch(vector){ 57 switch(vector){
58 case XTRACT_MAGNITUDE_SPECTRUM: 58
59 for(n = 0; n < M; n++){ 59 /* case XTRACT_MAGNITUDE_SPECTRUM:
60 result[n] = sqrt(XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) / N; 60 for(n = 1; n < M; n++){
61 result[M + n] = n * q; 61 result[n] = sqrt(XTRACT_SQ(rfft[n]) +
62 } 62 XTRACT_SQ(rfft[N - n - 1])) / N;
63 break; 63 result[M + n] = n * q;
64 }
65 break;
66 */
64 case XTRACT_LOG_MAGNITUDE_SPECTRUM: 67 case XTRACT_LOG_MAGNITUDE_SPECTRUM:
65 for(n = 0; n < M; n++){ 68 for(n = 1; n < M; n++){
66 if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) > XTRACT_LOG_LIMIT) 69 if ((temp = XTRACT_SQ(rfft[n]) +
70 XTRACT_SQ(rfft[N - n - 1])) > XTRACT_LOG_LIMIT)
67 temp = log(sqrt(temp) / N); 71 temp = log(sqrt(temp) / N);
68 else 72 else
69 temp = XTRACT_LOG_LIMIT_DB; 73 temp = XTRACT_LOG_LIMIT_DB;
70 /*Normalise*/ 74 /*Normalise*/
71 result[n] = (temp + XTRACT_DB_SCALE_OFFSET) / XTRACT_DB_SCALE_OFFSET; 75 result[n] =
72 result[M + n] = n * q; 76 (temp + XTRACT_DB_SCALE_OFFSET) / XTRACT_DB_SCALE_OFFSET;
73 } 77 result[M + n] = n * q;
74 break; 78 }
79 break;
80
75 case XTRACT_POWER_SPECTRUM: 81 case XTRACT_POWER_SPECTRUM:
76 for(n = 0; n < M; n++){ 82 for(n = 1; n < M; n++){
77 result[n] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) / NxN; 83 result[n] = (XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n - 1]))
78 result[M + n] = n * q; 84 / NxN;
79 } 85 result[M + n] = n * q;
80 break; 86 }
87 break;
88
81 case XTRACT_LOG_POWER_SPECTRUM: 89 case XTRACT_LOG_POWER_SPECTRUM:
82 for(n = 0; n < M; n++){ 90 for(n = 1; n < M; n++){
83 if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) > XTRACT_LOG_LIMIT) 91 if ((temp = XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n - 1])) >
92 XTRACT_LOG_LIMIT)
84 temp = log(temp / NxN); 93 temp = log(temp / NxN);
85 else 94 else
86 temp = XTRACT_LOG_LIMIT_DB; 95 temp = XTRACT_LOG_LIMIT_DB;
87 result[n] = (temp + XTRACT_DB_SCALE_OFFSET) / XTRACT_DB_SCALE_OFFSET; 96 result[n] = (temp + XTRACT_DB_SCALE_OFFSET) /
88 result[M + n] = n * q; 97 XTRACT_DB_SCALE_OFFSET;
89 } 98 result[M + n] = n * q;
90 break; 99 }
100 break;
101
91 default: 102 default:
92 /* MAGNITUDE_SPECTRUM */ 103 /* MAGNITUDE_SPECTRUM */
93 for(n = 0; n < M; n++){ 104 for(n = 1; n < M; n++){
94 result[n] = sqrt(XTRACT_SQ(rfft[n]) + XTRACT_SQ(rfft[N - n])) / N; 105 result[n] = sqrt(XTRACT_SQ(rfft[n]) +
95 result[M + n] = n * q; 106 XTRACT_SQ(rfft[N - n - 1])) / N;
96 } 107 result[M + n] = n * q;
97 break; 108 }
98 } 109 break;
99 110 }
100 /* result[0] = fabs(temp[0]) / N */ 111
112 /* Set the DC component to 0 */
113 result[0] = result[M] = 0.f;
114 /* Set the Nyquist */
101 result[N] = q * M; 115 result[N] = q * M;
102 116
103 fftwf_destroy_plan(plan); 117 fftwf_destroy_plan(plan);
104 fftwf_free(rfft); 118 fftwf_free(rfft);
105 free(input); 119 free(input);
180 return XTRACT_SUCCESS; 194 return XTRACT_SUCCESS;
181 } 195 }
182 196
183 #else 197 #else
184 198
185 int xtract_magnitude_spectrum(const float *data, const int N, const void *argv, float *result){ 199 int xtract_spectrum(const float *data, const int N, const void *argv, float *result){
186 200
187 XTRACT_NEEDS_FFTW; 201 XTRACT_NEEDS_FFTW;
188 return XTRACT_NO_RESULT; 202 return XTRACT_NO_RESULT;
189 203
190 } 204 }