Unified · Split
| @@ -26,10 +26,6 @@ · 25 unmodified lines | |||
| 26 | #include <string.h> | 26 | #include <string.h> |
| 27 | #include <math.h> | 27 | #include <math.h> |
| 28 | 28 | ||
| 29 | #if (HAVE_FFTW3 == 1) | ||
| 30 | |||
| 31 | #include <fftw3.h> | ||
| 32 | |||
| 33 | #define MAX_SPEC_LEN (1<<18) | 29 | #define MAX_SPEC_LEN (1<<18) |
| 34 | #define MAX_PEAKS 10 | 30 | #define MAX_PEAKS 10 |
| 35 | 31 | ||
| @@ -75,109 +71,6 @@ calculate_snr (float *data, int len) · 39 unmodified lines | |||
| 75 | return snr ; | 71 | return snr ; |
| 76 | } /* calculate_snr */ | 72 | } /* calculate_snr */ |
| 77 | 73 | ||
| 78 | |||
| 79 | |||
| 80 | |||
| 81 | /*- | ||
| 82 | |||
| 83 | The analogues to | ||
| 84 | |||
| 85 | rfftw_create_plan | ||
| 86 | rfftw_one | ||
| 87 | |||
| 88 | with | ||
| 89 | |||
| 90 | FFTW_REAL_TO_COMPLEX | ||
| 91 | FFTW_COMPLEX_TO_REAL | ||
| 92 | |||
| 93 | directions are | ||
| 94 | |||
| 95 | fftw_plan_r2r_1d | ||
| 96 | |||
| 97 | with kind | ||
| 98 | |||
| 99 | FFTW_R2HC or | ||
| 100 | FFTW_HC2R | ||
| 101 | |||
| 102 | followed by fftw_execute. The stride etcetera arguments of | ||
| 103 | rfftw are now in fftw_plan_many_r2r. | ||
| 104 | |||
| 105 | Instead of | ||
| 106 | |||
| 107 | rfftwnd_create_plan or | ||
| 108 | rfftw2d_create_plan or | ||
| 109 | rfftw3d_create_plan | ||
| 110 | |||
| 111 | followed by | ||
| 112 | |||
| 113 | rfftwnd_one_real_to_complex or | ||
| 114 | rfftwnd_one_complex_to_real | ||
| 115 | |||
| 116 | you now use | ||
| 117 | |||
| 118 | fftw_plan_dft_r2c or | ||
| 119 | fftw_plan_dft_r2c_2d or | ||
| 120 | fftw_plan_dft_r2c_3d or | ||
| 121 | fftw_plan_dft_c2r or | ||
| 122 | fftw_plan_dft_c2r_2d or | ||
| 123 | fftw_plan_dft_c2r_3d | ||
| 124 | |||
| 125 | respectively, followed by fftw_execute. As usual, the | ||
| 126 | strides etcetera of | ||
| 127 | |||
| 128 | rfftwnd_real_to_complex or | ||
| 129 | rfftwnd_complex_to_real | ||
| 130 | |||
| 131 | are no specified in the advanced planner routines, | ||
| 132 | |||
| 133 | fftw_plan_many_dft_r2c or | ||
| 134 | fftw_plan_many_dft_c2r | ||
| 135 | |||
| 136 | |||
| 137 | -*/ | ||
| 138 | |||
| 139 | |||
| 140 | |||
| 141 | |||
| 142 | static void | ||
| 143 | log_mag_spectrum (double *input, int len, double *magnitude) | ||
| 144 | { fftw_plan plan = NULL ; | ||
| 145 | |||
| 146 | double maxval ; | ||
| 147 | int k ; | ||
| 148 | |||
| 149 | if (input == NULL || magnitude == NULL) | ||
| 150 | return ; | ||
| 151 | |||
| 152 | plan = fftw_plan_r2r_1d (len, input, magnitude, FFTW_R2HC, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT) ; | ||
| 153 | if (plan == NULL) | ||
| 154 | { printf ("%s : line %d : create plan failed.\n", __FILE__, __LINE__) ; | ||
| 155 | exit (1) ; | ||
| 156 | } ; | ||
| 157 | |||
| 158 | fftw_execute (plan) ; | ||
| 159 | |||
| 160 | /* (k < N/2 rounded up) */ | ||
| 161 | maxval = 0.0 ; | ||
| 162 | for (k = 1 ; k < len / 2 ; k++) | ||
| 163 | { magnitude [k] = sqrt (magnitude [k] * magnitude [k] + magnitude [len - k - 1] * magnitude [len - k - 1]) ; | ||
| 164 | maxval = (maxval < magnitude [k]) ? magnitude [k] : maxval ; | ||
| 165 | } ; | ||
| 166 | |||
| 167 | memset (magnitude + len / 2, 0, len / 2 * sizeof (magnitude [0])) ; | ||
| 168 | |||
| 169 | /* Don't care about DC component. Make it zero. */ | ||
| 170 | magnitude [0] = 0.0 ; | ||
| 171 | |||
| 172 | /* log magnitude. */ | ||
| 173 | for (k = 0 ; k < len ; k++) | ||
| 174 | { magnitude [k] = magnitude [k] / maxval ; | ||
| 175 | magnitude [k] = (magnitude [k] < 1e-15) ? -200.0 : 20.0 * log10 (magnitude [k]) ; | ||
| 176 | } ; | ||
| 177 | |||
| 178 | return ; | ||
| 179 | } /* log_mag_spectrum */ | ||
| 180 | |||
| 181 | /*============================================================================== | 74 | /*============================================================================== |
| 182 | ** There is a slight problem with trying to measure SNR with the method used | 75 | ** There is a slight problem with trying to measure SNR with the method used |
| 183 | ** here; the side lobes of the windowed FFT can look like a noise/aliasing peak. | 76 | ** here; the side lobes of the windowed FFT can look like a noise/aliasing peak. |
| @@ -290,17 +183,94 @@ find_snr (const double *magnitude, int len) · 106 unmodified lines | |||
| 290 | return snr ; | 183 | return snr ; |
| 291 | } /* find_snr */ | 184 | } /* find_snr */ |
| 292 | 185 | ||
| 186 | /*============================================================================== | ||
| 187 | */ | ||
| 188 | |||
| 189 | #if (HAVE_FFTW3 == 1) | ||
| 190 | |||
| 191 | #include <fftw3.h> | ||
| 192 | |||
| 193 | static void | ||
| 194 | log_mag_spectrum (double *input, int len, double *magnitude) | ||
| 195 | { fftw_plan plan = NULL ; | ||
| 196 | |||
| 197 | double maxval ; | ||
| 198 | int k ; | ||
| 199 | |||
| 200 | if (input == NULL || magnitude == NULL) | ||
| 201 | return ; | ||
| 202 | |||
| 203 | plan = fftw_plan_r2r_1d (len, input, magnitude, FFTW_R2HC, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT) ; | ||
| 204 | if (plan == NULL) | ||
| 205 | { printf ("%s : line %d : create plan failed.\n", __FILE__, __LINE__) ; | ||
| 206 | exit (1) ; | ||
| 207 | } ; | ||
| 208 | |||
| 209 | fftw_execute (plan) ; | ||
| 210 | |||
| 211 | /* (k < N/2 rounded up) */ | ||
| 212 | maxval = 0.0 ; | ||
| 213 | for (k = 1 ; k < len / 2 ; k++) | ||
| 214 | { magnitude [k] = sqrt (magnitude [k] * magnitude [k] + magnitude [len - k - 1] * magnitude [len - k - 1]) ; | ||
| 215 | maxval = (maxval < magnitude [k]) ? magnitude [k] : maxval ; | ||
| 216 | } ; | ||
| 217 | |||
| 218 | memset (magnitude + len / 2, 0, len / 2 * sizeof (magnitude [0])) ; | ||
| 219 | |||
| 220 | /* Don't care about DC component. Make it zero. */ | ||
| 221 | magnitude [0] = 0.0 ; | ||
| 222 | |||
| 223 | /* log magnitude. */ | ||
| 224 | for (k = 0 ; k < len ; k++) | ||
| 225 | { magnitude [k] = magnitude [k] / maxval ; | ||
| 226 | magnitude [k] = (magnitude [k] < 1e-15) ? -200.0 : 20.0 * log10 (magnitude [k]) ; | ||
| 227 | } ; | ||
| 228 | |||
| 229 | return ; | ||
| 230 | } /* log_mag_spectrum */ | ||
| 231 | |||
| 293 | #else /* (! HAVE_FFTW3) */ | 232 | #else /* (! HAVE_FFTW3) */ |
| 294 | 233 | ||
| 295 | double | 234 | #include "kiss_fft.h" |
| 296 | calculate_snr (float *data, int len) | ||
| 297 | { double snr = 200.0 ; | ||
| 298 | 235 | ||
| 299 | data = data ; | 236 | static void |
| 300 | len = len ; | 237 | log_mag_spectrum (double *input, int len, double *magnitude) |
| 238 | { fftw_plan plan = NULL ; | ||
| 301 | 239 | ||
| 302 | return snr ; | 240 | double maxval ; |
| 303 | } /* calculate_snr */ | 241 | int k ; |
| 242 | |||
| 243 | if (input == NULL || magnitude == NULL) | ||
| 244 | return ; | ||
| 245 | |||
| 246 | plan = fftw_plan_r2r_1d (len, input, magnitude, FFTW_R2HC, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT) ; | ||
| 247 | if (plan == NULL) | ||
| 248 | { printf ("%s : line %d : create plan failed.\n", __FILE__, __LINE__) ; | ||
| 249 | exit (1) ; | ||
| 250 | } ; | ||
| 251 | |||
| 252 | fftw_execute (plan) ; | ||
| 253 | |||
| 254 | /* (k < N/2 rounded up) */ | ||
| 255 | maxval = 0.0 ; | ||
| 256 | for (k = 1 ; k < len / 2 ; k++) | ||
| 257 | { magnitude [k] = sqrt (magnitude [k] * magnitude [k] + magnitude [len - k - 1] * magnitude [len - k - 1]) ; | ||
| 258 | maxval = (maxval < magnitude [k]) ? magnitude [k] : maxval ; | ||
| 259 | } ; | ||
| 260 | |||
| 261 | memset (magnitude + len / 2, 0, len / 2 * sizeof (magnitude [0])) ; | ||
| 262 | |||
| 263 | /* Don't care about DC component. Make it zero. */ | ||
| 264 | magnitude [0] = 0.0 ; | ||
| 265 | |||
| 266 | /* log magnitude. */ | ||
| 267 | for (k = 0 ; k < len ; k++) | ||
| 268 | { magnitude [k] = magnitude [k] / maxval ; | ||
| 269 | magnitude [k] = (magnitude [k] < 1e-15) ? -200.0 : 20.0 * log10 (magnitude [k]) ; | ||
| 270 | } ; | ||
| 271 | |||
| 272 | return ; | ||
| 273 | } /* log_mag_spectrum */ | ||
| 304 | 274 | ||
| 305 | #endif | 275 | #endif |
| 306 | 276 | ||
| @@ -311,4 +281,3 @@ calculate_snr (float *data, int len) · 4 unmodified lines | |||
| 311 | ** | 281 | ** |
| 312 | ** arch-tag: 7ae937c5-53a5-45d8-8aa2-793de9c35f0a | 282 | ** arch-tag: 7ae937c5-53a5-45d8-8aa2-793de9c35f0a |
| 313 | */ | 283 | */ |
| 314 | |||