Merge pull request MTG#138 from dwsdolce/master

Python 3 lecture plot code changes

lectures/04-STFT/plots-code/spectrogram.py

@@ -19,15 +19,19 @@
 plt.figure(1, figsize=(9.5, 6))
 
 plt.subplot(211)
-numFrames = int(mX[:,0].size)
+numFrames = int(mX[:,0].size) + 1
+# Size of X and Y must be 1 larger than the size of mX for flat shading
 frmTime = H*np.arange(numFrames)/float(fs)                             
-binFreq = np.arange(N/2+1)*float(fs)/N                         
+binFreq = np.arange(N/2+1+1)*float(fs)/N
+binFreq2 = np.arange(mX[:,0].size)*float(fs)/N
 plt.pcolormesh(frmTime, binFreq, np.transpose(mX), shading = 'flat')
 plt.title('mX (piano.wav), Hamming window, M=1001, N=1024, H=256')
 plt.autoscale(tight=True)
 
 plt.subplot(212)
-numFrames = int(pX[:,0].size)
+# Size of X must be 1 larger than the size of np.diff(pX) (which has the Y axis size
+# reduced by 1) for flat shading
+numFrames = int(pX[:,0].size) + 1
 frmTime = H*np.arange(numFrames)/float(fs)                             
 binFreq = np.arange(N/2+1)*float(fs)/N                         
 plt.pcolormesh(frmTime, binFreq, np.diff(np.transpose(pX),axis=0), shading = 'flat')

lectures/04-STFT/plots-code/stft-system.py

@@ -23,15 +23,18 @@
 plt.axis([0,x.size/float(fs),min(x),max(x)])
 
 plt.subplot(412)
-numFrames = int(mX[:,0].size)
+# Size of X and Y must be 1 larger than the size of mX for flat shading
+numFrames = int(mX[:,0].size) + 1
 frmTime = H*np.arange(numFrames)/float(fs)                             
-binFreq = np.arange(mX[0,:].size)*float(fs)/N                         
+binFreq = np.arange(mX[0,:].size + 1)*float(fs)/N                         
 plt.pcolormesh(frmTime, binFreq, np.transpose(mX), shading = 'flat')
 plt.title('mX, Hamming window, M=1024, N=1024, H=512')
 plt.autoscale(tight=True)
 
 plt.subplot(413)
-numFrames = int(pX[:,0].size)
+# Size of X must be 1 larger than the size of np.diff(pX) (which has the Y axis size
+# reduced by 1) for flat shading
+numFrames = int(pX[:,0].size) + 1
 frmTime = H*np.arange(numFrames)/float(fs)                             
 binFreq = np.arange(pX[0,:].size)*float(fs)/N                         
 plt.pcolormesh(frmTime, binFreq, np.diff(np.transpose(pX),axis=0), shading = 'flat')

lectures/05-Sinusoidal-model/plots-code/hamming.py

@@ -1,6 +1,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from scipy.fftpack import fft
+eps = np.finfo(float).eps
 
 M = 64
 N = 512
@@ -18,7 +19,9 @@
 
 
 X = fft(fftbuffer)
-mX = 20*np.log10(abs(X)) 
+absX = abs(X)
+absX[absX < eps] = eps
+mX = 20*np.log10(absX) 
 mX1[:hN] = mX[hN:]
 mX1[N-hN:] = mX[:hN]      
 

lectures/05-Sinusoidal-model/plots-code/sine-analysis-synthesis.py

@@ -6,6 +6,7 @@
 sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../../software/models/'))
 import dftModel as DFT
 import utilFunctions as UF
+eps = np.finfo(float).eps
 
 (fs, x) = UF.wavread('../../../sounds/oboe-A4.wav')
 M = 601
@@ -21,13 +22,14 @@
 ploc = UF.peakDetection(mX, t)
 iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)
 freqs = iploc*fs/N 
