Commit 80e5eaea authored by thomas.forbriger's avatar thomas.forbriger
Browse files

ts/croposp [FEATURE]: provide PSD scaling test

parent 84f720a6
# this is <Makefile.psd>
# ----------------------------------------------------------------------------
#
# Copyright (c) 2020 by Thomas Forbriger (BFO Schiltach)
#
# test scaling of power spectral density
#
# Usage:
#
# make -f Makefile.psd clean; make -f Makefile.psd psdplot.psp TS_RMS=70.71
# make -f Makefile.psd clean; make -f Makefile.psd psdplot.psp TS_RMS=180.
# make -f Makefile.psd clean; make -f Makefile.psd psdplot.psp TS_RMS=100.
#
# other parameter (see below) may be changed as well
#
# REVISIONS and CHANGES
# 02/10/2020 V1.0 Thomas Forbriger
#
# ============================================================================
#
.PHONY: clean
clean:
/bin/rm -fv *.bin *.ps *.bak *.xxx
# ======================================================================
# create time series
# ==================
#
# Compute a time series of random Gaussian distributed samples (white noise)
# with a given rms amplitude.
#
# We use the libgsl random number generators to produce white noise time
# series, which are as incoherent as possible. The recommended random number
# generator is ranlxd
#
# run
# randomseries -xhelp=gslrng
# to see comments
#
# set time series parameters
# --------------------------
# sampling interval
TS_DT=0.01
# number of samples
TS_N=1000000
# rms amplitude
TS_RMS=100.
#
# Nyquist frequency
FNy=$(shell dc -e '7 k 0.5 $(TS_DT) / p')
# power spectral density
PSD=$(shell dc -e '7 k $(TS_RMS) $(TS_RMS) * 2. / $(FNy) / p')
# one-sided PSD
PSD2=$(shell dc -e '7 k 2. $(PSD) * p')
#
.PHONY: tsparameters
tsparameters:
@LANG=C; printf "sampling interval: %10.5f s\n" $(TS_DT)
@LANG=C; printf "number of samples: %10d\n" $(TS_N)
@LANG=C; printf "rms amplitude: %10.3f counts\n" $(TS_RMS)
@LANG=C; printf "Nyquist frequency: %10.3f Hz\n" $(FNy)
@LANG=C; printf "power spectral density: %10.2f counts**2/Hz\n" $(PSD)
@LANG=C; printf "one-sided PSD: %10.2f counts**2/Hz\n" $(PSD2)
#
# set random number generator type
# --------------------------------
# Choose:
RNGTYPE=ranlxd2
RNGSEED=0
# create to raw time series
whitenoise.bin:
randomseries -v -o -t bin -rngtype $(RNGTYPE) \
-nsamples $(TS_N) \
-dt $(TS_DT) -seed $(RNGSEED) \
-std $(TS_RMS) -mean 0. $@
# ----------------------------------------------------------------------
# compute PSD with croposp
# ------------------------
#
# set croposp options
CROPOSP_NSEGMENTS=20
CROPOSP_OVERLAP=0.5
CROPOSP_DIVISOR=100
CROPOSP_PADFACTOR=1
CROPOSP_AVGOUT=--avgout
CROPOSP_NLOG=20
CROPOSP_LOG=--log $(CROPOSP_NLOG) $(CROPOSP_AVGOUT)
CROPOSP_OPT=--DEBUG
CROPOSP_OPT=
# define string summarizing essential croposp options
# this will be included in the title of diagrams
CROPOSP_OPT=--nsegments=$(CROPOSP_NSEGMENTS) \
--divisor=$(CROPOSP_DIVISOR) \
--padfactor=$(CROPOSP_PADFACTOR) \
--overlap=$(CROPOSP_OVERLAP) \
$(CROPOSP_LOG)
# set standard (non-processing options)
CROPOSP_STDOPT=--verbose --overwrite --itype bin
# compute PSD for raw files (apparent instrumental noise and seismic
# background noise separately)
psdcroposp.xxx: whitenoise.bin
croposp $(CROPOSP_OPT) $(CROPOSP_STDOPT) \
--psd $@ $< $<
# ----------------------------------------------------------------------
# compute PSD with foutra (as a reference)
# ----------------------------------------
#
FOUTRA_OPT=-nsegments=$(CROPOSP_NSEGMENTS)
FOUTRA_STDOPT=-v -o -type bin -power -rbw=0.1 -demean -detrend \
-logascii=0.1 -avgascii
# compute PSD with foutra, which is used as a reference to confirm correct
# calibration of algorithms
psdfoutra.001.asc: whitenoise.bin
foutra $(FOUTRA_STDOPT) $(FOUTRA_OPT) \
-ASCII=psdfoutra junk.xxx $<
# convert foutra output to the file format expected by croposplot.py
psdfoutra.xxx: psdfoutra.001.asc
printf "# #1: foutra (reference): %s\n" $< > $@
cat $< >> $@
# ======================================================================
# plot diagram
# ------------
#
psdplot.xxx: psdfoutra.xxx psdcroposp.xxx
echo '#!/bin/gnuplot -f' > $@
echo 'set terminal postscript color enhanced' >> $@
echo 'set output "psdplot.ps"' >> $@
echo 'set logscale xy' >> $@
echo 'set grid' >> $@
echo 'set dashtype 2 (5,2,2,2)' >> $@
echo 'set style line 1 lt 1 lw 2 lc "#ff0000"' >> $@
echo 'set style line 2 lt 1 lw 2 lc "#0000ff"' >> $@
echo 'set style line 3 lt 3 lw 2 dt 2 lc "#8888ff"' >> $@
echo 'set ylabel "2*P(f) / counts^{2} Hz^{-1}"' >> $@
echo 'set xlabel "frequency / Hz"' >> $@
echo 'set title "PSD of white noise with rms-amplitude of $(TS_RMS) counts"' >> $@
printf 'set label "sampling interval: %10.5f s"\\\n' $(TS_DT) >> $@
printf ' at graph 0.4,0.4 left\n' >> $@
printf 'set label "number of samples: %10d"\\\n' $(TS_N) >> $@
printf ' at graph 0.4,0.35 left\n' >> $@
printf 'set label "rms amplitude: A_{rms}=%10.3f counts"\\\n' $(TS_RMS) >> $@
printf ' at graph 0.4,0.3 left\n' >> $@
printf 'set label "Nyquist frequency: f_{Ny}=%10.3f Hz"\\\n' $(FNy) >> $@
printf ' at graph 0.4,0.25 left\n' >> $@
printf 'set label "power spectral density: P(f)=%10.2f counts**2/Hz"\\\n' $(PSD) >> $@
printf ' at graph 0.4,0.2 left\n' >> $@
printf 'set label "one-sided PSD: 2*P(f)=%10.2f counts**2/Hz"\\\n' $(PSD2) >> $@
printf ' at graph 0.4,0.15 left\n' >> $@
printf 'set label "P(f)=A_{rms}^{2}/(2*f_{Ny}) A_{rms}=sqrt(2*P(f)*f_{Ny})"\\\n' $(PSD2) >> $@
printf ' at graph 0.4,0.1 left\n' >> $@
printf 'plot "$(word 1,$^)" u 1:2 w l ls 1 \\\n' >> $@
printf 't "$(patsubst psd%.xxx,%,$(word 1,$^))" \\\n' >> $@
printf ', \\\n' >> $@
printf '"$(word 2,$^)" u 1:2 w l ls 2\\\n' >> $@
printf 't "$(patsubst psd%.xxx,%,$(word 2,$^))" \\\n' >> $@
printf ', \\\n' >> $@
printf '$(PSD2) ls 3 t "theoretical value"\n' >> $@
psdplot.ps: psdplot.xxx psdfoutra.xxx psdcroposp.xxx
gnuplot $<
# ======================================================================
# previewer commands
# ------------------
%.pdp: %.pdf; evince $<; /bin/rm -fv $<
%.psp: %.ps; gv $<; /bin/rm -fv $<
# ----- END OF Makefile.psd -----
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment