<jupyter_start><jupyter_text>OMP-based CS vs. PWD of Early Room Reflections in 3D  & Time Domain<NAME>, <NAME> (2019): "On the detection quality of early room reflection directions using compressive sensing on rigid spherical microphone array data", In: Proc. of 23rd Intl Congress on Acoustics, Aachen, Germany, 2676-2683May/June/September 2019OMP...orthogonal matching pursuit, cf. [<NAME>, Sparse and Redundant Representations, Springer, 2010](https://link.springer.com/book/10.1007/978-1-4419-7011-4)CS...Compressive/Compressed Sensing, cf. [<NAME>, Sparse and Redundant Representations, Springer, 2010](https://link.springer.com/book/10.1007/978-1-4419-7011-4)PWD...Plane Wave Decomposition, cf. [<NAME>, Fundamentals of Spherical Array Processing, Springer, 2015](https://link.springer.com/book/10.1007/978-3-662-45664-4) System of Equations* $\mathbf{y}$...mic array temporal responses of 3D spherical or random array* $\mathbf{A}$...sensing matrix / dictionary for 3D plane wave directions, temporal responses* $\mathbf{x}$...coefficients for sparse directions of incoming plane wavesto solve$\mathbf{y} = \mathbf{A} \cdot \mathbf{x}$ for $\mathbf{x}$ with e.g. OMPfor a certain time window Coordinate System Convention* azimuth is denoted with phi in deg, $0\leq\phi < 2 \pi$* colatitude/**polar angle** is denoted with theta in deg, $0 \leq \theta \leq \pi$* radius is denoted with r in m, $0 < r < \infty$for spherical coordinates to cartesian $x = r \sin(\theta) \cos(\phi)$,$y = r \sin(\theta) \sin(\phi)$,$z = r \cos(\theta)$ Varirables and Dimensionswe utilize the following quantities* phimic, thetamic, Rmic...azimuth, polar angle, radius of microphone array* Nmic...number of microphones in array* phipw, thetapw...azimuth, polar angle of plane wave incidence* Npw...number of plane wave directions in PWD or dictionary/sensing matrix* Nt...number of samples for acoustic impulse response, time axis $0 \leq nt < Nt$ * Ntwin...number of samples for actual window / time frame on which $\mathbf{y} = \mathbf{A} \cdot \mathbf{x}$is solved for $\mathbf{x}$ by OMP, thus we get specific $\mathbf{x}$ for time instances at $0\leq nt \leq Nt-Ntwin$* temporal sampling frequency $f_s$ in Hz* Nf...number of FFT bins for frequencies $0 \leq f \leq \frac{f_s}{2}$ for even FFT size* c...speed of sound in m/s, we assume constant 343 m/ssee a test case at the very end of this notebook, this validates OMP-CS matrix algebra with the introduced variables Critical Notes / Remarks / Observations* evaluation of 3D data (phi, theta) over time cannot be visualized in a surface plot with at the same time easy interpretation. This is why we introduced the scatter plots phi/time over level and theta/time over level* sensing matrix for our acoustic problem has very high coherence, random mic array vs. random plane wave directions only slightly improves the situation. We consider random free field mic arrays as not practical, since the setup/validation time vs. the guarantee that this setup outperforms a PWD in the scenario under discussion make the effort ratio rather poor. * early room reflections (at least in the ideal ISM) exhibit reflections from different angles but (almost) exactly at same times (you can arrange such scenarios by source and receiver position tweaking). The OMP-CS sliding window might not be the best approach to identify these temporal coincident events, unless one searches for more than one x-coefficient. But then, for how many we should look, this is a highly arbitrary choice.  Import Modules / General Variables<jupyter_code>import sys
print('needs "compress" environment (see info below): ', sys.executable)
#/Users/xxxxx/anaconda3/envs/compress/bin/python

from functools import partial
import iear
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 
import micarray  # from sfa-numpy toolbox https://github.com/spatialaudio/sfa-numpy
import numpy as np
from numpy.linalg import norm
import scipy.signal as signal
from scipy.io import loadmat
import sfs  # https://github.com/sfstoolbox/sfs-python
from sklearn.linear_model import OrthogonalMatchingPursuit
from timeit import default_timer as timer

plt.rcParams.update({'font.size': 11})

plotPWD = True
plotOMP = True<jupyter_output><empty_output><jupyter_text>Plot Spherical Mic Array Positions<jupyter_code>def plot_positions_on_sphere(phi, theta, r):
    """Plot sphere mic positions."""
    x = np.array((r * np.sin(theta) * np.cos(phi),\
           r * np.sin(theta) * np.sin(phi),\
           r * np.cos(theta)))
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(x[0,:], x[1,:], x[2,:], c='C0', marker='o')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.axis('Square')
    plt.title(x.shape)

if False:  # plot test case
    mic_array_pos = (1,2,-3)
    phi, theta, weights = micarray.modal.angular.grid_lebedev(10)
    x, y, z = sfs.util.sph2cart(phi, theta, phi*0+1)
    phi, theta, r = sfs.util.cart2sph(x + mic_array_pos[0],
                                      y + mic_array_pos[1],
                                      z + mic_array_pos[2])
    plot_positions_on_sphere(phi, theta, r)<jupyter_output><empty_output><jupyter_text>Acoustic Impulse Response of Plane Wave onto Freefield Array<jupyter_code>def array_response_3d_pw(Nt, xmic, xpw, fs=48000, c=343, osf=4):
    """Acoustic impulse response of plane wave to freefield mic array."""
    # input:
    # Nt...samples for acoustic impuse response
    # xmic...mic positions -> 
    # np.array([[phi1,theta1,r1],[phi2,theta2,r2],...,[phiNmic,thetaNmic,rNmic]]).T
    # xmic.shape = (3,Nmic)
    # xpw...one incident plane wave angle (phiPW, thetaPW, rPW), rPW is ignored
    # all angles in rad, all radii in m
    # fs...sampling frequency in Hz
    # c...speed of sound in m/s
    # osf...oversampling factor, integer>=1
    # return:
    # ir as (Nmic,Nt)
    
    Nmic = xmic.shape[1]
    ir = np.zeros((Nmic, int(osf * Nt)))

    # specific plane wave angle, radius should be unity, but must be 0<=r<oo
    phi, theta, r = xpw
    xpw = (r * np.cos(phi) * np.sin(theta),\
           r * np.sin(phi) * np.sin(theta),\
           r * np.cos(theta))
    # r of xpw is lost in subsequent code, we don't need it
    # since we only evaluate the angle in the dot product
    
    for n in range(Nmic):  # do for all requested mic positions
        phi, theta, r = xmic[:, n]
        xmic_tmp = (r * np.cos(phi) * np.sin(theta),\
                    r * np.sin(phi) * np.sin(theta),\
                    r * np.cos(theta))
        # cos of angle between xmic and xpw
        cosphi = np.dot(xmic_tmp, xpw) / (norm(xmic_tmp)*norm(xpw))
        # cos angle * norm(xPW) projects to a distance along xMic
        dr = r - cosphi * r  # note: r is allowed to vary over mic positions  
        dt = int(np.round(dr / c * fs * osf))  # corresponding time
        ir[n, dt] += osf  # insert a Dirac with unit amplitude after resampling 
    
    return signal.resample(ir, int(ir.shape[1]/osf), axis=1)

if False:  # plot test case
    fs = 48000
    r = 343 * 1/1000 # corresponds to 1 ms sound propagation
    Nt = int(4/1000*fs)
    t = np.arange(0,Nt)/fs*1000 
    
    xmic = np.array(  # mic array
        [
            [0, np.pi/2, r], # phi, theta, r
            [np.pi/2, np.pi/2, r],
            [np.pi, np.pi/2, r],
            [np.pi/4,np.pi/4, r]
        ]
    ).T
    xpw = np.array([0, np.pi/2, 1])  # plane wave incidence
    
    ir = array_response_3d_pw(Nt, xmic, xpw, fs=48000, c=343, osf=4)  

    plt.plot(t, ir[0,:],label='Mic 1')
    plt.plot(t, ir[1,:],label='Mic 2')
    plt.plot(t, ir[2,:],label='Mic 3')
    plt.plot(t, ir[3,:],label='Mic 4')
    plt.legend()
    plt.xlim(0,3)
    plt.xlabel('t / ms')
    plt.ylabel('Mic Array Response due to Plane Wave')
    print('Nt:', Nt,', xmic:', xmic.shape, ', xpw:', xpw.shape, ', ir:', ir.shape)
    # Nt: 192 , xmic: (3, 4) , xpw: (3,) , ir: (4, 192)
def plot_reflection_as_scatter(xdb, phipw, thetapw, toffs):  # check proper argument handling!!!
    #print(xdb.shape, 'max:', np.amax(xdb))
    dBclip = -10
    xnorm = xdb
    print(xnorm.shape, 'max:', np.max(xnorm))
    print('for OMP max should be 0, for PWD max might be >0', 
          'when a later reflection sums up to more energy than direct sound')
    ts = []  # sparse time sparse all entries
    phis = []  # phi sparse all entries
    thetas = []  # theta sparse all entries
    xs = []  # x in dB sparse all entries for coloring
    for n in range(xdb.shape[1]):  # get all entries that fall into the dB range:
        idx = xnorm[:,n] > dBclip
        xs_tmp = xnorm[idx, n]
        phis_tmp = phipw[idx]
        thetas_tmp = thetapw[idx]
        ts_tmp = xs_tmp*0 + n
        xs = np.insert(xs, len(xs), xs_tmp, axis=0)
        phis = np.insert(phis, len(phis), phis_tmp, axis=0)
        thetas = np.insert(thetas, len(thetas), thetas_tmp, axis=0)
        ts = np.insert(ts, len(ts), ts_tmp, axis=0)
    phis = np.mod(phis, 2*np.pi)  # in range 0...360 deg!
    print(xs.shape, ts.shape, phis.shape, thetas.shape)

    # define reference circles for min/max visualization in scatter plot
    # in top/left corner:
    xs = np.insert(xs, len(xs), (-1, dBclip), axis=0)
    phis = np.insert(phis, len(phis), (352*np.pi/180, 357*np.pi/180), axis=0)
    thetas = np.insert(thetas, len(thetas), (352/360*np.pi, 357/360*np.pi), axis=0)
    ts = np.insert(ts, len(ts), (16.5/1000*fs, 15.5/1000*fs), axis=0)

    # define radius function, size and min clip:
    scatter_radius = (xs-dBclip)*15
    scatter_radius[scatter_radius<3] = 3

    cmap = plt.get_cmap(cmap_str,-dBclip+1)  # start from 1 dB since the brightest color
    # in cmap might be to bright for paper print
    fig, (axp, axt, cax) = plt.subplots(ncols=3,figsize=(30/2.54,(30/3)/2.54), 
                      gridspec_kw={'width_ratios':[1,1, 0.05]})
    fig.subplots_adjust(wspace=0.25)
    imp = axp.scatter(ts/fs*1000, phis*180/np.pi, scatter_radius, xs, cmap=cmap, vmin=dBclip, vmax=1)
    imt = axt.scatter(ts/fs*1000, thetas*180/np.pi, scatter_radius, xs, cmap=cmap, vmin=dBclip, vmax=1)
    axp.set_yticks(np.arange(0,360+45,45))
    axt.set_yticks(np.arange(0,180+22.5,22.5))
    axt.set_yticklabels(['0','','45','','90','','135','','180'])
    axp.set_xticks(np.arange(15,60,5))
    axt.set_xticks(np.arange(15,60,5))
    axp.grid(True, which='both', color=grid_col, lw=0.25)
    axt.grid(True, which='both', color=grid_col, lw=0.25)
    axp.set_xlim(15, 40)
    axt.set_xlim(15, 40)
    axp.set_ylim(0, 360)
    axt.set_ylim(0, 180)
    axp.set_xlabel('t / ms')
    axt.set_xlabel('t / ms')
    axp.set_ylabel(r'$\phi_\mathrm{PW}$ / deg')
    axt.set_ylabel(r'$\theta_\mathrm{PW}$ / deg')
    axt.set_facecolor('w')
    axp.set_facecolor('w')
    fig.colorbar(imp, cax=cax)
    cax.set_xlabel('dB')
    
    # plot known/reference image sources:
    for n in range(len(ref_azi_is)):  
        axp.plot(ref_r_is[n]/c*1000-toffs*Rmic/c*1000,
                 np.mod(ref_azi_is[n]*180/np.pi,360),
                 '.', color=mcol2, markersize=6, fillstyle='full', mew=1)
        axp.text(ref_r_is[n]/c*1000-toffs*Rmic/c*1000,
                 np.mod(ref_azi_is[n]*180/np.pi,360), n+1,
                 fontweight='bold', fontfamily='monospace', color=mcol1)

    for n in range(len(ref_elev_is)):
        axt.plot(ref_r_is[n]/c*1000-toffs*Rmic/c*1000,
                 np.mod(ref_elev_is[n]*180/np.pi,360),
                 '.', color=mcol2, markersize=6, fillstyle='full', mew=1)
        axt.text(ref_r_is[n]/c*1000-toffs*Rmic/c*1000,
                 np.mod(ref_elev_is[n]*180/np.pi,360), n+1,
                 fontweight='bold', fontfamily='monospace', color=mcol1)<jupyter_output><empty_output><jupyter_text>General Parameter<jupyter_code># plasma, magma, viridis, cividis, inferno
# are most suitable for inspecting the sky map:
#cmap = plt.get_cmap('cividis', 12)
#cmap_str = 'cividis'
#grid_col = 'darkred'
#mcol1 = 'darkred'  # marker and grid colors should not interfere with colormap
#mcol2 = 'aliceblue'

cmap = plt.get_cmap('inferno', 12)
cmap_str = 'inferno'
grid_col = 'darkgreen'
mcol1 = 'darkgreen'
mcol2 = 'aliceblue'

interpolation_method = 'None'  # only for imshow
#shading='gouraud'  # only for pcolormesh
shading='flat'  # only for pcolormesh

Flag_Mic2D = False  # True: horizontal circle, False: Lebedev Grid
Flag_PWD2D = False  # True: 2D PW dictionary, False: 3D PW dictionary 
# either:
Flag_Testy2D = False  # 2D test data, 45deg, 128 sample offset scenario
# or:
Flag_Testy3D = False  # 3D test data for y, i.e. the reference image sources
#or both False, then we take the ISM data!<jupyter_output><empty_output><jupyter_text>Image Source Model's Sound Field onto Rigid Sphere Arraycreate measurement vector $\mathbf{y}$ Setup<jupyter_code>fs = 16000  # sampling frequency in Hz
Nt = 2400  # number of time-domain samples
c = 343  # speed of sound in m/s
f = np.fft.rfftfreq(Nt, 1/fs)
f[0] = f[1] # DC treatment
k = 2*np.pi*f/c

# room characteristics
O_sf = 15  # modal order of incident sound field, ensure resolution up to fs/2
L = [10, 11, 3]  # rectangular room, e.g. length x, width y, height z 

#DAGA 2018 data:
#coeffs = [.5, .5, .95, .95, 0.0, 0.0]  # reflection coefficient per wall
#from DAGA
#coeffs = [[0.58, 0.79, 0.90, 0.92, 0.94, 0.94],  # wood paneling 1/4"
#          [0.58, 0.79, 0.90, 0.92, 0.94, 0.94],
#          [.95, .9, 1, 1, 1, 1],  # direct sound
#          [0.58, 0.79, 0.90, 0.92, 0.94, 0.94],
#          [0.58, 0.79, 0.90, 0.92, 0.94, 0.94]
#          ]

if False:
    print('no floor / no ceiling!')
    roomtype = '4walls'
    coeffs = [.85, .85, .85, .85, 0.0, 0.0]  # reflection coefficient per wall
    # floor and ceiling totally absorbing to avoid reflections outside horizontal
    # plane, cf. this situation to a very small yard
    # with snow-covered meadow
    # with open air ceiling
    # and enclosed by four small houses' walls of plain brick with very small
    # windows, thus highly frequency independent absorption of approx. 0.1-0.15
else:
    print('all walls are brick!')
    roomtype = '6walls'
    coeffs = [.85, .85, .85, .85, 0.85, 0.85]  # reflection coefficient per wall
    
### ISM ###
max_order = 4 # maximum order of image sources, make sure that all reflections
# are included that should be investigated

# point source characteristics
source_pos = [7, 8, 1.8]

# mic array characteristics
#Rmic = 0.0875  # human head
Rmic = 0.1  # radius of array
O_mic = 10  # order of microphone array
#O_mic=2,3,4,10  for Lebedev leads to 14,26,38,170 microphones
setup = 'rigid'  # configuration of array
mic_array_pos = [2.1, 4.5, 1.8]  # mic array origin position
print('distance source to receiver:',
      norm(np.array(mic_array_pos)-np.array(source_pos)), 'm')
SNR = None  # signal peak-to-noise level ratio of microphone signals

# reference positions of image sources
ref_x_is, ref_order = sfs.util.image_sources_for_box(source_pos, L, max_order)
ref_level_is = np.prod(coeffs**ref_order, axis=1)
xd = ref_x_is - mic_array_pos
ref_azi_is, ref_elev_is, ref_r_is = sfs.util.cart2sph(xd[:,0], xd[:,1], xd[:,2])
print('max r_is:', np.max(ref_r_is), 'm, ', np.max(ref_r_is)/c*1000, 'ms, ',
     np.max(ref_r_is)/c*fs+100, 'Samples required for ISM')

tmix = iear.mixing_time(L)[0]  # 50% or 95% mixing time Lindau/Kosanke
tcut = 40
print('tmix:', tmix, ', samples at least required:', tmix/1000*fs+100)
if True:  # sort image sources, delete image sources above mixing time
    idx = sorted(range(len(ref_r_is)), key=lambda k: ref_r_is[k])
    ref_r_is = ref_r_is[idx]
    ref_azi_is = ref_azi_is[idx]
    ref_elev_is = ref_elev_is[idx]

    print('### reference image sources: ###')
    print('phi_is in deg, theta_is in deg, r_is in m, t_is in ms:')
    cutidx = (ref_r_is < tcut/1000*c) 
    ref_azi_is = ref_azi_is[cutidx]   
    ref_elev_is = ref_elev_is[cutidx]
    ref_r_is = ref_r_is[cutidx]
    
    print(len(ref_azi_is), ' image sources fall into tcut=',tcut,'ms with ISM of order', max_order)
    for n in range(len(ref_azi_is)):
        print(n+1, ',',
              np.round(ref_azi_is[n]*180/np.pi), ',', 
              np.round(ref_elev_is[n]*180/np.pi), ',',
              np.round(ref_r_is[n],3), '',
              np.round(ref_r_is[n]/c*1000,3))
if roomtype=='4walls':  # consider only theta=90deg image sources, 
    # !!! see important remarks above !!!
    # use this if last two entries in coeff are zero, since this data is used 
    # to plot the reference positions of image sources
    idx90 = np.array([q for q,v in enumerate(ref_elev_is*180/np.pi==90.) if v == True])
    ref_azi_is = ref_azi_is[idx90]
    ref_elev_is = ref_elev_is[idx90]
    ref_r_is = ref_r_is[idx90]
    print(len(ref_azi_is), ' HORIZONTAL image sources fall into tmix with ISM of order', max_order)
    for n in range(len(ref_azi_is)):
        print(n+1, ',',
              np.round(ref_azi_is[n]*180/np.pi), ',', 
              np.round(ref_elev_is[n]*180/np.pi), ',',
              np.round(ref_r_is[n],3), '',
              np.round(ref_r_is[n]/c*1000,3))

Ntwin = int(np.ceil((2*Rmic+1/2*Rmic)/c*fs))  # length of templates (& sliding window)
print('Ntwin:', Ntwin, 'samples')

t = np.arange(0,Nt)/fs
tpwd = np.arange(0,Nt)/fs
tomp = np.arange(0,Nt-Ntwin+1)/fs<jupyter_output><empty_output><jupyter_text>2D or 3D Mic Array**Important**: arrays are **centered around (0,0,0), not centered around mic_array_pos**, so pleaseensure centering around mic_array_pos if required, iear.spherical_image_sources does this job inherently<jupyter_code>if Flag_Mic2D:  # 2D horizontal equi-angular circle
    Nmic = int(360/36)
    phimic = np.arange(0, 2*np.pi, 2*np.pi/Nmic) + 2*np.pi/16
    print(phimic*180/np.pi)
    thetamic = phimic*0 + np.pi/2
    rmic = phimic*0 + Rmic
    weightmic = phimic*0 + 2*np.pi/Nmic
    xmic = np.array([phimic, thetamic, rmic])
    print('2D Mic Array Horizontal Equi-Angular Circle, dim xmic: (3, Nmic)=', xmic.shape, ', Nmic:', Nmic)
else:  # Tetrahedron or 3D Lebedev sphere
    if O_mic == 0:  # tetrahedron
        print('!!! check rotation of tetra !!!')
        phi1, theta1, r = sfs.util.cart2sph(+1,+1,+1)
        phi2, theta2, r = sfs.util.cart2sph(+1,-1,-1)
        phi3, theta3, r = sfs.util.cart2sph(-1,+1,-1)
        phi4, theta4, r = sfs.util.cart2sph(-1,-1,+1)
        phimic = np.array([phi1, phi2, phi3, phi4]).T
        thetamic = np.array([theta1, theta2, theta3, theta4]).T
        rmic = np.array([Rmic, Rmic, Rmic, Rmic]).T
        weightmic = rmic
        xmic = np.array([phimic, thetamic, rmic])
        Nmic = phimic.shape[0]
        print('3D Mic Array Tetrahedron, dim xmic: (3, Nmic)=', xmic.shape, ', Nmic:', Nmic)
        #plot_positions_on_sphere(phimic, thetamic, [np.sqrt(3), np.sqrt(3), np.sqrt(3), np.sqrt(3)])
    else:
        phimic, thetamic, weightmic = micarray.modal.angular.grid_lebedev(O_mic)
        rmic = phimic*0 + Rmic
        Nmic = phimic.shape[0]
        xmic = np.array([phimic, thetamic, rmic])
        print('3D Mic Array Lebedev Sphere, dim xmic: (3, Nmic)=', xmic.shape, ', Nmic:', Nmic)<jupyter_output><empty_output><jupyter_text>Image Source Model: Point Source in a Rect Room to a Rigid Sphere Array<jupyter_code>Y = iear.spherical_image_sources(O_sf, k, phimic, thetamic, Rmic,
                                 source_pos, mic_array_pos, L,
                                 max_order, coeffs, setup)
if SNR is not None:  # add uncorrelated noise
    Y, _ = iear.add_wgn(Y, SNR)
y = np.transpose(np.fft.irfft(Y, axis=0))
print('dim Y: (frequencies, mics):', Y.shape)
print('dim y: (mics, samples)=', y.shape)
print('frequency limit of microphone array: {:5.1f} Hz'.format(iear.upper_frequency_limit(O_mic, Rmic)))
print('frequency limit of modal representation: {:5.1f} Hz'.format(iear.upper_frequency_limit(O_sf, Rmic)))
print('Schroeder frequency of room: {:5.1f} Hz'.format(iear.schroeder_frequency(L, coeffs)))
print('mixing time of room: {0:5.1f} ms (50%) {1:5.1f} ms (95%)'.format(iear.mixing_time(L)[0],iear.mixing_time(L)[1]))
# dim Y: (frequencies, mics): (513, 170)
# dim y: (mics, samples)= (170, 1024)
# frequency limit of microphone array: 5459.0 Hz
# frequency limit of modal representation: 8188.5 Hz
# Schroeder frequency of room:  43.1 Hz
# mixing time of room:  31.2 ms (50%)  54.0 ms (95%)<jupyter_output><empty_output><jupyter_text>or 3D Test Data<jupyter_code>if Flag_Testy3D: # test with given/known y/Y under ! freefield ! conditions:
    
    # shift mic array to mic_array_pos
    x_tmp, y_tmp, z_tmp = sfs.util.sph2cart(phimic, thetamic, rmic)
    phi_tmp, theta_tmp, r_tmp = sfs.util.cart2sph(x_tmp + mic_array_pos[0],
                                                  y_tmp + mic_array_pos[1],
                                                  z_tmp + mic_array_pos[2])
    xmic_tmp = np.array((phi_tmp, theta_tmp, r_tmp))
    # which we however do not deploy in this test case:
    
    y = 0 * array_response_3d_pw(Nt, xmic, (0, 0, 1), fs)  # alloc RAM
    for n in range(len(ref_azi_is)):  # images sources
        y_tmp = array_response_3d_pw(Nt, xmic,
                                     (ref_azi_is[n], ref_elev_is[n], 1), fs)
        
        # amplitude decay point source, direct sound with unity amplitude
        # to simulate ISM this is only meaningful for IS order 1 or highly
        # reflective boundaries, where multiple reflections are only attenutated
        # by 1/r law:
        y_tmp = ref_r_is[0] * 1/(ref_r_is[n]**1) *\
        np.roll(y_tmp, int(np.round((ref_r_is[n]-Rmic)/c*fs)), 1) # **1.68 simulates the specific absorption!!!
        
        # no amplitude decay
        #y_tmp = np.roll(y_tmp, int(np.round((ref_r_is[n]-Rmic)/c*fs)), 1)
        
        y += y_tmp
        
    Y = np.transpose(np.fft.rfft(y, axis=1))
    if SNR is not None:  # add uncorrelated noise
        Y, SNR_tmp = iear.add_wgn(Y, SNR)
        y = np.transpose(np.fft.irfft(Y, axis=0))
        print('SNR RMS=', SNR_tmp, 'dB')    
    print('### 3D test data in y, Y ###')

print('dim y:', y.shape)
print('dim Y:', Y.shape)
# dim y: (72, 1024)
# dim Y: (513, 72)<jupyter_output><empty_output><jupyter_text>or 2D Test Data<jupyter_code>if Flag_Testy2D: # test with given/known y/Y under ! freefield ! conditions:
    phi_tmp = np.arange(0, 360, 45)*np.pi/180
    print(np.round(phi_tmp*180/np.pi))
    # this could have been implemented more elegant with a for loop
    # however it's convenient for easy access to individual reflections:
    y1 = array_response_3d_pw(Nt, xmic, (phi_tmp[1], np.pi/2, 1), fs)
    y1 = 1/1*np.roll(y1, 128*1, 1)

    y2 = array_response_3d_pw(Nt, xmic, (phi_tmp[2], np.pi/2, 1), fs)
    y2 = 1/1*np.roll(y2, 128*2, 1)

    y3 = array_response_3d_pw(Nt, xmic, (phi_tmp[3], np.pi/2, 1), fs)
    y3 = 1/1*np.roll(y3, 128*3, 1)

    y4 = array_response_3d_pw(Nt, xmic, (phi_tmp[4], np.pi/2, 1), fs)
    y4 = 1/1*np.roll(y4, 128*4, 1)

    y5 = array_response_3d_pw(Nt, xmic, (phi_tmp[5], np.pi/2, 1), fs)
    y5 = 1/1*np.roll(y5, 128*5, 1)
    
    y6 = array_response_3d_pw(Nt, xmic, (phi_tmp[6], np.pi/2, 1), fs)
    y6 = 1/1*np.roll(y6, 128*6, 1)
    
    y7 = array_response_3d_pw(Nt, xmic, (phi_tmp[7], np.pi/2, 1), fs)
    y7 = 1/1*np.roll(y7, 128*7, 1)    

    y = y1 + y2 + y3 + y4 + y5 + y6 + y7
    Y = np.transpose(np.fft.rfft(y, axis=1))

    if SNR is not None:  # add uncorrelated noise
        Y, SNR_tmp = iear.add_wgn(Y, SNR)
        y = np.transpose(np.fft.irfft(Y, axis=0))
        print('SNR RMS=', SNR_tmp, 'dB')
    print('### 2D Test data in y, Y ###')
    
print('dim y:', y.shape)
print('dim Y:', Y.shape)
# dim y: (72, 1024)
# dim Y: (513, 72)
fig = plt.figure(figsize=(8, 5))
plt.plot(t*1000,y.T);
plt.xlim((15,55))
plt.grid(True)
plt.xlabel('t / ms')
plt.ylabel('Mic IRs h[k]');
plt.title('AIRs of Images Sources to Mic Array')<jupyter_output><empty_output><jupyter_text>PW Directions 2D or 3D Plane Wave Dictionary<jupyter_code>if Flag_PWD2D:
    Npw = 360
    print('TBD uniform vs. random!!!')
    phipw = np.arange(0, 2*np.pi, 2*np.pi/Npw)
    phipw = np.sort(np.random.uniform(low=0, high=2*np.pi, size=(Npw,)))
    plt.plot(phipw)
    plt.xlabel('pw number')
    plt.ylabel('phipw / rad')
    thetapw = phipw*0 + np.pi/2
    rpw = phipw*0 + 1
    xpw = np.array((phipw, thetapw, rpw))
    print('### 2D PW Dict ### dim xpw: (3, Npw)=', xpw.shape, 'Npw:', Npw)
    O_pwd = int(<PASSWORD>)
else:
    O_dict = 27 # needs an even number if Gauss grid is used for horizontal circle !!!
    
    # either:
    print('Lebedev PW Dict:') 
    phipw, thetapw, _ = micarray.modal.angular.grid_lebedev(O_dict)

    # or:
    #print('Gauss Grid')
    #phipw, thetapw, _ = micarray.modal.angular.grid_gauss(O_dict)
    
    xpw = np.array((phipw, thetapw, phipw*0+1))
    Npw = xpw.shape[1]
    print('### 3D PW Dict ### dim xpw: (3, Npw)=', xpw.shape, 'Npw:', Npw)
    if O_mic==0:  # we then use a tetrahedron, which has Ambisonics order 1
        O_pwd = 1
    else:
        O_pwd = <PASSWORD>  # (O+1)^2 = Nmics criterion, here conservatively for Lebedev
    
    print('frequency limit of modal representation: {:5.1f} Hz'.format(iear.upper_frequency_limit(O_pwd, Rmic))) 

# used if Gauss grid is applied for PWD surface plot of horizontal plane
idx90 = np.array([q for q,v in enumerate(thetapw*180/np.pi==90.) if v == True])
phipw90_deg = phipw[idx90]*180/np.pi
print('Npw:', len(phipw), ', NpwHor (i.e. theta=90deg):', idx90.shape[0],
      'phi res:', 360/idx90.shape[0], 'deg')
file_id = 'Room_'+str(roomtype)+'_ISOrder_'+str(max_order)+'_Nmic_'+str(Nmic)+'_Npw_'+str(Npw)+'_SNR_'+str(SNR)
print(file_id)<jupyter_output><empty_output><jupyter_text>Plane Wave Decomposition<jupyter_code>xpwd = np.zeros((Npw,Nt))
#phimic, thetamic, _ = xmic
start = timer()
if True:
    print('### matched filter beamformer ###')
    for n in range(Npw):
        if np.mod(n, int(Npw/10))==0:
               print(n, '/', Npw, 'plane waves')
        phipw_tmp, thetapw_tmp, _ = xpw[:,n]
        XPWD = iear.pwd_matched_filter(O_pwd, k,
                                       phimic, thetamic, Rmic, weightmic,
                                       phipw_tmp, thetapw_tmp, Y, setup)
        xpwd[n,:] = np.fft.irfft(XPWD)  # write temporal response to matrix
else:
    print('### modal filter beamformer ###')
    regularizer = partial(micarray.modal.radial.regularize, a0=10, method='Tikh')
    beta = 0  #0...rect, 8.6...approx blackman 
    window = np.kaiser(2*O_pwd+1, beta)[O_pwd:]
    for n in range(Npw):
        if np.mod(n, int(Npw/10))==0:
               print(n, '/', Npw, 'plane waves')        
        phipw_tmp, thetapw_tmp, _ = xpw[:,n]
        XPWD = iear.pwd_modal(O_pwd, k,
                              phimic, thetamic, Rmic, weightmic,
                              phipw_tmp, thetapw_tmp, Y, regularizer, setup, window)
        xpwd[n,:] = np.fft.irfft(XPWD)  # write temporal response to matrix    

end = timer()
print('PWD processing time:', end - start, 'seconds')
print('dim XPWD: (Nf, plane waves)=', XPWD.shape)
print('dim xpwd: (plane waves, Nt)=', xpwd.shape)<jupyter_output><empty_output><jupyter_text>Amplitude Decay Compensation<jupyter_code>Comp1r = True

gComp1r = (tpwd*c) / ref_r_is[0]
gComp1r[tpwd<ref_r_is[0]/c] = 1
plt.plot(tpwd*1000, 20*np.log10(gComp1r))
plt.xlabel('t / ms')
plt.ylabel('A / dB')
plt.title('Compensation for 1/r')  # works only for single reflections or very low absorption coefficients
print(gComp1r.shape)
if Comp1r:
    xpwd_tmp = xpwd*gComp1r
else:
    xpwd_tmp = xpwd<jupyter_output><empty_output><jupyter_text>Scatter Plot<jupyter_code>xdB = 20*np.log10(np.abs(xpwd_tmp)+10**(-200/20))
first_refl_lvl = np.amax(xdB[:,tpwd<(ref_r_is[0]+Rmic/2)/c])
# Important note: first_refl_lvl handling works only if the direct sound is
# really the loudest entry, which might fail in highly aliased data
xdB = xdB - first_refl_lvl
print(np.amax(xdB))<jupyter_output><empty_output><jupyter_text>Normalize to Direct Sound<jupyter_code>if plotPWD:
    plot_reflection_as_scatter(xdB, phipw, thetapw, toffs = False)
    plt.savefig('PWD_Scatter_'+file_id+'_CompPntSrcDecay_'+str(Comp1r), dpi=300)<jupyter_output><empty_output><jupyter_text>Dictionary / Sensing Matrixcreate matrix $\mathbf{A}$TBD: freefield vs. rigid, regularization of in A?!?!<jupyter_code>print('TBD: A freefield vs. rigid?!?!')
A = np.zeros((Nmic, Ntwin, Npw))
print('A is Nmic x Ntwin x Npw:', A.shape, 'as init for loop')

for n in range(Npw):
    A[:,:,n] = array_response_3d_pw(Ntwin, xmic, xpw[:,n], fs)

mutual_coherence = 0
for n in range(Ntwin):
    tmp = micarray.util.coherence_of_columns(A[:,n,:])  # tmp = Nmic x Npw for each sample
    if tmp > mutual_coherence:  # get max
        mutual_coherence = tmp    
print('max mutual coherence of A:', mutual_coherence)

A = np.reshape(A,(Nmic*Ntwin, Npw))
print('A is (Nmic_x_Ntwin) x Npw:', A.shape, 'for OMP solver')<jupyter_output><empty_output><jupyter_text>OMP-Solverfind $\mathbf{x}$ of $\mathbf{y} = \mathbf{A} \mathbf{x}$<jupyter_code>omp = OrthogonalMatchingPursuit(n_nonzero_coefs=3, tol=None)
x = np.zeros((Npw, Nt-Ntwin+1))
start = timer()
for kt in range(Nt-Ntwin+1):  # do for all samples
    if np.mod(kt, 100)==0:
        print(kt, '/', Nt-Ntwin, 'time windows')
    y_win = np.hstack(y[:, kt:kt+Ntwin])  # get Mic Array phi/theta x time window in stacked form
    omp.fit(A, y_win)  # get x from OMP for that time window (y = A x)
    x[:, kt] = omp.coef_
end = timer()
print('OMP solver processing time:', end - start, 'seconds')
print('y_win:', y_win.shape,'A:', A.shape, 'x:', x.shape[0])  # check dimension consistency
print('dim x:', x.shape)<jupyter_output><empty_output><jupyter_text>Amplitude Compensation<jupyter_code>Comp1r = True

#t = np.arange(0,Nt-Ntwin+1)/fs
gComp1r = (tomp*c) / ref_r_is[0]
gComp1r[tomp<ref_r_is[0]/c] = 1
plt.plot(tomp*1000, 20*np.log10(gComp1r))
plt.xlabel('t / ms')
plt.ylabel('A / dB')
plt.title('Compensation for 1/r') # works only for single reflections or very low absorption coefficients
print(gComp1r.shape)
if Comp1r:
    x_tmp = x*gComp1r
else:
    x_tmp = x<jupyter_output><empty_output><jupyter_text>Normalize to Direct Sound<jupyter_code>xdB = 20*np.log10(np.abs(x_tmp)+10**(-200/20))
first_refl_lvl = np.amax(xdB[:,tomp<(ref_r_is[0]+Rmic/2)/c])
# Important note: first_refl_lvl handling works only if the direct sound is
# really the loudest entry, which might fail in highly aliased data
xdB = xdB - first_refl_lvl
print(np.amax(xdB))<jupyter_output><empty_output><jupyter_text>Scatter Plot<jupyter_code>if plotOMP:
    plot_reflection_as_scatter(xdB, phipw, thetapw, toffs = True)
    plt.savefig('OMP_Scatter_'+file_id+'_CompPntSrcDecay_'+str(Comp1r), dpi=300)<jupyter_output><empty_output><jupyter_text>Appendix A: Additional PlotsNote the different surface plot methods `imshow` vs. `pcolormesh` and theirinterpolation vs. shading capabilities.**2D plots needs a Gauss Grid to properly cut out the horizontal plane of the Gauss grid**We have not maintained these plots for a long time, since the above used scatter plots do their job pretty well. So, there's no guarantee that they run properly.However, we left the other plots below for documentation and inspiration of the progress. OMP-CS Sky Map as pcolormesh, 2Dentries of the Hor 2D Plane only, this is convenient for horizontal image sources only, note that single pixels might be missed due to surface plot and if a coeff x is not located in the hor planemight lead to **misconception**, thus be careful<jupyter_code>if False:
    t = np.arange(0, x.shape[1])/fs*1000
    xdb = 20*np.log10(np.abs(x[idx90])+10**(-200/20))
    vmax = np.amax(xdb)

    fig = plt.figure(figsize=(12,6))
    plt.pcolormesh(t, phipw90_deg, xdb, cmap=cmap,
                   vmin=vmax-12, vmax=vmax, rasterized=True,
                  shading=shading)
    for n in range(len(ref_azi_is)):  # plot ref image src that exist in hor plane
        if np.abs(ref_elev_is[n] - np.pi/2) < 1e-3:
            plt.plot(ref_r_is[n]/c*1000-Rmic/c*1000,
                     np.mod(ref_azi_is[n]*180/np.pi,360), 'o', color=mcol1, markersize=10, fillstyle='none')
            plt.plot(ref_r_is[n]/c*1000-Rmic/c*1000,
                     np.mod(ref_azi_is[n]*180/np.pi,360), 'o', color=mcol2, markersize=8, fillstyle='none')
    plt.yticks(np.arange(0,360,45))
    plt.xlim((15, 55))
    plt.colorbar(label='dB (1dB/col, 12dB res)')
    plt.xlabel('t / ms')
    plt.ylabel('$\phi_\mathrm{PW}$ / deg')
    plt.title('OMP-CS')
    plt.grid(True, which='both', color=grid_col, lw=0.25)

    fig = plt.figure(figsize=(12,6))
    plt.pcolormesh(t, phipw90_deg, xdb, cmap=cmap,
                   vmin=vmax-36, vmax=vmax, rasterized=True,
                  shading=shading)
    for n in range(len(ref_azi_is)):  # plot ref image src that exist in hor plane
        if np.abs(ref_elev_is[n] - np.pi/2) < 1e-3:
            plt.plot(ref_r_is[n]/c*1000-Rmic/c*1000,
                     np.mod(ref_azi_is[n]*180/np.pi,360), 'o', color=mcol1, markersize=10, fillstyle='none')
            plt.plot(ref_r_is[n]/c*1000-Rmic/c*1000,
                     np.mod(ref_azi_is[n]*180/np.pi,360), 'o', color=mcol2, markersize=8, fillstyle='none')
    plt.yticks(np.arange(0,360,45))
    plt.xlim((15, 55))
    plt.colorbar(label='dB (3dB/col, 36dB res)')
    plt.xlabel('t / ms')
    plt.ylabel('$\phi_\mathrm{PW}$ / deg')
    plt.title('OMP-CS')
    plt.grid(True, which='both', color=grid_col, lw=0.25)<jupyter_output><empty_output><jupyter_text>OMP-CS Sky Map as scatter, 2Dentries of the Hor 2D Plane only, this is convenient for horizontal image sources onlysingle entries in the hor plane are not missed due to scatter plot, however entries from OMP that are not in hor plane are missedmight lead to **misconception**, be careful<jupyter_code>if False:
    dBclip = -10
    xnorm = xdb - np.max(xdb)
    print(xnorm.shape)

    ts = []  # sparse time all entries
    ps = []  # sparse phi all entries
    xs = []  # sparse x in dB all entries for coloring
    for n in range(len(phipw90_deg)):
        idx = xnorm[n,:]>dBclip
        xs_tmp = xnorm[n, idx]
        ts_tmp = t[idx]  # time sparse slice
        ps_tmp = ts_tmp*0 + phipw90_deg[n]  # phi sparse slice
        ts = np.insert(ts, len(ts), ts_tmp, axis=0)
        ps = np.insert(ps, len(ps), ps_tmp, axis=0)
        xs = np.insert(xs, len(xs), xs_tmp, axis=0)
    print(len(ts), len(ps))    

    cmap = plt.get_cmap('inferno', -dBclip)
    fig = plt.figure(figsize=(12,6))
    plt.scatter(ts, ps, (xs-dBclip)*30, xs, cmap=cmap, vmin=dBclip, vmax=0)
    for n in range(len(ref_azi_is)):  # plot ref image src that exist in hor plane
        if np.abs(ref_elev_is[n] - np.pi/2) < 1e-3:
            plt.plot(ref_r_is[n]/c*1000-Rmic/c*1000,
                     np.mod(ref_azi_is[n]*180/np.pi,360), 'o', color=mcol1, markersize=10, fillstyle='none')
            plt.plot(ref_r_is[n]/c*1000-Rmic/c*1000,
                     np.mod(ref_azi_is[n]*180/np.pi,360), 'o', color=mcol2, markersize=8, fillstyle='none')
    plt.yticks(np.arange(0,360,45))
    plt.xlim((15, 55))
    plt.colorbar(label='dB')
    plt.xlabel('t / ms')
    plt.ylabel('$\phi_\mathrm{PW}$ / deg')
    plt.title('OMP-CS entries in hor 2D plane')
    plt.grid(True, which='both', color=grid_col, lw=0.25)
    ax = plt.gca()
    ax.set_facecolor('w')
    print(np.max(xs-dBclip), np.min(xs-dBclip))<jupyter_output><empty_output><jupyter_text>PWD Sky Map as imshow, 3D this visualization is not convenient in terms of plane wave index (encoding phi / theta), but it potentially shows all 3D entries in xhowever, due to imshow single pixels might be missed in visualizationmight lead to **misconception**, be careful<jupyter_code>if False:
    xdb = 20*np.log10(np.abs(xpwd)+10**(-200/20))
    vmax = np.amax(xdb)
    fig = plt.figure(figsize=(12, 6))
    plt.imshow(xdb, cmap, 
               vmin=vmax-12, vmax=vmax, rasterized=True, aspect='auto',
              interpolation=interpolation_method);
    plt.xlim((-10+128*0, -10+128*3))
    plt.colorbar(label='dB')
    plt.xlabel('time index nt')
    plt.ylabel('plane wave index')
    plt.title('PWD')
    plt.grid(True, which='both', color=grid_col, lw=0.25)

    fig = plt.figure(figsize=(12, 6))
    plt.imshow(xdb, cmap, 
               vmin=vmax-12, vmax=vmax, rasterized=True, aspect='auto',
              interpolation=interpolation_method);
    plt.xlim((-10+128*3, -10+128*6))
    plt.colorbar(label='dB')
    plt.xlabel('time index nt')
    plt.ylabel('plane wave index')
    plt.title('PWD')
    plt.grid(True, which='both', color=grid_col, lw=0.25)

    fig = plt.figure(figsize=(12, 6))
    plt.imshow(xdb, cmap, 
               vmin=vmax-12, vmax=vmax, rasterized=True, aspect='auto',
              interpolation=interpolation_method);
    plt.xlim((-10+128*6, -10+128*9))
    plt.colorbar(label='dB')
    plt.xlabel('time index nt')
    plt.ylabel('plane wave index')
    plt.title('PWD')
    plt.grid(True, which='both', color=grid_col, lw=0.25)<jupyter_output><empty_output><jupyter_text>CS-OMP Sky Map as imshow, 3D this visualization is not convenient in terms of plane wave index (encoding phi / theta), but it potentially shows all 3D entries in xhowever, due to imshow single pixels might be missed in visualization, especially for CS-OMP very criticalmight lead to **misconception**, be careful<jupyter_code>if False:
    xdb = 20*np.log10(np.abs(x)+10**(-200/20))
    vmax = np.amax(xdb)
    fig = plt.figure(figsize=(12, 6))
    plt.imshow(xdb, cmap, 
               vmin=vmax-12, vmax=vmax, rasterized=True, aspect='auto',
              interpolation=interpolation_method);
    plt.xlim((-10+128*0, -10+128*3))
    plt.colorbar(label='dB')
    plt.xlabel('time index nt')
    plt.ylabel('pw index')
    plt.title('OMP-CS')
    plt.grid(True, which='both', color=grid_col, lw=0.25)

    fig = plt.figure(figsize=(12, 6))
    plt.imshow(xdb, cmap, 
               vmin=vmax-12, vmax=vmax, rasterized=True, aspect='auto',
              interpolation=interpolation_method);
    plt.xlim((-10+128*3, -10+128*6))
    plt.colorbar(label='dB')
    plt.xlabel('time index nt')
    plt.ylabel('pw index')
    plt.title('OMP-CS')
    plt.grid(True, which='both', color=grid_col, lw=0.25)

    fig = plt.figure(figsize=(12, 6))
    plt.imshow(xdb, cmap, 
               vmin=vmax-12, vmax=vmax, rasterized=True, aspect='auto',
              interpolation=interpolation_method);
    plt.xlim((-10+6*128, -10+9*128))
    plt.colorbar(label='dB')
    plt.xlabel('time index nt')
    plt.ylabel('pw index')
    plt.title('OMP-CS')
    plt.grid(True, which='both', color=grid_col, lw=0.25)<jupyter_output><empty_output><jupyter_text>CS-OMP Sky Map as scatter, 3Dthis is convenient in terms of showing the sparsity of the data, but interpretationin terms of phi/theta is very hard, thus two scatter plots for individual phi and thetais finally proposed in `plot_reflection_as_scatter(xdb, phipw, thetapw)`<jupyter_code>if False:
    xdb = 20*np.log10(np.abs(x)+10**(-200/20))
    dBclip = -10
    xnorm = xdb - np.max(xdb)
    print(xdb.shape, xnorm.shape)

    ts = []  # sparse time all entries
    ps = []  # sparse phi/theta all entries
    xs = []  # sparse x in dB all entries for coloring
    for n in range(xdb.shape[0]):
        idx = xnorm[n,:]>dBclip
        xs_tmp = xnorm[n, idx]
        ts_tmp = t[idx]  # sparse time slice
        ps_tmp = ts_tmp*0 + n  # sparse phi slice
        ts = np.insert(ts, len(ts), ts_tmp, axis=0)
        ps = np.insert(ps, len(ps), ps_tmp, axis=0)
        xs = np.insert(xs, len(xs), xs_tmp, axis=0)
    print(len(ts), len(ps))    

    cmap = plt.get_cmap('inferno', -dBclip)
    fig = plt.figure(figsize=(12,6))
    plt.scatter(ts, ps, (xs-dBclip)*30, xs, cmap=cmap, vmin=dBclip, vmax=0)
    plt.ylim((0,xdb.shape[0]))
    plt.xlim((15, 55))
    plt.colorbar(label='dB')
    plt.xlabel('t / ms')
    plt.ylabel('pw index')
    plt.title('OMP-CS')
    plt.grid(True, which='both', color=grid_col, lw=0.25)
    ax = plt.gca()
    ax.set_facecolor('w')
    print(np.max(xs-dBclip), np.min(xs-dBclip))<jupyter_output><empty_output><jupyter_text>Check Dictionary Entries w.r.t. Time Alignment and Level in Detailangular characteristics cannot be estimated in these plots, for that use the above surface/scatter plotswith the scatter plot routines, these plots are somewhat obsolete<jupyter_code>if False:  # PWD
    fig = plt.figure(figsize=(8, 5))
    plt.plot(np.abs(xpwd.T));
    plt.grid(True)
    plt.xlabel('time index nt')
    plt.ylabel('pwd |x| (linear)');
    plt.title('PWD');
if False:  # PWD
    xdb = 20*np.log10(np.abs(xpwd.T)+10**(-200/20))
    vmax = np.amax(xdb)
    fig = plt.figure(figsize=(8, 5))
    plt.plot(xdb);
    plt.ylim((vmax-dBres, vmax+1))
    plt.grid(True)
    plt.xlabel('time index nt')
    plt.ylabel('pwd x in dB');
    plt.title('PWD');    
if False:  # OMP
    fig = plt.figure(figsize=(8, 5))
    plt.plot(np.abs(x.T));
    plt.grid(True)
    plt.xlabel('time index nt')
    plt.ylabel('pw dict |x| (linear)')
    plt.title('OMP-CS');
if False:  # OMP
    xdb = 20*np.log10(np.abs(x.T)+10**(-200/20))
    vmax = np.amax(xdb)
    fig = plt.figure(figsize=(8, 5))
    plt.plot(xdb);
    plt.ylim((vmax-dBres, vmax+1))
    plt.grid(True)
    plt.xlabel('time index nt')
    plt.ylabel('pw dict |x| in dB')
    plt.title('OMP-CS');<jupyter_output><empty_output><jupyter_text>Appendix B: Example to Check Handling and Matrix Dimensions for CSchose dimensions are not related to the simulations given in the paper<jupyter_code>if False:
    c = 343  # speed of sound in m/s
    fs = 10000  # sampling frequency in Hz
    Rmic = c*1/1000  # radius mic array in m
    Nt = 1000  # samples
    Ntwin = int(np.ceil(2*Rmic/c*fs)+80)  # length of templates (& sliding window)
    print('fs:', fs,'Hz, Rmic:', Rmic, 'm, Nt:', Nt, 'samples, Ntwin:', Ntwin, 'samples')
    # fs: 10000 Hz, Rmic: 0.343 m, Nt: 1000 samples, Ntwin: 100 samples
if False:
    phipw, thetapw, _ = micarray.modal.angular.grid_lebedev(7)
    if False:  # test case
        phipw = np.arange(0,360,45)*np.pi/180
        thetapw = phipw*0 + np.pi/2
    rpw = phipw * 0 + 1
    xpw = np.array((phipw, thetapw, rpw))
    Npw = xpw.shape[1]
    print('dim xpw: (3, Npw)=', xpw.shape, ', ', Npw, 'plane wave directions in dictionary/sensing matrix')
    # dim xpw: (3, Npw)= (3, 86) ,  86 plane wave directions in dictionary/sensing matrix
if False:
    phimic, thetamic, _ = micarray.modal.angular.grid_lebedev(2)
    if False:  # test case
        phimic = np.arange(0,360,90)*np.pi/180
        thetamic = phimic*0 + np.pi/2
    rmic = phimic * 0 + 1
    xmic = np.array((phimic, thetamic, rmic))
    Nmic = xmic.shape[1]
    print('dim xmic: (3, Nmic)=', xmic.shape, ', ', Nmic, 'microphones')
    # dim xmic: (3, Nmic)= (3, 14) ,  14 microphones
if False:
    A = np.zeros((Nmic, Ntwin, Npw))
    print('dim A: (Nmic , Ntwin, Npw)=', A.shape, 'as init to get all dict PWs')
    for n in range(Npw):
        A[:,:,n] = array_response_3d_pw(Ntwin, xmic, xpw[:,n], fs)
    A = np.reshape(A, (Nmic*Ntwin, Npw))
    print('dim A: (Nmic * Ntwin, Npw)=', A.shape, 'reshaped for OMP solver')
    # dim A: (Nmic , Ntwin, Npw)= (14, 100, 86) as init to get all dict PWs
    # dim A: (Nmic * Ntwin, Npw)= (1400, 86) reshaped for OMP solver
if False:
    phi = 0
    theta = np.pi/2
    r = 1
    xtst = (phi, theta, r)
    y = array_response_3d_pw(Nt, xmic, xtst, fs)
    print('dim y: (Nmic, Nt)=', y.shape)
    # dim y: (Nmic, Nt)= (14, 1000)
if False:
    # do OMP for the first time frame
    k = 0
    omp = OrthogonalMatchingPursuit(n_nonzero_coefs=1, tol=None)
    x = np.zeros((Npw, Nt-Ntwin+1))
    y_win = np.hstack(y[:, k:k+Ntwin])
    omp.fit(A, y_win)
    x[:, k] = omp.coef_
    print('Nmic * Ntwin =', Ntwin*Nmic)
    print('dim y_win: (Nmic * Ntwin, 1)=', y_win.shape,
          '\ndim dim A: (Nmic * Ntwin, Npw)=', A.shape,
          '\ndim x slice : (Npw, 1)=(', x.shape[0], ', 1)')
    print('dim x: (Npw x Time Frames)=', x.shape)
    # Nmic * Ntwin = 1400
    # dim y_win: (Nmic * Ntwin, 1)= (1400,) 
    # dim dim A: (Nmic * Ntwin, Npw)= (1400, 86) 
    # dim x slice : (Npw, 1)=( 86 , 1)
    # dim x: (Npw x Time Frames)= (86, 901)<jupyter_output><empty_output><jupyter_text>References<jupyter_code>\begin{thebibliography}{99}

\bibitem{Jarrett2017}
<NAME>.; <NAME>.; <NAME>. Theory and Applications of Spherical Microphone Array Processing. Springer, (CH), 2017.

\bibitem{Rafaely2015}
<NAME>. Fundamentals of Spherical Array Processing. Springer, Berlin (GER), 2015.

\bibitem{Elad2010}
<NAME>. Sparse and Redundant Representations. Springer, NY (USA), 2010.

\bibitem{Gerstoft2018}
<NAME>.; <NAME>.; <NAME>.; <NAME>. Introduction to compressive sensing in acoustics. J. Acoust. Soc. Am., Vol 143 (6), 2018, pp 3731-3736.

\bibitem{Xenaki2014}
<NAME>.; <NAME>.; <NAME>. Compressive beamforming. J. Acoust. Soc. Am., Vol 136 (1), 2014, pp 260-271.

\bibitem{Spors2018}
<NAME>.; <NAME>. On the Estimation of Acoustic Reflection Coefficients from In-Situ Measurements using a Spherical Microphone Array. Proc. of the 44th Deutsche Jahrestagung für Akustik (DAGA), Munich (GER), March 2018, pp 1326-1329,

\bibitem{Mignot2013}
<NAME>.; <NAME>.; <NAME>. Room Reverberation Reconstruction: Interpolation of the Early Part Using Compressed Sensing. {IEEE} Audio, Speech, Language Process., Vol 21 (11), 2013, pp 2301-2312.

\bibitem{Xenaki2018}
<NAME>.; <NAME>.; <NAME>. Sound source localization and speech enhancement with sparse Bayesian learning beamforming.
J. Acoust. Soc. Am., Vol 143 (6), 2018, pp 3912-3921.

\bibitem{Verburg2018}
<NAME>.; <NAME>. Reconstruction of the sound field in a room using compressive sensing. J. Acoust. Soc. Am., Vol 143 (6), 2018, pp 3770-3779.

\bibitem{Lindau2012}
<NAME>.; <NAME>.; <NAME>. Perceptual Evaluation of Model- and Signal-Based Predictors of the Mixing Time in Binaural Room Impulse Responses.
J. Audio Eng. Soc., Vol 60 (11), 2012, pp 889-898.

\end{thebibliography}<jupyter_output><empty_output><jupyter_text>Python Environmentneeds `iear.py`sfstoolbox commit: ddbe35d3f6128e860a1fc8b3b2a7b0e857015064sfa-numpy commit: bff5737ef429f31228d20a9e1d0ce7d46d3080d3<jupyter_code>Package           Version  Location                           
----------------- -------- -----------------------------------
appdirs           1.4.3    
appnope           0.1.0    
attrs             19.1.0   
backcall          0.1.0    
bleach            3.1.0    
certifi           2019.3.9 
chardet           3.0.4    
cycler            0.10.0   
dask              1.2.2    
decorator         4.4.0    
defusedxml        0.5.0    
entrypoints       0.3      
idna              2.8      
imageio           2.5.0    
ipykernel         5.1.0    
ipython           7.5.0    
ipython-genutils  0.2.0    
jedi              0.13.3   
Jinja2            2.10.1   
jsonschema        3.0.1    
jupyter-client    5.2.4    
jupyter-core      4.4.0    
jupyterlab        0.35.6   
jupyterlab-server 0.2.0    
kiwisolver        1.1.0    
MarkupSafe        1.1.1    
matplotlib        3.0.3    
mistune           0.8.4    
mpmath            1.1.0    
nbconvert         5.5.0    
nbformat          4.4.0    
networkx          2.3      
notebook          5.7.8    
numpy             1.15.1   
olefile           0.46     
orthopy           0.5.3    
pandocfilters     1.4.2    
parso             0.4.0    
pexpect           4.7.0    
pickleshare       0.7.5    
Pillow            6.0.0    
pip               19.1     
pipdate           0.3.2    
prometheus-client 0.6.0    
prompt-toolkit    2.0.9    
ptyprocess        0.6.0    
pycparser         2.19     
Pygments          2.4.0    
pyparsing         2.4.0    
pyrsistent        0.15.2   
python-dateutil   2.8.0    
PyWavelets        1.0.3    
pyzmq             18.0.1   
quadpy            0.11.3   
requests          2.21.0   
scikit-image      0.13.1   
scikit-learn      0.19.1   
scipy             1.2.1    
Send2Trash        1.5.0    
setuptools        41.0.1   
sfa-numpy         0.0.1    /Users/xxxxx/git/sfa-numpy         
sfs               0.4.0    /Users/xxxxx/git/sfs/sfs-python    
six               1.12.0   
sounddevice       0.3.12   /Users/xxxxx/git/python-sounddevice
sympy             1.4      
terminado         0.8.2    
testpath          0.4.2    
toolz             0.9.0    
tornado           6.0.2    
traitlets         4.3.2    
urllib3           1.24.3   
wcwidth           0.1.7    
webencodings      0.5.1    
wheel             0.33.4   
appdirs==1.4.3
appnope==0.1.0
attrs==19.1.0
backcall==0.1.0
bleach==3.1.0
certifi==2019.3.9
chardet==3.0.4
cycler==0.10.0
dask==1.2.2
decorator==4.4.0
defusedxml==0.5.0
entrypoints==0.3
idna==2.8
imageio==2.5.0
ipykernel==5.1.0
ipython==7.5.0
ipython-genutils==0.2.0
jedi==0.13.3
Jinja2==2.10.1
jsonschema==3.0.1
jupyter-client==5.2.4
jupyter-core==4.4.0
jupyterlab==0.35.6
jupyterlab-server==0.2.0
kiwisolver==1.1.0
MarkupSafe==1.1.1
matplotlib==3.0.3
mistune==0.8.4
mpmath==1.1.0
nbconvert==5.5.0
nbformat==4.4.0
networkx==2.3
notebook==5.7.8
numpy==1.15.1
olefile==0.46
orthopy==0.5.3
pandocfilters==1.4.2
parso==0.4.0
pexpect==4.7.0
pickleshare==0.7.5
Pillow==6.0.0
pipdate==0.3.2
prometheus-client==0.6.0
prompt-toolkit==2.0.9
ptyprocess==0.6.0
pycparser==2.19
Pygments==2.4.0
pyparsing==2.4.0
pyrsistent==0.15.2
python-dateutil==2.8.0
PyWavelets==1.0.3
pyzmq==18.0.1
quadpy==0.11.3
requests==2.21.0
scikit-image==0.13.1
scikit-learn==0.19.1
scipy==1.2.1
Send2Trash==1.5.0
-e git+https://github.com/spatialaudio/sfa-numpy.git@<PASSWORD>#egg=sfa_numpy
-e git+https://github.com/sfstoolbox/sfs-python.git@c67db62da81d0e63515bbe86bdda548684851edd#egg=sfs
six==1.12.0
-e git+https://github.com/spatialaudio/python-sounddevice.git@cc5c8e39d40c38e69da727e14a33d6be48ff4132#egg=sounddevice
sympy==1.4
terminado==0.8.2
testpath==0.4.2
toolz==0.9.0
tornado==6.0.2
traitlets==4.3.2
urllib3==1.24.3
wcwidth==0.1.7
webencodings==0.5.1
name: compress
channels:
  - conda-forge
  - defaults
dependencies:
  - appnope=0.1.0=py36_1000
  - attrs=19.1.0=py_0
  - backcall=0.1.0=py_0
  - blas=1.1=openblas
  - bleach=3.1.0=py_0
  - ca-certificates=2019.3.9=hecc5488_0
  - certifi=2019.3.9=py36_0
  - cycler=0.10.0=py_1
  - dask-core=1.2.2=py_0
  - decorator=4.4.0=py_0
  - defusedxml=0.5.0=py_1
  - entrypoints=0.3=py36_1000
  - freetype=2.10.0=h24853df_0
  - imageio=2.5.0=py36_0
  - ipykernel=5.1.0=py36h24bf2e0_1002
  - ipython=7.5.0=py36h24bf2e0_0
  - ipython_genutils=0.2.0=py_1
  - jedi=0.13.3=py36_0
  - jinja2=2.10.1=py_0
  - jpeg=9c=h1de35cc_1001
  - jsonschema=3.0.1=py36_0
  - jupyter_client=5.2.4=py_3
  - jupyter_core=4.4.0=py_0
  - jupyterlab=0.35.6=py36_0
  - jupyterlab_server=0.2.0=py_0
  - kiwisolver=1.1.0=py36h770b8ee_0
  - libblas=3.8.0=5_hd44dcd8_netlib
  - libcblas=3.8.0=5_hd44dcd8_netlib
  - libcxx=8.0.0=2
  - libcxxabi=8.0.0=2
  - libffi=3.2.1=h6de7cb9_1006
  - libgfortran=3.0.1=0
  - liblapack=3.8.0=5_hd44dcd8_netlib
  - libpng=1.6.37=h2573ce8_0
  - libsodium=1.0.16=h1de35cc_1001
  - libtiff=4.0.10=h79f4b77_1001
  - markupsafe=1.1.1=py36h1de35cc_0
  - matplotlib=3.0.3=py36_1
  - matplotlib-base=3.0.3=py36habd5b94_1
  - mistune=0.8.4=py36h1de35cc_1000
  - nbconvert=5.5.0=py_0
  - nbformat=4.4.0=py_1
  - ncurses=6.1=h0a44026_1002
  - networkx=2.3=py_0
  - notebook=5.7.8=py36_0
  - numpy=1.15.1=py36_blas_openblashd3ea46f_1
  - olefile=0.46=py_0
  - openblas=0.2.20=8
  - openssl=1.1.1b=h01d97ff_2
  - pandoc=2.7.2=0
  - pandocfilters=1.4.2=py_1
  - parso=0.4.0=py_0
  - pexpect=4.7.0=py36_0
  - pickleshare=0.7.5=py36_1000
  - pillow=6.0.0=py36h7095ceb_0
  - pip=19.1=py36_0
  - prometheus_client=0.6.0=py_0
  - prompt_toolkit=2.0.9=py_0
  - ptyprocess=0.6.0=py_1001
  - pygments=2.4.0=py_0
  - pyparsing=2.4.0=py_0
  - pyrsistent=0.15.2=py36h01d97ff_0
  - python=3.6.7=h8dc6b48_1004
  - python-dateutil=2.8.0=py_0
  - pywavelets=1.0.3=py36h1e5eb4f_1
  - pyzmq=18.0.1=py36h2d07e9b_1
  - readline=7.0=hcfe32e1_1001
  - scikit-image=0.13.1=py36_0
  - scikit-learn=0.19.1=py36_blas_openblas_201
  - scipy=1.2.1=py36hbd7caa9_1
  - send2trash=1.5.0=py_0
  - setuptools=41.0.1=py36_0
  - six=1.12.0=py36_1000
  - sqlite=3.26.0=h1765d9f_1001
  - terminado=0.8.2=py36_0
  - testpath=0.4.2=py_1001
  - tk=8.6.9=ha441bb4_1001
  - toolz=0.9.0=py_1
  - tornado=6.0.2=py36h01d97ff_0
  - traitlets=4.3.2=py36_1000
  - wcwidth=0.1.7=py_1
  - webencodings=0.5.1=py_1
  - wheel=0.33.4=py36_0
  - xz=5.2.4=h1de35cc_1001
  - zeromq=4.3.1=h0a44026_1000
  - zlib=1.2.11=h1de35cc_1004
  - pip:
    - appdirs==1.4.3
    - chardet==3.0.4
    - idna==2.8
    - mpmath==1.1.0
    - orthopy==0.5.3
    - pipdate==0.3.2
    - quadpy==0.11.3
    - requests==2.21.0
    - sympy==1.4
    - urllib3==1.24.3<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1489160097028502, 'fraction_numerical': 0.05889933292904791, 'mean_word_length': 3.0786056577523575, 'pattern_counts': {'":': 0, '<': 161, '<?xml version=': 0, '>': 153, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 3}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '171140', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 22840, 'n_tokens_neox': 20974, 'n_words': 4793, 'path': 'estimate_early_reflections_OMP.ipynb', 'repo_name': 'spatialaudio/doa-early-room-reflections-pwd-vs-omp', 'size': '52947'}