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):
        radius = 0.3
        plt.plot(x+radius*c, y+radius*s,'k',alpha=0.35)  # draw unit circle
        plt.plot([x,x+radius*real(Z)],[y,y+radius*imag(Z)],'k',alpha=0.35) # draw radius
        plt.plot(x+radius*real(Z),y+radius*imag(Z),'o',markersize=markersize,color=colors[(j*k)%len(kvals)]) # mark 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 *

rate,x = wavfile.read('piano_low_f.wav')
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)