python for autozero

#!/usr/bin/env python3
from matplotlib import pyplot as plt
import numpy as np
import scipy as sci

np.seterr(divide='ignore',invalid='ignore')

fTs=np.arange(-10,10,0.01)
pi=np.pi

h0=((1-sci.sin(2*pi*fTs)/(2*pi*fTs))**2 +((1-sci.cos(2*pi*fTs))/(2*pi*fTs))**2)**0.5
hn=((sci.sin(2*pi*fTs)/(2*pi*fTs))**2 +((1-sci.cos(2*pi*fTs))/(2*pi*fTs))**2)**0.5

def az_white(fTs,fcTs,shft_num):
    # fTs is k*fs, fc is -3dB bw, shft_num is k of k*fs
    noise_pwr = []
    base = fTs # just to initialize base
    for i in range(-shft_num,shft_num+1):
        fTs_i = fTs-i
        if i == 0:
            base = 1/(1+(fTs_i/fcTs)**2)*h0**2
            noise_pwr.append(base)
            plt.plot(fTs,base,'b--',label="Baseband Component")
        else:
            noise_pwr.append(1/(1+(fTs_i/fcTs)**2)*hn**2)

    total_pwr = sum(noise_pwr)
    fold_pwr = total_pwr - base
    plt.plot(fTs,fold_pwr,'r--',label="Foldover Component")
    plt.plot(fTs,total_pwr,color="red",label="AZ output PSD")
    plt.plot(fTs,[1]*len(list(fTs)),'k--',label="AZ input PSD")
    plt.grid(True)
    plt.legend(loc="best")
    plt.xlim(-1,1)
    plt.ylim(0,16)
    plt.xlabel("fTs")
    plt.ylabel("Normalized White Noise PSD of AZ")


def az_1_f(fTs,fcTs,fkTs,shft_num):
    # fTs is k*fs, fc is -3dB bw, fk is 1/f noise corner, shft_num is k of k*fs
    noise_pwr = []
    base = fTs # just to initialize base
    for i in range(-shft_num,shft_num+1):
        fTs_i = fTs-i
        if i == 0:
            base = (fkTs/abs(fTs_i))*1/(1+(fTs_i/fcTs)**2)*h0**2
            noise_pwr.append(base)
            plt.plot(fTs,base,'b--',label="Baseband Component")
        else:
            noise_pwr.append(fkTs/abs(fTs_i)*1/(1+(fTs_i/fcTs)**2)*hn**2)

    total_pwr = sum(noise_pwr)
    fold_pwr = total_pwr - base
    plt.plot(fTs,fold_pwr,'r--',label="Foldover Component")
    plt.plot(fTs,total_pwr,color="red",label="AZ output PSD")
    plt.plot(fTs,(fkTs/abs(fTs)*1/(1+(fTs/fcTs)**2)),'k--',label="AZ input PSD")
    plt.grid(True)
    plt.legend(loc="best")
    plt.xlim(-1,1)
    plt.ylim(0,16)
    plt.xlabel("fTs")
    plt.ylabel("Normalized 1/f Noise PSD of AZ")

plt.figure()
plt.subplot(1,2,1)
az_white(fTs,5,1000)
plt.subplot(1,2,2)
az_1_f(fTs,5,1,1000)
plt.show()