-Y = UF.genSpecSines(freqs, ipmag, ipphase, Ns, fs)       
-mY = 20*np.log10(abs(Y[:hNs]))
+Y = UF.genSpecSines(freqs, ipmag, ipphase, Ns, fs)
+absY = abs(Y[:hNs])
+absY[absY < eps] = eps
+mY = 20*np.log10(absY)
 pY = np.unwrap(np.angle(Y[:hNs]))
 y= fftshift(ifft(Y))*sum(blackmanharris(Ns))
 
 plt.figure(1, figsize=(9, 6))
-
 plt.subplot(4,1,1)
 plt.plot(np.arange(-M/2,M/2), x1, 'b', lw=1.5)
 plt.axis([-M/2,M/2, min(x1), max(x1)])
@@ -44,9 +46,10 @@
 plt.axis([0, hNs,-90,max(mY)+2])
 plt.title("mY; Blackman-Harris; Ns = 512")
 
+yReal = np.real(y)
 plt.subplot(4,1,4)
-plt.plot(np.arange(Ns), y, 'b', lw=1.5)
-plt.axis([0, Ns,min(y),max(y)])
+plt.plot(np.arange(Ns), yReal, 'b', lw=1.5)
+plt.axis([0, Ns,min(yReal),max(yReal)])
 plt.title("y; Ns = 512")
 
 plt.tight_layout()

lectures/05-Sinusoidal-model/plots-code/sineModelAnal-bendir-phase.py

@@ -22,10 +22,10 @@
 
 plt.figure(1, figsize=(9.5, 7))
 maxplotfreq = 800.0
-maxplotbin = int(N*maxplotfreq/fs)
+maxplotbin = int(np.ceil(N*maxplotfreq/fs))
 numFrames = int(mX[:,0].size)
-frmTime = H*np.arange(numFrames)/float(fs)                             
-binFreq = np.arange(maxplotbin+1)*float(fs)/N 
+frmTime = H*np.arange(numFrames)/float(fs)                        
+binFreq = np.arange(maxplotbin)*float(fs)/N 
 plt.pcolormesh(frmTime, binFreq, np.transpose(np.diff(pX[:,:maxplotbin+1],axis=1)))                        
 plt.autoscale(tight=True)
 

lectures/05-Sinusoidal-model/plots-code/spec-sine-synthesis-lobe.py

@@ -4,6 +4,7 @@
 from scipy.fftpack import fft, ifft
 import math
 import sys, os, functools, time
+eps = np.finfo(float).eps
 
 sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../../software/models/'))
 
@@ -32,7 +33,8 @@
 
 powerX = sum(2*mX[0:hN]**2)/N
 
