Makefile

# this is <Makefile>
# ----------------------------------------------------------------------------
# 
# Copyright (c) 2019 by Thomas Forbriger (BFO Schiltach) 
# 
# compile croposp
# 
# ----
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version. 
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# ----
#
# REVISIONS and CHANGES 
#    06/02/2019   V1.0   Thomas Forbriger
# 
# ============================================================================
#

alldata: cpsd1.bin cpsd2.bin cpsd3.bin cpsd4.bin

.PHONY: clean
clean: ; 
	-find . -name \*.bak | xargs --no-run-if-empty /bin/rm -v
	-/bin/rm -vf flist *.bin *.o *.xxx *.ps *.pdf *~ 
	-/bin/rm -vf *.asc

# ======================================================================
# create time series
# ==================
#
# set time series parameters
# --------------------------
CPSD_DT=0.01
CPSD_DUR=8640.
#CPSD_DUR=864.

# ----------------------------------------------------------------------
# create raw time series
# ----------------------
# Several time series of normally distributed random samples with zero mean
# are computed with siggenx. To make them incoherent, each run requires its
# own specific seed value for the random number generator. We use the
# nanosecond field of the system clock.
#
CPSD_SIGGENOPT=-v -o -d $(CPSD_DT) -T $(CPSD_DUR) -ot bin

cpsd1raw.bin:
	/bin/bash -c 'export GSL_RNG_SEED=$$(date +"%N"); \
	  siggenx 14 $@ $(CPSD_SIGGENOPT); \
	  echo seed: $$GSL_RNG_SEED'
cpsd2raw.bin:
	/bin/bash -c 'export GSL_RNG_SEED=$$(date +"%N"); \
	  siggenx 14 $@ $(CPSD_SIGGENOPT); \
	  echo seed: $$GSL_RNG_SEED'
cpsd3raw.bin:
	/bin/bash -c 'export GSL_RNG_SEED=$$(date +"%N"); \
	  siggenx 14 $@ $(CPSD_SIGGENOPT); \
	  echo seed: $$GSL_RNG_SEED'
cpsd4raw.bin:
	/bin/bash -c 'export GSL_RNG_SEED=$$(date +"%N"); \
	  siggenx 14 $@ $(CPSD_SIGGENOPT); \
	  echo seed: $$GSL_RNG_SEED'
# create time series for common noise signal
cpsdbgraw.bin:
	/bin/bash -c 'export GSL_RNG_SEED=$$(date +"%N"); \
	  siggenx 14 $@ $(CPSD_SIGGENOPT); \
	  echo seed: $$GSL_RNG_SEED'

# provide a time series plot to confirm incoherency
rawsig.ps: cpsd1raw.bin cpsd2raw.bin cpsd3raw.bin cpsd4raw.bin cpsdbgraw.bin
	stuplox -d $@/cps -s x -i -ty bin $^

# ----------------------------------------------------------------------
# scale and filter time series
# ----------------------------
#
#  The raw random sequences are adjusted to represent four sets of incoherent
#  noise and one random sequence ('bg') which will be common in all final data
#  files. We understand the first four as a simulation of instrumental noise
#  and the bg-sequence as a representation of background ground noise to be
#  common in all reacordings.
#
# set amplitudes
CPSD_AMP1=0.2
CPSD_AMP2=4.
CPSD_AMP3=100.
CPSD_AMP4=1.e-2
CPSD_AMPBG=10.

# modify white noise to channel specific frequency characteristic
cpsd1mod.bin: cpsd1raw.bin
	printf "fac  $(CPSD_AMP1)\nhpb  1.,2\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
cpsd2mod.bin: cpsd2raw.bin
	printf "fac  $(CPSD_AMP2)\nhpb 4.,4\nlpb  2.,4\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
cpsd3mod.bin: cpsd3raw.bin
	printf "fac  $(CPSD_AMP3)\nlp2  100.,3\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
cpsd4mod.bin: cpsd4raw.bin
	printf "fac  $(CPSD_AMP4)\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
