R wavelets

In [1]:

    

genrick <- function(duration, dt, f) {
    # RICKER WAVELET GENERATOR
    # Written by Leonard Lisapaly (leonardl@fisika.ui.ac.id)
    # Converted to R by J Lees, 2004
    #
    # INPUTS
    # freq = wavelet dominant frequency [Hz]
    # dt   = sampling interval [sec]
    # nw   = duration of wavelet [sec]
    nw = duration / dt
    a  = f * sqrt(pi) / 2
    nc = (nw + 1) / 2
    tc = (nc - 1) * dt
    t  = seq(from=0, length=nw-1) * dt
    b  = pi * f * (t - tc)
    w  = a * (1 - 2*b^2) * exp(-b^2)
    return(w)
}


    

In [2]:

    

imp <- matrix(rep(1, len=51), nrow = 51) * 2550 * 2650

# Impedance, imp: VP    RHO
imp[c(10:15),] =  2700 * 2750
imp[c(15:27),] =  2400 * 2450
imp[c(27:35),] =  2800 * 3000


    

In [3]:

    

wavelet <- genrick(0.040, 0.001, 100)


    

In [4]:

    

matplot(wavelet)


    

In [5]:

    

# Python equivalents
#          imp[1:]         imp[:-1]          imp[1:]          imp[:-1]
m <- (tail(imp, -1) - head(imp, -1)) /  (tail(imp, -1) + head(imp, -1))


    

In [6]:

    

x <- c(1, -1, 0, 0, 0, 0)
head(toeplitz(x), -1)


    

    







	
 1
-1
 0
 0
 0
 0
	
-1
 1
-1
 0
 0
 0
	
 0
-1
 1
-1
 0
 0
	
 0
 0
-1
 1
-1
 0
	
 0
 0
 0
-1
 1
-1

In [8]:

    

convmtx <- function(h, n) {
    col_1 = np.r_[h[0], np.zeros(n-1)]
    row_1 = np.r_[h, np.zeros(n-1)]
    return(toeplitz(col_1, row_1))
}


    

In [ ]:

    

convmtx()


    

In [ ]:

    




    

In [44]:

    

# Or...
solve(A, b)


    

    
Out[44]:





	
1
	
2

In [43]:

    

A %*% x - b


    

    
Out[43]:






	
0
	
0

In [ ]: