# 538 Day -3: Trigonometric interpolation and spectral differentiation¶

## Visualization of the matrix of the discrete Fourier transform¶

In [1]:
from numpy import *
from matplotlib import rcdefaults
#rcdefaults()  # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'
#import seaborn as sns  # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt  # use matplotlib functionality directly
%matplotlib inline
#sns.set()

In [2]:
from numpy import linspace,ones_like,real,imag,cos,sin,pi,exp,arange
import matplotlib.pyplot as plt
from matplotlib import cm

def fft_matrix_viz(n):
L = 1.
vsep = -1.
markersize = 75./float(n)
i = complex(0.,1.)
t = linspace(0,L,501,endpoint=True)
kvals = range(0,n)#-n//2,n)
colormap = cm.get_cmap('hsv') # a cyclic colormap
nc = len(kvals)
colors = colormap( linspace(0,1,nc,endpoint=False) )  # get nc colors from the colormap
# pastelize colors
p = 0.5
colors = (1-p)*colors + p*ones_like(colors)
# darken colors
d = 0.9
colors *= [d,d,d,1]
fontsize = 20*4/n
plt.figure(figsize=(12, 12))

theta = linspace(0,2*pi,100)
c = cos(theta)
s = sin(theta)
def plotZ(x,y,Z):

plt.subplot(1,1,1,aspect=1,frameon=False)
for k in kvals:
zd = exp(-2*pi*i*k*arange(n)/float(n))
for j,Z in enumerate(zd): plotZ(j,vsep*k,Z)
plt.text(-0.5,vsep*k,str(k),va='center',ha='right',fontsize=fontsize)

for j in range(n):
plt.text(j,vsep*(min(kvals)-0.45),str(j),ha='center',fontsize=fontsize)
plt.xlim(-.5,n+.5)
plt.ylim(vsep*(max(kvals)+1),vsep*(min(kvals)-1))
plt.xticks([])
plt.yticks([])

plt.savefig('dft_visualization_'+'n'+str(n).zfill(4)+'.svg',bbox_inches='tight')

fft_matrix_viz(4)

In [3]:
fft_matrix_viz(8)


## Small example of DFT and IDFT¶

In [7]:
import numpy as np
from scipy.fftpack import fft, ifft
x = [1,2,3,5]
y = fft(x)
y

Out[7]:
array([11.+0.j, -2.+3.j, -3.+0.j, -2.-3.j])
In [8]:
ifft(y)

Out[8]:
array([1.+0.j, 2.+0.j, 3.+0.j, 5.+0.j])
In [9]:
## the matrix?
fft(np.eye(4))

Out[9]:
array([[ 1.+0.j,  1.+0.j,  1.+0.j,  1.-0.j],
[ 1.+0.j,  0.-1.j, -1.+0.j,  0.+1.j],
[ 1.+0.j, -1.+0.j,  1.+0.j, -1.-0.j],
[ 1.+0.j,  0.+1.j, -1.+0.j,  0.-1.j]])

## Another visualization of the DFT matrix¶

In [10]:
n=16
def fft_matrix_viz2(n):
L = 1.
vsep = -4.
markersize = 150./float(n)
i = complex(0.,1.)
t = linspace(0,L,501,endpoint=True)
kvals = range(n)#-n//2,n)

plt.figure(figsize=(12, 12))

recolor,imcolor = '#4444ff','#ff4444'
plt.subplot(1,1,1)
for k in kvals:
zd = exp(-2*pi*i*k*arange(n)/float(n))
zc = exp(-2*pi*i*k*t)
plt.plot([0,n],[vsep*k,vsep*k],'#dddddd')
plt.plot(t*n,vsep*k+real(zc),recolor,alpha=0.25)
plt.plot(t*n,vsep*k+imag(zc),imcolor,alpha=0.25)
plt.plot(arange(n),vsep*k+real(zd),'.',color=recolor,markersize=markersize)
plt.plot(arange(n),vsep*k+imag(zd),'.',color=imcolor,markersize=markersize)
plt.text(-0.75,vsep*k,str(k),va='center')

plt.title('$F_{'+str(n)+'}$ : real part blue, imag part red')
plt.xlim(-.5,n+.5)
plt.ylim(vsep*(max(kvals)+1),vsep*(min(kvals)-1))
plt.xticks(range(n))
plt.yticks([])

plt.savefig('dft_visualization2_'+'n'+str(n).zfill(4)+'.png')

fft_matrix_viz2(4)

In [11]:
fft_matrix_viz2(8)


## Example of DFT of a bigger vector¶

The vector is a recording of a piano note piano_low_f.wav (link fixed)

In [12]:
from scipy.io import wavfile
from numpy import *

print('sampling rate',rate)
print('x.shape',x.shape)

n = x.shape[0]
d = n/float(rate)  # duration of recording
print('duration',d,'secs')
#assert(0)
import matplotlib.pyplot as plt
#plt.plot(x)

from scipy.fftpack import fft

y = fft(x)

plt.figure(figsize=(20,5))
# mark the nominal frequency and its harmonics (integer multiples)
f = 43.6535
for i in range(20):
#plt.plot([i*f,i*f],[0,1.e9],'r')
plt.axvline(i*f,color='r')

plt.plot( arange(n)/d , abs(y),'o-' )

plt.xlim(0,3000/d)
plt.xlabel('frequency (Hz) (proportional to index k)');
plt.ylabel('$|y_k|$');

sampling rate 44100
x.shape (162030,)
duration 3.6741496598639456 secs


## Trigonometric interpolation¶

### Visualization of IDFT matrix¶

In [13]:
n=16
def fft_matrix_viz2(n):
L = 1.
vsep = -4.
markersize = 150./float(n)
i = complex(0.,1.)
t = linspace(0,L,501,endpoint=True)
kvals = range(-n//2+1,n//2)

plt.figure(figsize=(12, 12))

recolor,imcolor = '#4444ff','#ff4444'
plt.subplot(1,1,1)
for k in kvals:
zd = exp(-2*pi*i*k*arange(n)/float(n))
zc = exp(-2*pi*i*k*t)
plt.plot([vsep*k,vsep*k],[0,n],'#dddddd')
plt.plot(vsep*k+real(zc),t*n,recolor,alpha=0.25)
plt.plot(vsep*k+imag(zc),t*n,imcolor,alpha=0.25)
plt.plot(vsep*k+real(zd),arange(n),'.',color=recolor,markersize=markersize)
plt.plot(vsep*k+imag(zd),arange(n),'.',color=imcolor,markersize=markersize)
plt.text(vsep*k,-0.75,str(k),ha='center')

plt.title('$F_{'+str(n)+'}$ : real part blue, imag part red')
plt.ylim(n+.5,-.5)
plt.xlim(vsep*(min(kvals)-1),vsep*(max(kvals)+1))
plt.yticks(range(n))
plt.xticks([])

plt.savefig('idft_visualization2_'+'n'+str(n).zfill(4)+'.png')

fft_matrix_viz2(6)