cpsdbgmod.bin: cpsdbgraw.bin
	printf "fac  $(CPSD_AMPBG)\nlp1  2.\navg  0\nhpb  100.,4\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<

# add common noise signal to channel specific noise floor
cpsd%addraw.bin: cpsd%mod.bin cpsdbgmod.bin
	teseco -v -o -a -type bin -Type bin $@ $^

# modify sensor response
cpsd1add.bin: cpsd1addraw.bin
	printf "fac  1.\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
cpsd2add.bin: cpsd2addraw.bin
	printf "fac  1.\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
cpsd3add.bin: cpsd3addraw.bin
	printf "fac  1.\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
cpsd4add.bin: cpsd4addraw.bin
	printf "lpb  10.,2\nend\n" \
	  | tidofi -type bin -Type bin -v -o -cs $@ $<
	
# ----------------------------------------------------------------------
# adjust time series headers
# --------------------------
#
#  recorded data is expected to be taken with slightly different sampling;
#  croposp is expected to be able to handle this
#
# set header fields to reasonable time an channel identifiers
CPSD_DATE=2018/12/16_14:24:16
CPSD_DATE1=$(CPSD_DATE).010980
CPSD_DATE2=$(CPSD_DATE).011050
CPSD_DATE3=$(CPSD_DATE).010760
CPSD_DATE4=$(CPSD_DATE).010760

CPSD_HDOPTS=-verbose -overwrite -itype bin -otype bin -ss TST -sa XX -sc BHZ
cpsd1hd.bin: cpsd1add.bin
	sehefixx $< $@ $(CPSD_HDOPTS) -st $(CPSD_DATE1) -si SYN1
cpsd2hd.bin: cpsd2add.bin
	sehefixx $< $@ $(CPSD_HDOPTS) -st $(CPSD_DATE2) -si SYN2
cpsd3hd.bin: cpsd3add.bin
	sehefixx $< $@ $(CPSD_HDOPTS) -st $(CPSD_DATE3) -si SYN3
cpsd4hd.bin: cpsd4add.bin
	sehefixx $< $@ $(CPSD_HDOPTS) -st $(CPSD_DATE4) -si SYN4

# extract inconsistent time windows
CPSD_TOPT1=-samplesf 3 -samplesl 7
CPSD_TOPT2=-samplesf 25
CPSD_TOPT3=-samplesl 81
CPSD_TOPT4=-samplesl 1
cpsd1.bin: cpsd1hd.bin
	tijerasxx -verbose -overwrite -iformat bin -oformat bin \
	  $(CPSD_TOPT1) $@ $<
cpsd2.bin: cpsd2hd.bin
	tijerasxx -verbose -overwrite -iformat bin -oformat bin \
	  $(CPSD_TOPT2) $@ $<
cpsd3.bin: cpsd3hd.bin
	tijerasxx -verbose -overwrite -iformat bin -oformat bin \
	  $(CPSD_TOPT3) $@ $<
cpsd4.bin: cpsd4hd.bin
	tijerasxx -verbose -overwrite -iformat bin -oformat bin \
	  $(CPSD_TOPT4) $@ $<

# ======================================================================
# analyze time series
# ===================
#
# set croposp options
CPSD_NSEGMENTS=50
CPSD_OVERLAP=0.5
CPSD_DIVISOR=100
CPSD_PADFACTOR=1
CPSD_NLOG=20
CPSD_LOG=--log $(CPSD_NLOG) --avgout
CPSD_OPT=--DEBUG
CPSD_OPT=
OUTPUTS=psd npsd coherency transfer phase

# define string summarizing essential croposp options
# this will be included in the title of diagrams
CPSD_CROPOSP_OPT=--nsegments=$(CPSD_NSEGMENTS) \
		     --divisor=$(CPSD_DIVISOR) \
		     --padfactor=$(CPSD_PADFACTOR) \
		     --overlap=$(CPSD_OVERLAP) $(CPSD_LOG)

