Commit 0a26fcba authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

proceeding

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 2384
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 5ec0c0aa
# this is <Makefile>
# ----------------------------------------------------------------------------
# $Id: Makefile,v 1.41 2007-09-12 07:43:04 tforb Exp $
# $Id: Makefile,v 1.42 2007-09-12 15:11:14 tforb Exp $
#
# Copyright (c) 2005 by Thomas Forbriger (BFO Schiltach)
#
......@@ -67,7 +67,8 @@ susei evelo tesiff teswf: %: %.o
newprog $@
phasedsignals hamres siggen smoos dise: %: %.o
$(CC) $(CFLAGS) $< -o $@ -ltf -lsffu -ltime -lsff $(F2CLIB) \
$(CC) $(CFLAGS) $< -o $@ -ltf -lsffu -ltime -lsff -lgsl -lm \
-lgslcblas $(F2CLIB) \
-L$(LOCLIBDIR)
newprog $@
......
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: foutra.cc,v 1.6 2007-09-12 10:58:08 tforb Exp $
* $Id: foutra.cc,v 1.7 2007-09-12 15:11:14 tforb Exp $
* \author Thomas Forbriger
* \date 25/07/2006
*
......@@ -20,13 +20,16 @@
#define FOUTRA_VERSION \
"FOUTRA V1.2 Fourier transforms"
#define FOUTRA_CVSID \
"$Id: foutra.cc,v 1.6 2007-09-12 10:58:08 tforb Exp $"
"$Id: foutra.cc,v 1.7 2007-09-12 15:11:14 tforb Exp $"
#include <iostream>
#include <tfxx/commandline.h>
#include <aff/series.h>
#include <aff/iterator.h>
#include <aff/dump.h>
#include <aff/seriesoperators.h>
#include <aff/functions/avg.h>
#include <aff/subarray.h>
#include <tsxx/tsxx.h>
#include <tsxx/tapers.h>
#include <tsxx/sffheaders.h>
......@@ -48,6 +51,9 @@ struct Options {
bool verbose, overwrite, debug;
std::string inputformat;
bool amplitudespectrum, powerspectrum, boxcartaper;
bool avgconstbw, avgrelbw;
double decades;
int avgsamples;
}; // struct Options
// values type to be used for samples
......@@ -75,7 +81,7 @@ int main(int iargc, char* argv[])
FOUTRA_VERSION "\n"
"usage: foutra [-v] [-o] [-type type] [-D]" "\n"
" [-amplitude] [-power]" "\n"
" [-boxcar]" "\n"
" [-boxcar] [-avg[=n]] [-rbw[=n]]" "\n"
" outfile infile [t:T] [infile [t:T] ...]" "\n"
" or: foutra --help|-h" "\n"
};
......@@ -96,6 +102,9 @@ int main(int iargc, char* argv[])
"-amplitude calculate amplitude spectrum" "\n"
"-power calculate power spectrum" "\n"
"-type type select input file type" "\n"
"-type type select input file type" "\n"
"-avg[=n] smooth power spectrum by averaging over n samples" "\n"
"-rbw[=n] smooth power spectrum by averaging over n decades" "\n"
};
// define commandline options
......@@ -118,6 +127,10 @@ int main(int iargc, char* argv[])
{"power",arg_no,"-"},
// 7: apply boxcar taper
{"boxcar",arg_no,"-"},
// 8: average over constant number of samples
{"avg",arg_opt,"21"},
// 9: average over constant relative bandwidth (in decades)
{"rbw",arg_opt,"0.167"},
{NULL}
};
......@@ -154,6 +167,10 @@ int main(int iargc, char* argv[])
opt.amplitudespectrum=cmdline.optset(5);
opt.powerspectrum=cmdline.optset(6);
opt.boxcartaper=cmdline.optset(7);
opt.avgconstbw=cmdline.optset(8);
opt.avgsamples=cmdline.int_arg(8);
opt.avgrelbw=cmdline.optset(9);
opt.decades=cmdline.double_arg(9);
if (opt.verbose)
{ cout << FOUTRA_VERSION << endl << FOUTRA_CVSID << endl; }
......@@ -254,6 +271,7 @@ int main(int iargc, char* argv[])
// apply taper
if (!opt.boxcartaper) { taper.apply(series); }
// call FFT for amplitude or power spectrum
if (opt.amplitudespectrum || opt.powerspectrum)
{
TFXX_debug(opt.debug, "main", "call FFT");
......@@ -277,6 +295,8 @@ int main(int iargc, char* argv[])
DUMP(coeff);
}
}
// take sqrt for amplitude spectrum
if (opt.amplitudespectrum)
{
aff::Iterator<Tseries> S(series);
......@@ -287,6 +307,98 @@ int main(int iargc, char* argv[])
}
}
// apply appropriate scaling and averaging for power spectrum
if (opt.powerspectrum)
{
// length of time series
double T=wid2.dt*wid2.nsamples;
// frequency sampling interval
double df=1/T;
// scaling factor to adjust for taper effect
double tapscaling;
if (opt.boxcartaper)
{
tapscaling=1;
}
else
{
// scaling factor for Hanning
// see below for derivation of this factor
tapscaling=1.6976527;
}
// we have an energy spectrum so far
// adjust scaling factor to obtain signal power
double scalingfactor=tapscaling/T;
series *= scalingfactor;
// apply smoothing
if (opt.avgconstbw)
{
if (opt.verbose)
{
cout
<< "apply smoothing by averaging over "
<< opt.avgsamples
<< " samples"
<< endl;
}
// create a copy
Tseries p(series.size());
p.copyin(series);
int d=opt.avgsamples/2+1;
for (int k=p.f(); k<=p.l(); ++k)
{
int k1=k-d;
int k2=k+d;
k1 = k1 < p.f() ? p.f() : k1;
k2 = k2 > p.l() ? p.l() : k2;
// cout << k1 << " " << k << " " << k2 << endl;
TFXX_assert(k2>k1, "negative bandwidth for averaging")
Tseries sub=aff::subarray(p)(k1,k2);
series(k)=aff::func::avg(sub);
series(k)=0;
for (int j=sub.f(); j<=sub.l(); ++j)
{
series(k) += sub(j);
}
series(k) /= double(sub.size());
}
}
else if (opt.avgrelbw)
{
if (opt.verbose)
{
cout
<< "apply smoothing by averaging over "
<< opt.decades
<< " decades"
<< endl;
}
// create a copy
Tseries p(series.size());
p.copyin(series);
double bwfactor=std::sqrt(std::pow(10.,opt.decades));
// cout << bwfactor << endl;
for (int k=p.f(); k<=p.l(); ++k)
{
// center frequency
double fm=k*df;
double f1=fm/bwfactor;
double f2=fm*bwfactor;
// cout << fm << " " << f1 << " " << f2 << endl;
int k1=int(f1/df);
int k2=int(f2/df)+1;
k1 = k1 < p.f() ? p.f() : k1;
k2 = k2 > p.l() ? p.l() : k2;
TFXX_assert(k2>k1, "negative bandwidth for averaging")
Tseries sub=aff::subarray(p)(k1,k2);
series(k)=aff::func::avg(sub);
}
}
}
// adjust sampling interval to frequency interval
wid2.dt=1./(wid2.dt*wid2.nsamples);
wid2.nsamples=series.size();
......@@ -325,4 +437,49 @@ int main(int iargc, char* argv[])
} // main
/*
* Further discussion of processing concepts
* =========================================
*
* Scaling factor for Hanning taper
* --------------------------------
*
* Tapering stationary noise with a Hanning taper reduces total signal energy
* an thus the signal power obtained through an FFT. From Walter's Matlab
* scripts I took:
*
* If data is the time series vector and y is the Hanning taper of appropriate
* length, Walter calculates
*
* w = sqrt(sum(y.*y)/ndata);
*
* and
*
* y=y/w;
*
* where ndata is the length of data and y.
*
* If k is the sample index with 0 <= k < N, then
*
* w = \sqrt{\sum\limits_{k=0}^{(N-1)} \sin^4(k \pi / N) / N}
*
* = \sqrt{\int\limits_0^T \sin^4(\pi t / T) d t / T}
*
* = \sqrt{\int\limits_0^\pi \sin^4(\phi) d \phi / \pi}
*
* = \sqrt{\frac{3 \pi}{16}} = 0.58904862
*
* by using (Gradshteyn and Ryzhik, eq. 2.513 7)
*
* \int \sin^4 x d x = \frac{3 x}{8} - \frac{\sin 2 x}{4} +
* \frac{\sin 4 x}{32}
*
* and thus
*
* \int\limits_0^\pi \sin^4(\phi) d \phi = 3 \pi^2 / 16
*
* T = N \Delta t, where \Delta t is the sampling interval.
*
*/
/* ----- END OF foutra.cc ----- */
Supports Markdown
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