Thread: Moving average?

Dear All -

This is something that is easy using the features available in MacSAC but otherwise difficult. First the easy solution. It uses the SUM ON/OFF feature of CHANGESTACK, and the HDRNUM function that returns the trace numbers of all traces whose BAZ header field is within PM degrees of the value given by AZ:

-----
#cat dostack
* azstack az x pm y
$keys az pm
sss; changestack all sum off
do nf list (hdrnum baz az $az$ pm $pm$)
changestack $nf$ sum on
enddo
sumstack
qs
-----

This lists a macro (called dostack) that will do a stack of traces available in memory (it assumes you've used SAC's ADDSTACK in SSS to add the traces in with the proper static delays). Then (in SAC)

m dostack az 10 pm 5
m dostack az 20 pm 5
...
gets you a succession of stacks. Note that the data files are only read once, and manipulated in memory. See HELP EXPRESSIONS for a description of HDRNUM.

The more difficult way re-reads the files for every stack at a different azimuth. It should work with any version of SAC, but relies on an external program to read a SAC file header and get a value from it (the program here is called "sachdrinfo" and it is one packaged with MacSAC; there are are probably others available on the web). It assumes your files are all selected by a shell command "ls *.bhz" which you need to change to select the proper set of your own files. It uses a shell script to calculate the azimuth difference from the target (the azimuth difference is calculated by the dot product to avoid branch cuts in azimuth) and to build up a set of addstack commands for each stack. The shell script is called "dostack-other"

------
#cat dostack-other
#!/bin/sh
#The first shell script arg is the target azimuth, and second arg is the range
# around the target azimuth.
az=${1:-0} pm=${2:-10}
ls *.bhz | sachdrinfo -s -input baz |
# sachdrinfo output is file name followed by azimuth in header
awk 'func cosd(x) {return cos(x*pi/180)}
func sind(x) {return sin(x*pi/180)}
BEGIN{pi=4*atan2(1,1); az='"${az}"'; un=0+'"${pm}"'
vx=sind(az); vy=cosd(az)
print "zerostack"
}
{dx=sind($2); dy=cosd($2); dotp=vx*dx+vy+dy; ang=atan2(sqrt(1-dotp*dotp),dotp)
if (ang*180/pi <= un) print "addstack",$1,"delay 0 s"
}'
-----

Next, in SAC:

sc sh dostack-other 10 5 > /tmp/sac.in
sss; m /tmp/sac.in; sumstack; qs
sc sh dostack-other 20 5 > /tmp/sac.in
sss; m /tmp/sac.in; sumstack; qs
...

I tested neither of these -- they are only quickly written down. You will have to be careful to check that they do what you want. However, it shows you the basic ideas and the difference between the complexity of the two solutions.

On 5 Jul 2012, at 08:25, Alireza Alinaghi wrote:

Dear Colleagues:

Does anyone know how the moving average is done in SAC?

I mean how can I move a window of lets say 10 (degrees, in azimuth) in steps of 5 (degrees) across the read in files and
calculate the stacked trace for each window?

Is there any macro for doing this?

Thank you in advance for your reply.

Alireza Alinaghi


_______________________________________________
sac-help mailing list
sac-help<at>iris.washington.edu
http://www.iris.washington.edu/mailman/listinfo/sac-help

George Helffrich
george.helffrich<at>bris.ac.uk