Import necessary libraries


In [1]:
import numpy as np, cmath
import scipy.io
from matplotlib import pyplot as plt

from numpy import pi, sin, cos, exp, log, random  #import basic functions from numpy that we'll need

%matplotlib inline

Import optional library Seaborn for pretty plots.


In [2]:
import seaborn as sns
sns.set_palette('muted')
sns.set_style('darkgrid')

Figure 11.1


In [3]:
srate = 1000. #sampling rate of 1kHz
time = np.arange(-1,1 +1/srate,1/srate) 
freq = 10 #in Hz
amp = 2 #amplitude of sine wave

sine_wave = amp*sin(2*pi*freq*time)

plt.figure()
plt.plot(time,sine_wave)
plt.axis([0,1,-3,3])
_=plt.title("My first sine wave")


Figure 11.2


In [4]:
srate = 500. #sampling rate in Hz

#create arrays of frequencies, amplitudes, and phases to plot
frex = np.array([3, 10, 5 ,15, 35])
amplit = np.array([5,15,10,5,7])
phases = np.pi*np.array([1/7.,1/8.,1.,1/2.,-1/4.])

time = np.arange(-1,1 +1/srate,1/srate) 


sine_waves = np.zeros([len(frex),len(time)])
for fi in xrange(len(frex)):
    sine_waves[fi,:] = amplit[fi] * sin(2*pi*frex[fi]*time +phases[fi])

#plot each sine wave individually
plt.figure()

for fi in xrange(len(frex)):
    plt.subplot(len(frex),2,2*(fi+1)-1)
    plt.plot(sine_waves[fi,:],linewidth=2)
        
#plot the sum of all sum waves
plt.subplot(1,2,2)
plt.plot(np.sum(sine_waves,axis=0))
plt.tight_layout()    
_=plt.title("Sum of sines")


Figure 11.3


In [5]:
plt.figure()
noise = 5 * random.randn(np.shape(sine_waves)[1])
plt.plot( np.sum(sine_waves,axis=0) + noise )
_=plt.title("sum of sine waves plus white noise")


Figure 11.4


In [6]:
time = np.arange(-1,1+1/srate,1/srate)

#create three sine waves
s1 = sin(2*pi*3*time)
s2 = 0.5 * sin(2*pi*8*time)
s3 = s1 + s2

#plot the sine waves
f = plt.figure()

for ii in xrange(3):    
    plt.subplot(2,3,ii+1)
    
    #plot the sinusoids using the eval command, which evaluates the input string
    eval("plt.plot(time,s"+str(ii+1)+")")
    plt.axis([0,1,-1.6,1.6])
    plt.yticks(np.arange(-1.5,2,.5))
    
    plt.subplot(2,3,ii+1+3)
    
    f = eval("np.fft.fft(s" + str(ii+1) +")/float(len(time))")
    hz = np.linspace(0,srate/2.,np.floor(len(time)/2.)+1) # we only have resolution up to SR/2 (Nyquist theorem)
    plt.bar(hz,np.absolute(f[:len(hz)]*2))
    plt.axis([0,11,0,1.2])
    plt.xticks(np.arange(0,11))
    
plt.tight_layout()


Figure 11.5


In [8]:
from mpl_toolkits.mplot3d import Axes3D
N = 10 #length of sequence
data = random.randn(N) #create random numbers, sampled from normal distribution
srate = 100 #sampling rate in Hz
nyquist = srate/2 #Nyquist frequency -- highest frequency you can measure the data

#initialize matrix for Fourier output

frequencies = np.linspace(0,nyquist,N/2+1)
time = np.arange(N)/float(N)

#Fourier transform is dot product between sine wave and data at each frequency
fourier = np.zeros(N)*1j #create complex matrix
for fi in xrange(N):
    sine_wave = np.exp(-1j *2 *pi*fi*time)
    fourier[fi] = np.sum(sine_wave*data)
    
fourier = fourier/float(N)


fig=plt.figure(figsize=(8,8))

plt.subplot(221)
plt.plot(data,"-o")
plt.xlim([0,N+1])
plt.title("Time domain representation of data")

ax = fig.add_subplot(222, projection='3d')
ax.plot(frequencies,np.angle(fourier[:N/2+1]),zs=np.absolute(fourier[:N/2+1])**2)
ax.view_init(20, 20)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Phase')
ax.set_zlabel('Power')
plt.title("3D representation of the Fourier transform")


