This tutorial will demonstrate: how to use the forward and inverse operations of the SCDT in the the PyTransKit package.
Parameters: Inputs: ---------------- reference : 1d array Reference signal used to calculate SCDT. x0 : 1d array Independent axis variable of reference signal. Functions:
Forward transform: Iposhat, Ineghat, Imasspos, Imassneg = stransform(I,x)
Inputs:
----------------
I : 1d array
Input signal to transform.
x : 1d array
Independent axis variable of the input signal to transform.
Outputs:
----------------
Iposhat : 1d array
CDT of positive values of input signal I.
Ineghat : 1d array
CDT of negative values of input signal I.
Imasspos : 1d array
Total mass of the positive part of the input signal I.
Imassneg : 1d array
Total mass of the negative part of the input signal I.
Inverse transform: I_recon = istransform(Iposhat, Ineghat, Imasspos, Imassneg)
Inputs:
----------------
Iposhat : 1d array
CDT of positive values of input signal I.
Ineghat : 1d array
CDT of negative values of input signal I.
Imasspos : 1d array
Total mass of the positive part of the input signal I.
Imassneg : 1d array
Total mass of the negative part of the input signal I.
Outputs:
----------------
I_recon : 1d array
Reconstructed signal.
The examples will cover the following operations:
By default, maps the SCDT values between $[0,1]$
*Import necessary libraries*
import numpy as np
from numpy import interp
import math
import matplotlib.pyplot as plt
from scipy.linalg import lstsq
from scipy import signal
import sys
sys.path.append('../')
from pytranskit.optrans.continuous import SCDT
*Reference signal $s_0$*
N = 300
## Gaussian distribution
t0 = np.linspace(0,1,N) # Domain of the reference
#s0 = np.exp(-10*(t0-0.5)**2)/np.sqrt(2*np.pi)
s0 = np.ones(N)
s0 = s0/s0.sum()
# scdt = SCDT(reference=s0,x0=t0)
plt.plot(t0,s0)
plt.show()
Target signal $s$
for example $\quad $ $s(t) = \pi^{-0.25}e^{4\pi j \omega (t-t_{\mu})}e^{-120(t-t_{\mu})^2} $
t = np.linspace(0,0.5,N) # Domain of the input signal
def signal_gabor(t):
w = 10
tm = t[-1]/2
s = np.real(math.pi**-0.25 * np.exp(4*math.pi*1j*w*(t-tm)) * np.exp(-120*((t-tm)**2)))
return s
s = signal_gabor(t)
s_CDF = np.cumsum(s)
# Plotting
plt.plot(t, s)
plt.title("signal",
fontsize=16, color='black')
plt.xlabel('t: time')
plt.ylabel('s(t)')
plt.show()
Compute SCDT transform of $s$ & the inverse transform
scdt = SCDT(reference=s0)
#Calculating SCDT of the target signal
##spos,sneg contain transform of the positve and negative parts of s
### smasspos, smassneg contain the masses of the positive and negative parts of s
spos, sneg, smasspos, smassneg = scdt.stransform(s)
fig, ax = plt.subplots(1, 4, sharex=False, sharey=False, figsize=(15,5))
ax[0].plot(t, np.cumsum(s), 'b-')
ax[0].set_xlabel('t: time')
ax[0].set_title('Original integral of $s$')
ax[0].legend(['$F_s(x)=\int_{-\infty}^x s(t)dt$'])
ax[1].plot(t0, spos, 'g-')
ax[1].set_title('CDT of positive part')
ax[1].legend(['$\hat{spos}$'])
ax[2].plot(t0, sneg, 'g-')
ax[2].set_title('CDT of negative part')
ax[2].legend(['$\hat{sneg}$'])
ax[3].plot(t, np.cumsum(s), 'b-')
ax[3].plot(t, scdt.istransform(spos, sneg, smasspos, smassneg), 'r-')
ax[3].set_title('Reconstruction')
ax[3].legend(['$F_s$','Inverse of SCDT of $s$'])
plt.show()
Maps the SCDT values to the domain $t$ of the input signal $s$
Target signal $s$
t = np.linspace(0,0.5,N) # Domain of the input signal
def signal_gabor(t):
w = 10
tm = t[-1]/2
s = np.real(math.pi**-0.25 * np.exp(4*math.pi*1j*w*(t-tm)) * np.exp(-120*((t-tm)**2)))
return s
s = signal_gabor(t)
s_CDF = np.cumsum(s)
# Plotting
plt.plot(t, s)
plt.title("signal",
fontsize=16, color='black')
plt.xlabel('t: time')
plt.ylabel('s(t)')
plt.show()
Compute SCDT transform of $s$ & the inverse transform
scdt = SCDT(reference=s0)
#Calculating SCDT of the input signal s
##spos,sneg contain transform of the positve and negative parts of s
### smasspos, smassneg contain the masses of the positive and negative parts of s
#### SCDT values will be mapped to the domain t of input signal s
spos, sneg, smasspos, smassneg = scdt.stransform(s, x=t)
fig, ax = plt.subplots(1, 4, sharex=False, sharey=False, figsize=(15,5))
ax[0].plot(t, np.cumsum(s), 'b-')
ax[0].set_xlabel('t: time')
ax[0].set_title('Original integral of $s$')
ax[0].legend(['$F_s(x)=\int_{-\infty}^x s(t)dt$'])
ax[1].plot(t0, spos, 'g-')
ax[1].set_title('CDT of positive part')
ax[1].legend(['$\hat{spos}$'])
ax[2].plot(t0, sneg, 'g-')
ax[2].set_title('CDT of negative part')
ax[2].legend(['$\hat{sneg}$'])
ax[3].plot(t, np.cumsum(s), 'b-')
ax[3].plot(t, scdt.istransform(spos, sneg, smasspos, smassneg), 'r-')
ax[3].set_title('Reconstruction')
ax[3].legend(['$F_s$','Inverse of SCDT of $s$'])
plt.show()
Test linear separability of three classes, i.e. class 𝑘∈{1,2,3} k ∈ { 1 , 2 , 3 } , Gabor, Sawtooth, Square apodized wave. Use only translation and scaling.
## Define a set of signals
t = np.linspace(-0.5,0.5,N)
def signal_gabor(t):
w = 30
s = np.real(math.pi**-0.25 * np.exp(2*math.pi*1j*w*(t)) * np.exp(-200*((t)**2)))
return s
def signal_sawtooth(t):
s2 = signal.sawtooth(2 * np.pi * 15 * t)* np.exp(-250*((t)**2))
return s2
def signal_square(t):
s3 = signal.square(2 * np.pi * 15 * t)*np.exp(-250*((t)**2))
return s3
s1 = signal_gabor(t)
s2 = signal_sawtooth(t)
s3 = signal_square(t)
## Plotting
fig, ax = plt.subplots(1, 3, sharex=False, sharey=False, figsize=(25,3))
#c = ['b*', 'ro', 'kx'][i]
ax[0].plot(t,s1,'b')
ax[0].set_title("Gabor wave")
ax[0].set_yticks([])
ax[1].plot(t,s2,'r')
ax[1].set_title("Sawtooth wave")
ax[1].set_yticks([])
ax[2].plot(t,s3,'k')
ax[2].set_title("Square wave")
ax[2].set_yticks([])
plt.show()
K = 3 # number of classes
L_slope = 20 # number of slope variations
L_shift = 30 # number of shifts
L = L_slope*L_shift
# N number of points in each signal
I = np.zeros((K,L,N))
Ihat = np.zeros((K,L,2*(N+1))) #we are going to store both positive and negative parts of the SCDT together by concatenating them
slope_array = np.linspace(0.75,2,L_slope)
shift_array = np.linspace(-0.25,0.25,L_shift)
#populate gabor class
k = 0
v = 0
for i_slope in range(L_slope):
for i_shift in range(L_shift):
b = slope_array[i_slope]
c = shift_array[i_shift]
g = b*(t) + c
g_prime = b
s = g_prime*signal_gabor(g)
s = s/1
noise = np.random.normal(0,0.02,N)
I[k,v,:] = s + noise
Ipos, Ineg, Imasspos, Imassneg = scdt.stransform(s+noise)
Ipos_mass = np.array([Imasspos])
Ineg_mass = np.array([Imassneg])
Ipos_final = np.concatenate((Ipos,Ipos_mass),axis=0)
Ineg_final = np.concatenate((Ineg,Ineg_mass),axis=0)
Ihat[k,v,:]= np.concatenate((Ipos_final,Ineg_final),axis=0)
#Ihat[k,v,:] = np.concatenate([Ipos,Ineg])
v = v+1
#populate sawthoot class
k = 1
v = 0
for i_slope in range(L_slope):
for i_shift in range(L_shift):
b = slope_array[i_slope]
c = shift_array[i_shift]
g = b*(t) + c
g_prime = b
s = g_prime*signal_sawtooth(g)
s = s/1
noise = np.random.normal(0,0.02,N)
I[k,v,:] = s + noise
Ipos, Ineg, Imasspos, Imassneg = scdt.stransform(s+noise)
Ipos_mass = np.array([Imasspos])
Ineg_mass = np.array([Imassneg])
Ipos_final = np.concatenate((Ipos,Ipos_mass),axis=0)
Ineg_final = np.concatenate((Ineg,Ineg_mass),axis=0)
Ihat[k,v,:]= np.concatenate((Ipos_final,Ineg_final),axis=0)
#Ihat[k,v,:] = np.concatenate([Ipos,Ineg])
v = v+1
#populate square wave class
k = 2
v = 0
for i_slope in range(L_slope):
for i_shift in range(L_shift):
b = slope_array[i_slope]
c = shift_array[i_shift]
g = b*(t) + c
g_prime = b
s = g_prime*signal_square(g)
s = s/1
noise = np.random.normal(0,0.02,N)
I[k,v,:] = s + noise
Ipos, Ineg, Imasspos, Imassneg = scdt.stransform(s+noise)
Ipos_mass = np.array([Imasspos])
Ineg_mass = np.array([Imassneg])
Ipos_final = np.concatenate((Ipos,Ipos_mass),axis=0)
Ineg_final = np.concatenate((Ineg,Ineg_mass),axis=0)
Ihat[k,v,:]= np.concatenate((Ipos_final,Ineg_final),axis=0)
#Ihat[k,v,:] = np.concatenate([Ipos,Ineg])
v = v+1
import random
N_plots = 7
index_array = np.zeros((3,N_plots),dtype=int)
index_array[0,:] = random.sample(range(0, L_slope*L_shift), N_plots)
index_array[1,:] = random.sample(range(0, L_slope*L_shift), N_plots)
index_array[2,:] = random.sample(range(0, L_slope*L_shift), N_plots)
## Plotting
fig, axs = plt.subplots(N_plots, 3, sharex=True, sharey=True, figsize=(25,3*N_plots))
color_array = ['b','r','k']
for k in range(3):
for p in range(N_plots):
v = index_array[k,p]
s = I[k,v,:]
axs[p,k].plot(t,s,color_array[k])
plt.show()
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.svm import LinearSVC
from sklearn.model_selection import train_test_split
import warnings
warnings.filterwarnings('ignore')
fontSize=14
X=np.reshape(I,(K*L,N)) #Combine the signals into a features vector X
Xhat=np.reshape(Ihat,(K*L,2*(N+1))) #Combine the transformed signals into a features vector Xhat
label=np.concatenate((np.zeros(L,),np.ones(L,),-1*np.ones(L,))) # Define the labels as -1,0,1 for the three classes
label = label.astype(np.int)
#split traning and testing data
data_train, data_test, label_train, label_test = train_test_split(X, label, test_size=0.5, random_state=42)
lda=LinearDiscriminantAnalysis(n_components=2)
svm=LinearSVC()
fig,ax=plt.subplots(1,2,figsize=(15,5))
dataLDA=lda.fit_transform(data_train,label_train)
dataLDA_test = lda.transform(data_test)
print("Test accuracy in signal space:",lda.score(data_test, label_test) )
for i in range(3):
l = [0, 1, -1][i]
c = ['b*', 'ro', 'kx'][i]
ax[0].plot(dataLDA_test[label_test==l,0],dataLDA_test[label_test==l,1],c)
x_min, x_max = dataLDA_test.min(axis=0)[0], dataLDA_test.max(axis=0)[0]
y_min, y_max = dataLDA_test.min(axis=0)[1], dataLDA_test.max(axis=0)[1]
nx, ny = 400, 200
xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),np.linspace(y_min, y_max, ny))
svm.fit(dataLDA,label_train)
Z = svm.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z[:].reshape(xx.shape)
ax[0].pcolormesh(xx, yy, Z,cmap='OrRd')
ax[0].contour(xx, yy, Z, linewidths=.5, colors='k')
ax[0].set_xticks([])
ax[0].set_yticks([])
ax[0].set_title("LDA projections in signal space",fontsize=fontSize)
lda=LinearDiscriminantAnalysis(n_components=2)
svm=LinearSVC()
data_train_h, data_test_h, label_train, label_test = train_test_split(Xhat, label, test_size=0.5, random_state=42)
dataLDAh=lda.fit_transform(data_train_h,label_train)
dataLDA_test_h = lda.transform(data_test_h)
print("Test accuracy in transform space:",lda.score(data_test_h, label_test) )
for i in range(3):
l = [0, 1, -1][i]
c = ['b*', 'ro', 'kx'][i]
ax[1].plot(dataLDA_test_h[label_test==l,0],dataLDA_test_h[label_test==l,1],c)
x_min, x_max = dataLDA_test_h.min(axis=0)[0], dataLDA_test_h.max(axis=0)[0]
y_min, y_max = dataLDA_test_h.min(axis=0)[1], dataLDA_test_h.max(axis=0)[1]
nx, ny = 400, 200
xx, yy = np.meshgrid(np.linspace(x_min, x_max, nx),np.linspace(y_min, y_max, ny))
svm.fit(dataLDAh,label_train)
Z = svm.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z[:].reshape(xx.shape)
ax[1].pcolormesh(xx, yy, Z,cmap='OrRd')
ax[1].contour(xx, yy, Z, linewidths=.5, colors='k')
ax[1].set_xticks([])
ax[1].set_yticks([])
ax[1].set_title("LDA projections in transform space",fontsize=fontSize)
plt.show()
Test accuracy in signal space: 0.32 Test accuracy in transform space: 0.9988888888888889