Cacher is the code snippet organizer for pro developers

import numpy as np
import matplotlib.pyplot as plt
import struct
import wave
import math


def iexp(n):
    return complex(math.cos(n), math.sin(n))

def is_pow2(n):
    return False if n == 0 else (n == 1 or is_pow2(n >> 1))

def fft_(xs, n, start=0, stride=1):
    "cooley-turkey fft"
    if n == 1: return [xs[start]]
    hn, sd = n // 2, stride * 2
    rs = fft_(xs, hn, start, sd) + fft_(xs, hn, start + stride, sd)
    for i in range(hn):
        e = iexp(-2 * math.pi * i / n)
        rs[i], rs[i + hn] = rs[i] + e * rs[i + hn], rs[i] - e * rs[i + hn]
    return rs


# Digital Fourier Transform (離散傅立葉轉換)
# 離散フーリエ変換
# In summary, DFT calculates the digital form of Fourier transform
# which brings a signal represented in time-domain to frequency-domain
# where each frequency bin (axis) has a magnitude and phase coefficient
def dft (N, g):
    G = [0.0] * N
    # for each frequency bin
    for k in range(N):
        # is the summation of the evaluated n-th sample complex values! (or Sin/Cos)
        for n in range(N):
            # As compared to its continuious version of FT
            # k/N is F (frequency)
            # n is t (time)
            real = np.cos(2 * np.pi * k * n / N)
            imag = - np.sin(2 * np.pi * k * n / N)
            G[k] += g[n] * complex(real, imag)
    return G

# sine wave generator 
# --------------------
# amplitude, base frequency, num of samples, sampling duration
def create_wave(A,f0,fs,t,num):
    data = []
    for n in np.arange(t * fs): # arange(): return evently spaced values
        y = (A / num) * np.sin(2 * np.pi * f0 * num * n / fs)
        # clip the amplitude to [-1 1]
        if y > 1.0: y = 1.0
        if y < -1.0: y = -1.0
        
        # C++ equivalence 'data.push_back(y)'
        data.append(y)
    return data

def main():
    
    # Generate square wave from sine waves
    data = create_wave(1, 261.63, 10000, 1, 1)
    for i in range(3,301,2):
        data = np.add(data, create_wave(1, 261.63, 10000, 1, i)) 

    data = data/np.max(data,0)
    N = 256                                       
    #G = dft(n0, N, data) 
    G = fft_(data,N)
    amp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in G]      
    phase = [np.arctan2(int(c.imag), int(c.real)) for c in G]   
    flist = [k * 10000 / N for k in range(N)]
   
   
    plt.figure(0)
    plt.plot(range(0, N), data[0:N])
    plt.axis([0, N, -1.0, 1.0])
    plt.xlabel("Time [sample]")
    plt.ylabel("Amplitude")


    #plt.subplot(312)
    plt.figure(1)
    plt.plot(flist, amp, marker='o', linestyle='-')
    #plt.axis([0, 10000/2, 0, 15])
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Amplitude spectrum")


    #plt.subplot(313)
    plt.figure(2)
    plt.plot(flist, phase, marker='o', linestyle='-')
    plt.axis([0, 10000/2, -np.pi, np.pi])
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Phase spectrum")
    plt.show()

if __name__ == '__main__':
    main()