-mask = np.zeros(N//2)
+# add eps so that log10 does not error
+mask = np.zeros(N//2) + eps
 mask[int(N*f0/fs-2*N/float(M)):int(N*f0/fs+3*N/float(M))] = 1.0
 mY = mask*mX[0:hN]
 powerY = sum(2*mY[0:hN]**2)/N
@@ -55,13 +57,13 @@
 plt.axis([0,hN,-120,0])
 
 plt.subplot(3,2,3)
-plt.plot(y[0:M], 'b', lw=1.5)
+plt.plot(np.real(y[0:M]), 'b', lw=1.5)
 plt.axis([0,M,-1,1])
 plt.title ('y (synthesis of main lobe)')
 
 plt.subplot(3,2,5)
 yerror = xw - y
-plt.plot(yerror, 'k', lw=1.5)
+plt.plot(np.real(yerror), 'k', lw=1.5)
 plt.axis([0,M,-.003,.003])
 plt.title ("error function: x-y; SNR = ${%d}$ dB" %(SNR1))
 
@@ -75,7 +77,7 @@
 
 powerX = sum(2*mX[0:hN]**2)/N
 
-mask = np.zeros(N//2)
+mask = np.zeros(N//2) + eps
 mask[int(N*f0/fs-4*N/float(M)):int(N*f0/fs+5*N/float(M))] = 1.0
 mY = mask*mX[0:hN]
 powerY = sum(2*mY[0:hN]**2)/N
@@ -94,13 +96,13 @@
 plt.axis([0,hN,-120,0])
 
 plt.subplot(3,2,4)
-plt.plot(y[0:M], 'b', lw=1.5)
+plt.plot(np.real(y[0:M]), 'b', lw=1.5)
 plt.axis([0,M,-1,1])
 plt.title ('y (synthesis of main lobe)')
 
 plt.subplot(3,2,6)
 yerror2 = xw - y
-plt.plot(yerror2, 'k', lw=1.5)
+plt.plot(np.real(yerror2), 'k', lw=1.5)
 plt.axis([0,M,-.003,.003])
 plt.title ("error function: x-y; SNR = ${%d}$ dB" %(SNR2))
 

lectures/05-Sinusoidal-model/plots-code/spectral-sine-synthesis.py

@@ -4,6 +4,7 @@
 from scipy.fftpack import fft, ifft, fftshift
 import math
 import sys, os, functools, time
+eps = np.finfo(float).eps
 
 sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../../software/models/'))
 
@@ -23,7 +24,9 @@
 ypphase = phases
 
 Y = UF.genSpecSines(freqs, ypmag, ypphase, Ns, fs)       
-mY = 20*np.log10(abs(Y[:hNs]))
+absmY = abs(Y[:hNs])
+absmY[absmY < eps] = eps
+mY = 20*np.log10(absmY)
 pY = np.unwrap(np.angle(Y[:hNs]))
 y= fftshift(ifft(Y))*sum(blackmanharris(Ns))
 
@@ -40,8 +43,8 @@
 plt.title("pY, phases (radians) = .5, 1.2, 2.3")
 
 plt.subplot(3,1,3)
-plt.plot(np.arange(-hNs, hNs), y, 'b', lw=1.5)
-plt.axis([-hNs, hNs,min(y),max(y)])
+plt.plot(np.arange(-hNs, hNs), np.real(y), 'b', lw=1.5)
+plt.axis([-hNs, hNs,min(np.real(y)),max(np.real(y))])
 plt.title("y")
 
 plt.tight_layout()

lectures/05-Sinusoidal-model/plots-code/synthesis-window-2.py

@@ -6,6 +6,7 @@
 sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../../software/models/'))
 import dftModel as DFT
 import utilFunctions as UF
+eps = np.finfo(float).eps
 
 (fs, x) = UF.wavread('../../../sounds/oboe-A4.wav')
 M = 601
@@ -23,7 +24,9 @@
 iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)
 freqs = iploc*fs/N 
 Y = UF.genSpecSines(freqs, ipmag, ipphase, Ns, fs)       
-mY = 20*np.log10(abs(Y[:hNs]))
+absY = abs(Y[:hNs])
+absY[absY < eps] = eps
+mY = 20*np.log10(absY)
 pY = np.unwrap(np.angle(Y[:hNs]))
 y= fftshift(ifft(Y))*sum(blackmanharris(Ns))
 sw = np.zeros(Ns) 
@@ -37,9 +40,9 @@
 plt.figure(1, figsize=(9, 6))
 
 plt.subplot(3,1,1)
-plt.plot(np.arange(-hNs,hNs), y, 'b', lw=1.5)
-plt.plot(np.arange(-hNs,hNs), max(y)*bh/max(bh), 'k', alpha=.5, lw=1.5)
-plt.axis([-hNs, hNs,min(y),max(y)+.1])
+plt.plot(np.arange(-hNs,hNs), np.real(y), 'b', lw=1.5)
+plt.plot(np.arange(-hNs,hNs), max(np.real(y))*bh/max(bh), 'k', alpha=.5, lw=1.5)
+plt.axis([-hNs, hNs,min(np.real(y)),max(np.real(y))+.1])
 plt.title("y; size = Ns = 512 (Blackman-Harris)")
 
 plt.subplot(3,3,4)
@@ -57,10 +60,10 @@
 plt.axis([-hNs, hNs,0,1])
 plt.title("triangular / Blackman-Harris")
 
-yw = y * sw / max(sw)
+yw = np.real(y * sw / max(sw))
 plt.subplot(3,1,3)
 plt.plot(np.arange(-hNs,hNs), yw, 'b', lw=1.5)
-plt.plot(np.arange(-hNs/2,hNs/2), max(y)*ow/max(ow), 'k', alpha=.5, lw=1.5)
+plt.plot(np.arange(-hNs/2,hNs/2), max(np.real(y))*ow/max(ow), 'k', alpha=.5, lw=1.5)
 plt.axis([-hNs, hNs,min(yw),max(yw)+.1])
 plt.title("yw = y * triangular / Blackman Harris; size = Ns/2 = 256")
 

lectures/05-Sinusoidal-model/plots-code/synthesis-window.py

@@ -6,6 +6,7 @@
 sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), '../../../software/models/'))
 import dftModel as DFT
 import utilFunctions as UF