# select sets to be analyzed
# four are available, three must be used at least
CPSD_SET=1 2 3
CPSD_SET=1 2 3 4
# compile list of file names as input to croposp
CPSD_BINFILES=$(addprefix cpsd,$(addsuffix .bin,$(CPSD_SET)))
# compile list of raw files (apparent instrumental noise and seismic
# background noise separately)
CPSD_RAWFILES=$(addprefix cpsd,$(addsuffix mod.bin,$(CPSD_SET) bg)) \

# run complete analysis on the apparent seismometer output
$(addsuffix _out.xxx,$(OUTPUTS)): $(CPSD_BINFILES)
	croposp $(CPSD_OPT) --verbose --overwrite --itype bin --trim \
	  --pattern="%I" \
	  $(CPSD_CROPOSP_OPT) \
	  --datetolerance 0.1 \
	  --psd psd_out.xxx \
	  --npsd npsd_out.xxx \
	  --coherency coherency_out.xxx \
	  --transfer transfer_out.xxx \
	  --phase phase_out.xxx \
	  $^

# compute PSD for raw files (apparent instrumental noise and seismic
# background noise separately)
psdraw.xxx: $(CPSD_RAWFILES)
	croposp $(CPSD_OPT) --verbose --overwrite --itype bin --trim \
	  --pattern="%F" \
	  $(CPSD_CROPOSP_OPT) \
	  --datetolerance 0.1 \
	  --psd $@ $^

# compute PSD with foutra, which is used as a reference to confirm correct
# calibration of algorithms
psd%.001.asc: cpsd%.bin
	foutra -v -o -type bin -ASCII=$(patsubst cpsd%.bin,psd%,$<) \
	  -power -rbw=0.1 -demean -detrend -logascii=0.1 -avgascii \
	  -nsegments=20 junk.xxx $^

# convert foutra output to the file format expected by croposplot.py
psd%foutra.xxx: psd%.001.asc
	printf "# #1: foutra (reference): %s\n" $< > $@
	cat $< >> $@

# ======================================================================
# display results
# ===============
#
# location of plot tool
CROPOSPLOT=../croposplot.py

# ----------------------------------------------------------------------
# reference plots
# ---------------
#
#  Display PSD together with results from foutra. 
#
# compile list of file names as analyzed with foutra
CPSD_REFFILES=$(addprefix psd,$(addsuffix foutra.xxx,$(CPSD_SET)))
# compile list of raw files (apparent instrumental noise and seismic
# background noise separately) analyzed with foutra
CPSD_RAWREFFILES=$(addprefix psd,$(addsuffix modfoutra.xxx,$(CPSD_SET) bg))

# display PSD of raw data (apparent instrumental noise and seismic
# background noise separately)
PSDraw.pdf: psdraw.xxx $(CPSD_RAWREFFILES)
	$(CROPOSPLOT) --title="$@ \n $(CPSD_CROPOSP_OPT)" \
	  --ylim=1.e-14:100. \
	  --titfontsize='small' -v --grid -o $@ $^

# display PSD of apparent seismometer output in comparison with foutra
# analysis
PSD.pdf: psd_out.xxx $(CPSD_REFFILES)
	$(CROPOSPLOT) --title="$@ \n $(CPSD_CROPOSP_OPT)" \
	  --titfontsize='small' -v --grid -o $@ $^

# ----------------------------------------------------------------------
# pattern rules for quick plots
# -----------------------------
%.pdf: %_out.xxx
	$(CROPOSPLOT) --title='$@ \n $(CPSD_CROPOSP_OPT)' \
	  --titfontsize='small' -v --grid -o $@ $<

%_lin.pdf: %_out.xxx
	$(CROPOSPLOT) --title="$@ \n $(CPSD_CROPOSP_OPT)" \
	  --titfontsize='small' -v --grid --nology -o \
	  $@ $<

# ----------------------------------------------------------------------
#  previewer commands
#  ------------------
%.pdp: %.pdf; evince $<; /bin/rm -fv $<
%.psp: %.ps; gv $<; /bin/rm -fv $<

# ----- END OF Makefile -----