0

I am trying to rotate some seismic data held in a numpy nd-array. This array has dimensions (N-receiver, M-sources,3-source_channels, 3-receiver channels, K-time channels).

I know how to set up the rotation if applied to a single time stamp (t_i) over a single receiver & single source station. The actual Z, R, T, N, E, notation isn't important for the general problem, just know that the transformation is defined like so:

In python, for a single timestamp I'd probably code up something like:

import numpy as np
a = 50.0 # example alpha
b = 130  # example beta
a_rotation = np.asarray([[1,0,0],[0,np.cos(a),np.sin(a)],[0,-np.sin(a),np.cos(a)]])
b_rotation = np.asarray([[1,0,0],[0,-np.cos(b),-np.sin(b)],[0,np.sin(b),-np.cos(b)]])

# pretend the zn's are actual float data
single_timeslice_data = np.asarray([[zz,zn,ze],[nz,nn,ne],[ez,en,ee]])

# operation w numpy matrix rotation

rotated_channels = a_rotation @ single_timeslice_data @ b_rotation

So my question is twofold:

  1. How do I apply this matrix product with numpy over all timesteps?

ex: (3 x 3) * (3 x 3 x K) * (3 x 3)

  1. How do I do this matrix product with numpy when there may be an arbitrary number of other dimensions?

ex: (3 x 3) * (N x M x 3 x 3 x K) * (3 x 3)

KEVIN MENDOZA
  • 75
  • 7
  • The basics of `@` is it does `dot`, matrix multiplication on the last 2 dimensions, and treats all the earlier ones as 'batch' dimensions, obeying the normal numpy broadcasting rules. – hpaulj Nov 01 '19 at 19:53

1 Answers1

1

1)
(3 x 3) * (3 x 3 * K) * (3 * 3) = (3 x 3 x K) 

np.einsum('ab,bcK,cd->adK', Arr1, Arr2, Arr3)

2)
(3 x 3) * (N x M x 3 x 3 x K) * (3 x 3) = (N x M x 3 x 3 x K)

np.einsum('ab,NMbcK,cd->NMadK', Arr1, Arr2, Arr3)
Attack68
  • 4,437
  • 1
  • 20
  • 40