Scientific
############### SOUND ##############
from scipy.io.wavfile import read
from scipy.io.wavfile import write
(fs,x) = read("../sms-tools/sounds/piano.wav")
print "Sample rate: " + str(fs)
print x
!play
#np.write(filename,fs,array)
############### MATH (and COMPLEX) ##############
np.conjugate(2+j)
sum([1,2,3,4,5])
sum([1+1j,2+1j,3+2j,4-1j,5])
#np.max(array)
############### Sinusoids ##############
# x[n] = A cos(2πf nT + φ)
# where T is the time between samples, equals 1/fs
# and n is the sample number
# --> x[n] = A cos(n*2πf/fs + φ)
# x[n] = A e j(ωnT + φ) = A cos(ωnT + φ) + jA sin(ωnT + φ)
#
# x[n] = e j2πkn/N = cos(2πkn/N) + jsin(2πkn/N)
# where N is DFT size
############### A2. Sinusoids ##############
############### Real Sinusoid ##############
import numpy as np
import matplotlib.pyplot as plt
A = .8 # amp
f0 = 1000 # freq
phi = np.pi/2 # starting phase
fs = 44100
# define an array with a range of values
t = np.arange(-.002, .002, 1.0/fs)
# create an array of sinusoid values
x = A * np.cos(t*2*np.pi*f0 + phi)
plt.plot(t,x)
plt.axis([-.002,.002,-.8,.8])
plt.xlabel("time")
plt.ylabel("amplitude")
plt.show()
############### Complex Sinusoid ##############
N = 32 # DFT size ???
k = 5 # frequency
n = np.arange(-N/2,N/2) # 'time' index a.k.a. sample number
s = np.exp(1j * 2 * np.pi * k * n / N)
plt.plot(n, np.real(s))
plt.axis([-N/2,N/2-1,-1,1])
plt.xlabel("n")
plt.ylabel("amplitude")
plt.show()
############### DFT ##############
# X[k] = SIGMA {n=0; N-1} x[n] e -j2πkn/N
# where
# k is frequency index in range 0..N-1
# n is sample number of time-domain signal
X = np.array([])
N=64
k0 = 7
x = np.exp(1j*2*np.pi*k0/N*np.arange(N))
# NB: range does not create an array, like np.arange does
# loop over all frequencies, which gives the spectral component
# for every frequency index
for k in range(N):
# inner product (thus complex conjugate!) of original signal
# with complex sinusoid with frequency k
# create an array of samples of the complex sinusoid
# I would create an incremental sum in an inner loop,
# but hey whatever...
s = np.exp(1j * 2 * np.pi * k * np.arange(N) / N)
#X = sum(x*np.conjugate(s))
# append value to array X, so we don't need k as an index into X
X = np.append(X, sum(x*np.conjugate(s)))
plt.plot(np.arange(N), np.real(x))
plt.axis([0,N-1,-1,1])
plt.xlabel("n")
plt.ylabel("amplitude")
plt.show()
plt.plot(np.arange(N), abs(X))
plt.axis([0,N-1,0,N])
plt.xlabel("n")
plt.ylabel("frequency")
plt.show()