Changeset 40219 for branches/eam_branches/ohana.20170822
- Timestamp:
- Nov 28, 2017, 10:18:12 AM (9 years ago)
- Location:
- branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test
- Files:
-
- 2 edited
-
periodogram-fm.sh (modified) (15 diffs)
-
periodogram.sh (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram-fm.sh
r40146 r40219 55 55 56 56 periodogram_fm t f df 5 50 period power 57 #periodogram t f 5 50 period power58 57 59 58 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f … … 88 87 89 88 periodogram_fm t f df 1 10 period power 90 #periodogram t f 1 10 period power91 89 92 90 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f … … 97 95 if (abs ($peakpos - $P) > 0.05) 98 96 $PASS = 0 97 end 98 99 if ($PLOT) 100 lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line 99 101 end 100 102 end … … 117 119 118 120 periodogram_fm t f df 2 30 period power 119 #periodogram t f 2 30 period power120 121 121 122 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f … … 126 127 if (abs ($peakpos - $P) > 0.05) 127 128 $PASS = 0 129 end 130 131 if ($PLOT) 132 lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line 128 133 end 129 134 end … … 145 150 set df = 0.01 + zero(f) 146 151 147 periodogram_fm t f 2 30 period power152 periodogram_fm t f df 2 30 period power 148 153 149 154 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f … … 155 160 $PASS = 0 156 161 end 157 end 158 159 # test using random samples, offset start, non-zero DC, some noise 162 163 if ($PLOT) 164 lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line 165 end 166 end 167 168 # test using 300 random samples, offset start, non-zero DC, some noise 160 169 macro test6 161 170 $PASS = 1 … … 178 187 set df = 0.01 + zero(f) 179 188 180 periodogram_fm t f 2 30 period power189 periodogram_fm t f df 2 30 period power 181 190 182 191 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f … … 188 197 $PASS = 0 189 198 end 190 end 191 192 # test using fewer random samples, offset start, non-zero DC, some noise 199 200 if ($PLOT) 201 lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line 202 end 203 end 204 205 # test using 100 fewer random samples, offset start, non-zero DC, some noise 193 206 macro test7 194 207 $PASS = 1 … … 211 224 set df = 0.01 + zero(f) 212 225 213 periodogram_fm t f 2 30 period power226 periodogram_fm t f df 2 30 period power 214 227 215 228 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f … … 221 234 $PASS = 0 222 235 end 223 end 224 225 # test using fewer random samples, high frequency, non-zero DC, some noise 236 237 if ($PLOT) 238 lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line 239 end 240 end 241 242 # test using Ndays random samples, RR Lyrae-sized light curves (0.7 mag), 243 # optional noise level 226 244 macro test8 227 if ($0 != 2)228 echo "USAGE: test8: Ndays")245 if ($0 != 4) 246 echo "USAGE: test8: Period Ndays (df)" 229 247 break 230 248 end 231 249 232 250 local Ndays 233 $Ndays = $1 234 235 $PASS = 1 236 break -auto off 237 238 local P PI 239 $PI = 3.14159265359 240 $P = 0.8*rnd(0) + 0.2 251 $P = $1 252 $Ndays = $2 253 $dM = $3 254 255 $PASS = 1 256 break -auto off 257 258 local PI 259 $PI = 3.14159265359 241 260 $trueP = $P 242 261 … … 246 265 247 266 # t is a time in days, but we always have 4 within 1 hour: 248 set t 0 = int(100 * rnd(x))249 set dtx = (3/24) * rnd(x) 250 set t0 = t 0+ dtx267 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100 268 set dtx = (3/24) * rnd(x); # choose a starting time within that night 269 set t0 = tday + dtx 251 270 252 271 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440 … … 260 279 set tmp = t0 + dt3; concat tmp t 261 280 262 set fraw = sin(2*$PI*t/$P) + 0.5281 set fraw = 0.75*sin(2*$PI*t/$P) 263 282 264 283 # 0.05 : peakpos = 14.95 265 284 # 0.10 : peakpos = 15.04 ( 266 gaussdev df t[] 0.0 0.25285 gaussdev df t[] 0.0 $dM 267 286 set f = fraw + df 268 set df = 0.01+ zero(f)269 270 periodogram_fm t f 0.1 2.0 period power287 set df = $dM + zero(f) 288 289 periodogram_fm t f df 0.1 20.0 period power 271 290 272 291 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f … … 278 297 $PASS = 0 279 298 end 280 end 299 if ($PLOT) 300 lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line 301 end 302 end 303 304 # test using Ndays random samples, RR Lyrae-sized light curves (0.7 mag), 305 # optional noise level 306 # compare periodogram and periodogram_fm 307 macro test9 308 if ($0 != 4) 309 echo "USAGE: test8: Period Ndays (df)" 310 break 311 end 312 313 local Ndays 314 $P = $1 315 $Ndays = $2 316 $dM = $3 317 318 $PASS = 1 319 break -auto off 320 321 local PI 322 $PI = 3.14159265359 323 $trueP = $P 324 325 delete -q x t f period power 326 327 create x 0 $Ndays 328 329 # t is a time in days, but we always have 4 within 1 hour: 330 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100 331 set dtx = (3/24) * rnd(x); # choose a starting time within that night 332 set t0 = tday + dtx 333 334 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440 335 set dt2 = (15.0 / 1440) * rnd(x) + (15 + 7.5) / 1440 336 set dt3 = (15.0 / 1440) * rnd(x) + (30 + 7.5) / 1440 337 338 delete -q t 339 concat t0 t 340 set tmp = t0 + dt1; concat tmp t 341 set tmp = t0 + dt2; concat tmp t 342 set tmp = t0 + dt3; concat tmp t 343 344 set fraw = 0.75*sin(2*$PI*t/$P) 345 346 # 0.05 : peakpos = 14.95 347 # 0.10 : peakpos = 15.04 ( 348 gaussdev df t[] 0.0 $dM 349 set f = fraw + df 350 set df = $dM + zero(f) 351 352 periodogram_fm t f df 0.1 20.0 period_fm power_fm 353 periodogram t f 0.1 20.0 period power 354 355 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f 356 # lim -n 1 period power; clear; box; plot period power 357 358 peak -q period_fm power_fm 359 $peakval_fm = $peakval 360 361 peak -q period power 362 vstat -q power 363 set power = power / $MAX 364 365 # if (abs ($peakpos - $P) > 0.05) 366 # $PASS = 0 367 # end 368 369 set freq = 1 / period 370 set freq_fm = 1 / period_fm 371 $Freq = 1 / $P 372 373 if ($PLOT) 374 if (1) 375 lim period power; clear; box 376 line -c red70 -lw 3 $P 0 to $P $peakval_fm; 377 plot period power -x line -c grey70 -lw 2 378 plot period_fm power_fm -x line -c black 379 else 380 lim freq power; clear; box 381 line -c red70 -lw 3 $Freq 0 to $Freq 1.0 382 plot freq power -x line -c grey70 -lw 2 383 plot freq_fm power_fm -x line -c black 384 end 385 end 386 end 387 388 389 # test using Ndays random samples, RR Lyrae-sized light curves (0.7 mag), 390 # optional noise level 391 # compare periodogram and periodogram_fm 392 macro test10 393 if ($0 != 4) 394 echo "USAGE: test8: Period Ndays (df)" 395 break 396 end 397 398 local Ndays 399 $P = $1 400 $Ndays = $2 401 $dM = $3 402 403 $PASS = 1 404 break -auto off 405 406 local PI 407 $PI = 3.14159265359 408 $trueP = $P 409 410 delete -q x t f period power 411 412 create x 0 $Ndays 413 414 # t is a time in days, but we always have 4 within 1 hour: 415 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100 416 set dtx = (3/24) * rnd(x); # choose a starting time within that night 417 set t0 = tday + dtx 418 419 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440 420 set dt2 = (15.0 / 1440) * rnd(x) + (15 + 7.5) / 1440 421 set dt3 = (15.0 / 1440) * rnd(x) + (30 + 7.5) / 1440 422 423 delete -q t 424 concat t0 t 425 set tmp = t0 + dt1; concat tmp t 426 set tmp = t0 + dt2; concat tmp t 427 set tmp = t0 + dt3; concat tmp t 428 429 set fraw = 0.75*sin(2*$PI*t/$P) 430 431 # 0.05 : peakpos = 14.95 432 # 0.10 : peakpos = 15.04 ( 433 gaussdev df t[] 0.0 $dM 434 set f = fraw + df 435 set df = $dM + zero(f) 436 437 periodogram_fm t f df 0.05 20.0 period_fm power_fm 438 439 gaussdev df t[] 0.0 $dM 440 set Fo = df 441 periodogram_fm t Fo df 0.05 20.0 period power 442 443 # lim -n 0 t f; clear; box; plot -x 2 -pt 2 t f 444 # lim -n 1 period power; clear; box; plot period power 445 446 peak -q period_fm power_fm 447 $peakval_fm = $peakval 448 449 peak -q period power 450 451 # if (abs ($peakpos - $P) > 0.05) 452 # $PASS = 0 453 # end 454 455 set freq = 1 / period 456 set freq_fm = 1 / period_fm 457 $Freq = 1 / $P 458 459 if ($PLOT) 460 if (1) 461 lim period_fm power_fm; clear; box 462 line -c red70 -lw 3 $P 0 to $P $peakval_fm; 463 plot period power -x line -c grey70 -lw 2 464 plot period_fm power_fm -x line -c black 465 else 466 lim freq power; clear; box 467 line -c red70 -lw 3 $Freq 0 to $Freq 1.0 468 plot freq power -x line -c grey70 -lw 2 469 plot freq_fm power_fm -x line -c black 470 end 471 end 472 end 473 474 # we have time (MJD) and mag 475 # we generate the folded lightcure and measure sigma relative to the smoothed version (bins of 0.1 period) 476 macro fold.one.period 477 if ($0 != 5) 478 echo "USAGE: fold.one.period (time) (mag) (magErr) (period)" 479 break 480 end 481 482 local myTime myMag myMagErr myPeriod 483 $myTime = $1 484 $myMag = $2 485 $myMagErr = $3 486 $myPeriod = $4 487 488 set phi = $myTime / $myPeriod - int($myTime / $myPeriod) 489 490 if ($PLOT_FOLD) 491 lim -n phi phi $myMag; clear; box; 492 end 493 494 delete -q magResid 495 496 $dPhi = 0.05; # half of bin size 497 create nphi $dPhi {1 + $dPhi} {2*$dPhi} 498 set magR = zero(nphi) 499 set magS = zero(nphi) 500 for i 0 nphi[] 501 subset tmp_mag_sub = $myMag where (phi >= nphi[$i] - $dPhi) && (phi < nphi[$i] + $dPhi) 502 vstat -q tmp_mag_sub 503 magR[$i] = $MEDIAN 504 magS[$i] = $SIGMA 505 506 set magDelta = tmp_mag_sub - $MEDIAN 507 concat magDelta magResid 508 509 if ($PLOT_FOLD) 510 subset tmp_phi_sub = phi where (phi >= nphi[$i] - $dPhi) && (phi < nphi[$i] + $dPhi) 511 if ($i % 2) 512 plot tmp_phi_sub tmp_mag_sub -pt 7 -sz 3 -c blue -lw 2 513 else 514 plot tmp_phi_sub tmp_mag_sub -pt 7 -sz 3 -c red -lw 2 515 end 516 end 517 end 518 519 if ($PLOT_FOLD) 520 plot -pt 10 -sz 1.5 phi $myMag -dy $myMagErr 521 plot -pt 2 -sz 2.0 -c red nphi magR -dy magS 522 end 523 524 vstat -q magResid 525 end 526 527 281 528 282 529 # Memory test … … 298 545 299 546 for i 0 100 300 periodogram_fm t f 2 30 period power547 periodogram_fm t f df 2 30 period power 301 548 end 302 549 -
branches/eam_branches/ohana.20170822/src/opihi/cmd.data/test/periodogram.sh
r40146 r40219 209 209 # test using fewer random samples, high frequency, non-zero DC, some noise 210 210 macro test8 211 if ($0 != 2)212 echo "USAGE: test8: Ndays")211 if ($0 != 4) 212 echo "USAGE: test8: Period Ndays df" 213 213 break 214 214 end 215 215 216 216 local Ndays 217 $Ndays = $1 218 219 $PASS = 1 220 break -auto off 221 222 local P PI 223 $PI = 3.14159265359 224 $P = 0.8*rnd(0) + 0.2 217 $P = $1 218 $Ndays = $2 219 $dM = $3 220 221 $PASS = 1 222 break -auto off 223 224 local PI 225 $PI = 3.14159265359 225 226 $trueP = $P 226 227 … … 230 231 231 232 # t is a time in days, but we always have 4 within 1 hour: 232 set t 0 = int(100 * rnd(x))233 set dtx = (3/24) * rnd(x) 234 set t0 = t 0+ dtx233 set tday = int(100 * rnd(x)); # choose Ndays random days between 0 and 100 234 set dtx = (3/24) * rnd(x); # choose a starting time within that night 235 set t0 = tday + dtx 235 236 236 237 set dt1 = (15.0 / 1440) * rnd(x) + ( 0 + 7.5) / 1440 … … 244 245 set tmp = t0 + dt3; concat tmp t 245 246 246 set fraw = sin(2*$PI*t/$P) + 0.5247 set fraw = 0.75*sin(2*$PI*t/$P) 247 248 248 249 # 0.05 : peakpos = 14.95 249 250 # 0.10 : peakpos = 15.04 ( 250 gaussdev df t[] 0.0 0.25251 gaussdev df t[] 0.0 $dM 251 252 set f = fraw + df 252 253 … … 260 261 if (abs ($peakpos - $P) > 0.05) 261 262 $PASS = 0 263 end 264 if ($PLOT) 265 lim period power; clear; box; line -c red70 -lw 3 $P 0 to $P $peakval; plot period power -x line 262 266 end 263 267 end
Note:
See TracChangeset
for help on using the changeset viewer.
