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

works

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: 1615
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 0cb7165e
# ---------------------------------------
#
# $Id: Makefile,v 1.21 2005-01-11 15:11:57 tforb Exp $
# $Id: Makefile,v 1.22 2005-01-12 15:37:19 tforb Exp $
#
# Makefile fuer tools /src/ts/wf
#
......@@ -64,8 +64,8 @@ tidofi fredofi sigfit: %: %.o
newprog $@
phasetest%.sff: phasedsignals
phasedsignals $@ -o -p $(patsubst phasetest%.sff,%,$@) -n 16 \
-t 10. -b 20.
phasedsignals $@ -v -o -p $(patsubst phasetest%.sff,%,$@) -n 16 \
-t 10. -b 3. -s 0. -O 2 -N 6 -h 0.5
P%.sff: %.sff; evelo $< $@ -P
%.s.ps: %.sff; splot $< 0. $@/ps 1.e-2
......
c this is <phasedsignals.f>
c ----------------------------------------------------------------------------
c ($Id: phasedsignals.f,v 1.4 2005-01-11 18:36:45 tforb Exp $)
c ($Id: phasedsignals.f,v 1.5 2005-01-12 15:37:20 tforb Exp $)
c
c Copyright (c) 2004 by Thomas Forbriger (BFO Schiltach)
c
......@@ -19,25 +19,28 @@ c
& //'compute synthetic signals with given phase')
character*(*) PHASEDSIGNALS_CVS_ID
parameter(PHASEDSIGNALS_CVS_ID=
& '$Id: phasedsignals.f,v 1.4 2005-01-11 18:36:45 tforb Exp $')
& '$Id: phasedsignals.f,v 1.5 2005-01-12 15:37:20 tforb Exp $')
c
c
logical overwrite
logical overwrite, arg_usedamp
character arg_model
double precision arg_duration, arg_delay, arg_bandwidth
integer arg_npower, arg_nord
double precision arg_damp
integer arg_npower, arg_nord, arg_npo
c
integer nsamples, maxsamples, maxfreq, nfreq, nord
integer nsamples, maxsamples, maxfreq, nfreq, nord, npo
parameter(maxsamples=100000,maxfreq=((maxsamples/2)+1))
double precision dt, t0, df, pi, pi2, f ,t, bw, fac, v
double precision dt, t0, df, pi, pi2, f, bw, fac
parameter(pi=3.141592653589793d0,pi2=2.d0*pi)
double precision amp(maxfreq), phase(maxfreq)
double precision sigtospec, sigtotime
double precision sigtospec, sigtotime, calfac, damp
parameter (sigtospec=-1.d0, sigtotime=1.d0)
real fdata(maxsamples), tf_rand
integer idata(maxsamples)
equivalence (fdata,idata)
double complex spectrum(maxsamples), ime, pfac
double complex calspec(maxsamples), lp
complex tf_lpbut
parameter(ime=(0.d0,1.d0))
integer lu, i
parameter(lu=12)
......@@ -47,14 +50,15 @@ c debugging
c commandline
integer maxopt, lastarg, iargc
character*80 argument
parameter(maxopt=9)
parameter(maxopt=11)
character*2 optid(maxopt)
character*40 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c here are the keys to our commandline options
data optid/2h-d, 2h-v, 2h-o, 2h-p, 2h-t, 2h-n, 2h-s, 2h-b, 2h-O/
data opthasarg/3*.FALSE.,6*.TRUE./
data optarg/3*1h-,1h0,2h1.,2h10,2h0.,2h1.,1h2/
data optid/2h-d, 2h-v, 2h-o, 2h-p, 2h-t, 2h-n, 2h-s, 2h-b, 2h-O,
& 2h-N, 2h-h/
data opthasarg/3*.FALSE.,8*.TRUE./
data optarg/3*1h-,1h0,2h1.,2h10,2h0.,2h1.,1h4,1h1,2h.7/
c
c------------------------------------------------------------------------------
c basic information
......@@ -66,6 +70,7 @@ c
print *,version
print *,'Usage: phasedsignals outfile [-o] [-p 0|m|r] [-O o]'
print *,' [-t d] [-n p] [-s t] [-b bw]'
print *,' [-N o] [-h h]'
print *,' or: phasedsignals -help'
if (argument(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
......@@ -82,6 +87,9 @@ c
print *,'-s t shift data by t'
print *,'-b bw set bandwidth'
print *,'-O o use low pass repsonse of order ''o'' '
print *,'-N o take filter coeff. to the power of ''o'' '
print *,'-h h take second order filter with damping'
print *,' ''h'' of critical'
print *,' '
print *,PHASEDSIGNALS_CVS_ID
stop
......@@ -102,6 +110,9 @@ c
read(optarg(7), *) arg_delay
read(optarg(8), *) arg_bandwidth
read(optarg(9), *) arg_nord
read(optarg(10), *) arg_npo
read(optarg(11), *) arg_damp
arg_usedamp=optset(11)
c
c------------------------------------------------------------------------------
c go
......@@ -114,65 +125,79 @@ c go
nfreq=(nsamples/2)+1
bw=arg_bandwidth
nord=arg_nord
npo=arg_npo
c
if (verbose) then
print *,'prepare spectrum of Butterworth Lowpass'
print *,'at ',bw,' Hz with order ',nord
endif
c
do i=1,nfreq
f=(i-1)*df
amp(i)=max(1.e8*exp(-(f/bw)**2),1.)
spectrum(i)=
if (arg_usedamp) then
spectrum(i)=lp(f,bw,arg_damp)**npo
else
spectrum(i)=tf_lpbut(sngl(f),sngl(bw),nord)**npo
endif
amp(i)=abs(spectrum(i))
calspec(i)=spectrum(i)
enddo
c
if (verbose) print *,'prepare calibration factor'
call ccspectrum(calspec,nsamples)
call tf_dfork(nsamples,calspec,sigtotime)
calfac=0.d0
do i=1,nsamples
calfac=max(calfac,abs(calspec(i)))
enddo
c
if (arg_model.eq.'0') then
if (verbose) print *,'prepare zero phase'
do i=1,nfreq
phase(i)=0.d0
enddo
elseif (arg_model.eq.'r') then
if (verbose) print *,'prepare random phase'
call tf_tsrand()
do i=1,nfreq
phase(i)=pi2*tf_rand()
enddo
elseif (arg_model.eq.'m') then
do i=1,nfreq
spectrum(i)=log(abs(amp(i)))
enddo
do i=2,nfreq-1
spectrum(nsamples+2-i)=log(abs(amp(i)))
enddo
call tf_dfork(nsamples,spectrum,sigtospec)
do i=2,nfreq-1
spectrum(i)=ime*spectrum(i)
enddo
spectrum(nfreq)=(0.d0,0.d0)
do i=nfreq+1,nsamples
spectrum(i)=-ime*spectrum(i)
enddo
call tf_dfork(nsamples,spectrum,sigtotime)
endif
c
if ((arg_model.eq.'0').or.(arg_model.eq.'r')) then
if (verbose) print *,'apply new phase'
do i=1,nfreq
phase(i)=real(spectrum(i))
pfac=exp(ime*phase(i))
spectrum(i)=amp(i)*pfac
enddo
else
stop 'ERROR: illegal model selection!'
endif
c
c
if (verbose) print *,'apply delay'
do i=1,nfreq
f=(i-1)*df
phase(i)=phase(i)-f*pi2*t0
phase(i)=-f*pi2*t0
pfac=exp(ime*phase(i))
spectrum(i)=amp(i)*pfac
if (abs(abs(pfac)-1.d0).gt.1.e-8) then
print *,i,amp(i),phase(i),pfac,ime
endif
spectrum(i)=spectrum(i)*pfac
enddo
c
if (verbose) print *,'complement spectrum'
call ccspectrum(spectrum,nsamples)
if (verbose) print *,'transform spectrum to time domain'
call tf_dfork(nsamples,spectrum,sigtotime)
fac=sqrt(float(nsamples))*df*pi2
fac=sqrt(float(nsamples))*df*pi2/calfac
fac=1.d0/calfac
if (verbose) print *,'normalize signal'
do i=1,nsamples
fdata(i)=real(spectrum(i))*fac
enddo
if (overwrite) call sff_New(lu,filename,i)
if (verbose) print *,'open file ',
& filename(1:index(filename,' ')-1)
call sffu_simpleopen(lu,filename)
if (verbose) print *,'write file'
call sffu_simplewrite(lu, .true., fdata, nsamples,
& sngl(dt), sngl(t0))
if (verbose) print *,'finished'
c
stop
end
......@@ -191,5 +216,20 @@ c
enddo
return
end
c----------------------------------------------------------------------
double complex function lp(omega,omegan,h)
double precision omega,omegan,h
double complex res, ime
parameter (ime=(0.d0,1.d0))
double precision q,rt
q=omega/omegan
rt=sqrt(1-h**2)
res=-1.d0/((q-ime*h-rt)*(q-ime*h+rt))
lp=res
return
end
c
c ----- END OF phasedsignals.f -----
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