plt.subplot(223)
plt.bar(frequencies,np.absolute(fourier[:N/2+1])**2)
plt.xlim([-5, 105])
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power")
plt.title("Power spectrum")

plt.subplot(224)
plt.bar(frequencies,np.angle(fourier[:N/2+1]))
plt.xlim([-5, 105])
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase angle")
plt.yticks(np.arange(-pi,pi,pi/2.))

plt.title("Phase spectrum")

plt.tight_layout()


Figure 11.6


In [9]:
reconstructed_data = np.zeros(N)

for fi in xrange(N):
    #scale sine wave by fourier coefficient
    sine_wave = fourier[fi] * exp(1j *2 * pi * fi * time)
    
    
    #sum the sine waves together, and take only the real part
    reconstructed_data += np.real(sine_wave)
    
plt.figure()
plt.plot(data,'ko',linewidth=4)
plt.plot(reconstructed_data,'r-*')
_=plt.legend(["original data","inverse Fourier transform"])


Figure 11.7


In [12]:
fft_data = np.fft.fft(data)/N

plt.figure()

plt.subplot(131)
plt.plot(frequencies,np.absolute(fourier[:N/2+1])**2,'ko',markersize=8,linewidth=5)
plt.plot(frequencies,np.absolute(fft_data[:N/2+1])**2,'r*-',markersize=1)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power")
plt.title("Power spectrum")

plt.subplot(132)

#wrap the computed phase angles because +pi = -pi
phase_manual = (np.angle(fourier[:N/2+1])  + pi) % (2 * pi ) - pi
phase_fft = (np.angle(fft_data[:N/2+1])  + pi) % (2 * pi ) - pi

plt.plot(frequencies,phase_manual,'ko',markersize=8,linewidth=5)
plt.plot(frequencies,phase_fft,'r*-',markersize=1)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase Angle")
plt.title("Phase spectrum")

plt.tight_layout()

#below line is a little verbose, but it is just real(ifft(fft(data)))
ifftData = np.real(np.fft.ifft(np.fft.fft(data))) 

plt.subplot(133)
plt.plot(reconstructed_data,'ko',markersize=8,linewidth=5)
plt.plot(ifftData,'r*-',markersize=1)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power")
plt.title("Data reconstruction")

plt.legend(["Manual Fourier transform","Numpy FFT"], bbox_to_anchor = (2, 1))

plt.tight_layout()


Figure 11.9


In [13]:
#create array of frequencies, amplitudes, and phases
frex = np.array([3,10,5,7])
amplit = np.array([5,15,10,5])
phases = pi*np.array([1/7.,1/8.,1,1/2.])

#create a time series of sequenced sine waves

srate = 500.
time = np.arange(-1,1+1/srate,1/srate)
stationary = np.zeros(len(time)*len(frex))
nonstationary = np.zeros(len(time)*len(frex))

for fi in xrange(len(frex)):
    #compute sine wave
    temp_sine_wave = amplit[fi] * sin(2*pi*frex[fi]*time+phases[fi])
    
    #enter into stationary time series
    stationary = stationary + np.tile(temp_sine_wave,(1,len(frex)))
    
    #optional change of amplitude over time
    temp_sine_wave *=time+1
    
    #start and stop indices for insertion of sine wave
    start_idx = fi * len(time)
    stop_idx = fi *len(time) + len(time)
    
    #enter into non-stationary time series
    nonstationary[start_idx:stop_idx] = temp_sine_wave

    
plt.figure()
plt.subplot(221)
plt.plot(stationary[0],'r')
plt.axis([0,stationary.shape[1],-30,30])
plt.title("stationary signal")
_,xticks=plt.xticks(np.arange(0,4000,500))
plt.setp(xticks,rotation=-45)

plt.subplot(222)
plt.plot(nonstationary)
plt.xlim([1,len(nonstationary)])
plt.title("non-stationary signal")
_,xticks=plt.xticks(np.arange(0,4000,500))
plt.setp(xticks,rotation=-45)

frequencies = np.linspace(0,srate/2,len(nonstationary)/2+1)
fft_nonstationary = np.fft.fft(nonstationary)/len(nonstationary)
fft_stationary = np.fft.fft(stationary[0])/stationary.shape[1]

plt.subplot(212)
plt.plot(frequencies,np.absolute(fft_stationary[:len(frequencies)]*2),'r')
plt.plot(frequencies,np.absolute(fft_nonstationary[:len(frequencies)]*2))
plt.xlim([0,np.max(frex)*2])
plt.legend(["Power stationary","Power non-stationary"])