+eps = np.finfo(float).eps
 
 (fs, x) = UF.wavread('../../../sounds/oboe-A4.wav')
 M = 601
@@ -23,9 +24,11 @@
 iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc)
 freqs = iploc*fs/N 
 Y = UF.genSpecSines(freqs, ipmag, ipphase, Ns, fs)       
-mY = 20*np.log10(abs(Y[:hNs]))
+absY = abs(Y[:hNs])
+absY[absY < eps] = eps
+mY = 20*np.log10(absY)
 pY = np.unwrap(np.angle(Y[:hNs]))
-y= fftshift(ifft(Y))*sum(blackmanharris(Ns))
+y= np.real(fftshift(ifft(Y))*sum(blackmanharris(Ns)))
 sw = np.zeros(Ns) 
 ow = triang(2*H);    
 sw[hNs-H:hNs+H] = ow  

lectures/05-Sinusoidal-model/plots-code/two-sines-spectrum-blackman.py

@@ -10,7 +10,7 @@
 (fs, x) = UF.wavread('../../../sounds/sine-440-490.wav')
 w = np.blackman(5291)
 N = 16384
-pin = .11*fs
+pin = int(.11*fs)
 x1 = x[pin:pin+w.size]
 mX, pX = DF.dftAnal(x1, w, N)
 
@@ -27,7 +27,7 @@
 
 plt.subplot(3,1,3)
 plt.plot(fs*np.arange(mX.size)/float(N), pX, 'c', lw=1.5)
-plt.axis([100,900,-1,18])
+plt.axis([100,900,7,32])
 plt.title ('pX')
 
 plt.tight_layout()

notebooks/E1-Python-and-sounds.ipynb

@@ -24,7 +24,9 @@
    "metadata": {},
    "source": [
     "## Part 1 - Reading in an audio file\n",
-    "The `read_audio_samples()` function bellow should read an audio file and return a specified number of consecutive samples of the file starting at a given sample. \n",
+    "Write a function that reads an audio file and returns 10 consecutive samples of the file starting from \n",
+    "the 50001th sample. This means that the output should exactly contain the 50001th sample to the 50010th \n",
+    "sample (10 samples). \n",
     "\n",
     "The input to the function is the file name (including the path), plus the location of first sample and the number of consecutive samples to take, and the output should be a numpy array.\n",
     "\n",
@@ -57,7 +59,7 @@
    "source": [
     "# E1 - 1.1: Complete the read_audio_samples() function\n",
     "\n",
-    "def read_audio_samples(input_file, first_sample=50001, num_samples=10):\n",
+    "def readAudio(input_file, first_sample=50001, num_samples=10):\n",
     "    \"\"\"Read num_samples samples from an audio file starting at sample first_sample\n",
     "    \n",
     "    Args:\n",
@@ -75,7 +77,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "You can use as input the sound files from the sounds directory, thus using a relative path to it. If you run the `read_audio_samples()` function using the `piano.wav` sound file as input, with the default arguments, it should return the following samples:\n",
+    "You can use as input the sound files from the sounds directory, thus using a relative path to it. If you run the `readAudio()` function using the `piano.wav` sound file as input, with the default arguments, it should return the following samples:\n",
     "```\n",
     "array([-0.06213569, -0.04541154, -0.02734458, -0.0093997, 0.00769066, 0.02319407, 0.03503525, 0.04309214, 0.04626606,  0.0441908], dtype=float32)\n",
     "```"
@@ -299,7 +301,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -313,7 +315,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.10"
+   "version": "3.10.4"
   }
  },
  "nbformat": 4,