plt.tight_layout()


Figure 11.10


In [14]:
data=scipy.io.loadmat("sampleEEGdata.mat")


eegdat4convol=np.squeeze(data["EEG"][0,0]["data"][46,:,0])
srate = float(data["EEG"][0,0]["srate"][0,0])

#create a Gaussian
time = np.arange(-1,1+1/srate,1/srate)
s = 5/float(2*pi*30) #standard deviation
gaussian = 1/30.*exp((-time**2)/(2*s**2)) 


plt.figure()
plt.subplot(231)
plt.plot(eegdat4convol)

plt.subplot(234)
plt.plot(gaussian)

plt.subplot(232)
convolvedSig=np.convolve(eegdat4convol,gaussian,mode='same')
plt.plot(convolvedSig)

plt.subplot(235)
plt.plot(np.absolute(np.fft.fft(convolvedSig)))


plt.subplot(233)
plt.plot(np.absolute(np.fft.fft(eegdat4convol)))


plt.subplot(236)
plt.plot(np.absolute(np.fft.fft(gaussian)))


plt.tight_layout()


Figure 11.11


In [15]:
srate = 1000.
time = np.arange(-0.5,0.5,1/srate)
f = 20
fg = np.array([15, 5])
s = sin(2*pi*f*time)

def plotConvolutionTheorem(ii=0):
    assert(ii<3 and ii>=0)
    
    g = exp((-time**2)/(2*(4/((2*pi*fg[ii])**2))))/fg[ii]
    plt.figure()

    plt.subplot(411)
    plt.plot(time,s)
    plt.title("Sine wave (signal)")

    plt.subplot(412)
    plt.plot(time,g)
    plt.title("Gaussian (kernel)")

    plt.subplot(413)
    plt.plot(time,np.convolve(s,g,mode="same"))
    plt.title("convolved signal")

    plt.subplot(427)
    fft_s = np.absolute(np.fft.fft(s))
    fft_s = fft_s[:np.floor(len(fft_s)/2)+1]/np.max(fft_s[:np.floor(len(fft_s)/2)+1])
    plt.plot(np.arange(0,501),fft_s,'r')

    fft_g = np.absolute(np.fft.fft(g))
    fft_g = fft_g[:np.floor(len(fft_g)/2)+1] / np.max(fft_g[:np.floor(len(fft_g)/2)+1])
    plt.plot(np.arange(0,501),fft_g)
    plt.axis([0,40,0,1.05])
    plt.title("individual power spectra")


    plt.subplot(428)
    plt.bar(np.arange(0,501),fft_g*fft_s)
    plt.axis([0,40,0,0.035])
    plt.title("multiplied power spectra")

    plt.tight_layout()

plotConvolutionTheorem(ii=0)



In [16]:
plotConvolutionTheorem(ii=1)


Figure 11.12


In [43]:
EEGtimes = data["EEG"][0,0]["times"][0]
srate = float(data["EEG"][0,0]["srate"][0,0])

time = np.arange(-1,1+1/srate,1/srate)
s = 5/(2*pi*30)
gaussian = exp((-time**2)/(2*s**2))/30.


plt.figure()

#plot EEG data
plt.subplot(411)
plt.plot(EEGtimes,eegdat4convol,'r')

#plot Gaussian
plt.subplot(412)
plt.plot(time,gaussian)

#plot convolved
plt.subplot(413)
plt.plot(EEGtimes,np.convolve(eegdat4convol,gaussian,mode="same"))

# Plot power spectrum of signal
plt.subplot(427)
nfft = len(eegdat4convol)
fft_s = np.absolute(np.fft.fft(eegdat4convol,nfft))
fft_s = fft_s[:np.floor(nfft/2)+1]
f = np.linspace(0,srate/2,np.floor(nfft/2)+1)
plt.plot(f,fft_s/np.max(fft_s),'r')

#plot power spectrum of Gaussian
fft_g = np.absolute(np.fft.fft(gaussian,nfft))
fft_g = fft_g[:np.floor(nfft/2)+1]
plt.plot(f,fft_g/np.max(fft_g))
plt.xlim([0,60])


#plot the convolved power spectrum
plt.subplot(428)
plt.plot(f,fft_s*fft_g)
plt.xlim([0,60])

plt.tight_layout()