Using Python Numpy FFT

[steve_bank]

Familiarizing myself with Python Numpy.

Numpy FFT.

fft() returns a mirror spectrum. for a signal amplitude of 1 there are two lines about zero frequency amplitude .5. Displaying one half the spectrum, 500 Hz in the example, multiply amplitude by 2.

Some fft routines scale the outputs, some do not. The DFT is a correlation, an integral sum for each frequency. To scale the transform to the input amplitude divide by n.

The fft output is complex, Magnitude is the absolute value. Phase arctan(imag/real).

Feed the complex output of fft() into ifft() and you get the original signal back.

the fft is used to clean up audio signals. For example filtering out a 60Hz interference signal. Zero out the offending signal in the frequency domain and take the inverse transform to get the filtered time domain signal..

Frequency and time are linked. Period of a sine in seconds is t = 1/frequency. The sampling frequency from the Nyquist Theorem must be greater than twice the highest input frequency to avoid aliasing. At a sample frequency of 1000 Hz the unaliased frequency range is < 500 Hz.

You can see aliasing by increasing the signal frequency beyond 500 Hz out past 1000 Hz.

The transform is blind. Numbers in and numbers out. Meaning is from how the inputs are scaled. For a sampling frequency of 1000 Hz the signal sis ampled every .001 seconds. That translates to a frequency resolution of (sampling frequency)/(number transform points). Increasing transform size increases frequency resolution. Which is why fast algorithms are important.

Zero frequency is the average value of the signal.

import numpy as np
import warnings
import matplotlib.pyplot as plt
#----------------------------------------------------
def plotxy(xstart,xstop,x,y):
        [fig, p1] = plt.subplots(1)
        p1.plot(x,y,linewidth=2.0,color="k")
        p1.grid(which='major', color='k',linestyle='-', linewidth=1)
        p1.set_xlim(xstart,xstop)
        plt.show()
#--------------------------------------------------------------------------------
def plot2x(xlo,xhi,x,y1,y2,title1,title2,xlabel1,xlabel2):
        font1 = {'family': 'arial',
        'color':  'black',
        'weight': 'heavy',
        'size': 10,
        }
        [fig, p] = plt.subplots(1,2)
        p[0].set_xlim(xlo,xhi)
        p[0].plot(x, y1,linewidth=2.0,color="k")
        p[0].grid(color='k', linestyle='-', linewidth=1)
        p[0].set_xlabel(xlabel1)
        p[0].set_title(title1, fontdict = font1)

        p[1].set_xlim(xlo,xhi)
        p[1].plot(x, y2,linewidth=2.0,color="k")
        p[1].grid(color='k', linestyle='-', linewidth=1)
        p[1].set_xlabel(xlabel2)
        p[1].set_title(title2, fontdict = font1)
        plt.show()
#--------------------------------------------------------------------------------
def plotfft(reqs,mag,phase,fl0,fhi):
        font1 = {'family': 'arial',
        'color':  'black',
        'weight': 'heavy',
        'size': 10,
        }
        [fig, p] = plt.subplots(1,2)

        p[0].set_xlim(flo,fhi)
        p[0].plot(freqs, mag,linewidth=2.0,color="k")
        p[0].grid(color='k', linestyle='-', linewidth=1)
        p[0].set_xlabel("Frenquency Hertz", fontdict = font1)
        p[0].set_title("Magnitude", fontdict = font1)
          
        p[1].set_xlim(flo,fhi)
        p[1].set_ylim(-90,90)
        p[1].plot(freqs, phase,linewidth=2.0,color="k")
        p[1].grid(color='k', linestyle='-', linewidth=1)
        p[1].set_ylabel("Degrees", fontdict = font1)
        p[1].set_xlabel("Frenquency Hertz", fontdict = font1)
        p[1].set_title("Phase", fontdict = font1)
        plt.show()
#----------------------------------------------------------------------------------
        
def make_sin(f1,f2,a1,a2,os1,os2,phi1,phi2 ,y,t):
    n = len(y)
    for i in range(n):
            y1 = a1*np.sin(2*np.pi*f1*t[i]+phi1)+os1
            y2 = a2*np.sin(2*np.pi*f2*t[i]+phi2)+os2
            y[i] = y1 + y2
#-----------------------------------------------------------------------------------
def make_time_freq(fs,fr,t):
    n = len(t)
    dt = 1/fs
    df = fs/float(n)
    freq = 1.
    for i in range(n):
        t[i] = i*dt
        fr[i] = freq
        freq += df
#---------------------------------------------------------------------------------
n = pow(2,16)
y = np.ndarray(shape=(n),dtype=np.double)
freqs = np.ndarray(shape=(n),dtype=np.double)
yc = np.ndarray(shape=(n),dtype=np.cdouble)
ym = np.ndarray(shape=(n),dtype=np.double)
yp = np.ndarray(shape=(n),dtype=np.double)
t = np.ndarray(shape=(n),dtype=np.double)

fs = 1000.
f1 = 60
f2 = 200
a1 = 1
a2 = 1
os1 = 0
os2 = 0
phi1 = 0
phi2 = 0
make_time_freq(fs,freqs,t)
make_sin(f1,f2,a1,a2,os1,os2,phi1,phi2 ,y,t)
ytf = np.fft.fft(y)/n  # scale by n

for i in range(n):
        if(i==0):ym[i] = ytf[i].real  # first bin is average value
        else:ym[i] = abs(ytf[i])*2
        yp[i] = np.arctan(ytf[i].imag/ytf[i].real)*180/np.pi

yinv = np.fft.ifft(ytf*n)   # frequncy domain to time domain                                                                         

flo = 0  #Hertz
fhi = 500  #make 1000 to see full spectrum
tlo = 0. # seconds
thi=  .05
plotfft(freqs,ym,yp,flo,fhi)
plot2x(tlo,thi,t,y,yinv,"Original Sgnal","Recovered Signal","secomds","seconds")

#zero out 60Hz
for i in range(n):
        if(freqs[i] < 100   or freqs[i] > 900):
                ytf[i] = 0 + 0.j

# reconstruct signal
yinv = np.fft.ifft(ytf*n)   # frequncy domain to time domain     
plotfft(freqs,abs(ytf),yp,0,500)           
plot2x(tlo,thi,t,y,yinv,"Original Sgnal","Filtered Signal","secomds","seconds")

 

[steve_bank]

One more example. Picking out a signal corrupted with noise. Zero out any frequency in the spectrum to zero. In a real system there could be a dynamic algorithm to minimize noise and maximize signal to noise ratio. The FFT is everywhere.

A little cleaner code.

import numpy as np
import warnings
import matplotlib.pyplot as plt

#----------------------------------------------------
def plotxy(xstart,xstop,x,y,title,xlabel):
        [fig, p1] = plt.subplots(1)
        p1.plot(x,y,linewidth=2.0,color="k")
        p1.grid(which='major', color='k',linestyle='-', linewidth=1)
        p1.set_xlim(xstart,xstop)
        p1.set_title(title)
        p1.set_xlabel(xlabel)
        plt.show()
#--------------------------------------------------------------------------------
def plot2x(xlo,xhi,x,y1,y2,title1,title2,xlabel1,xlabel2):
        font1 = {'family': 'arial',
        'color':  'black',
        'weight': 'heavy',
        'size': 10,
        }
        [fig, p] = plt.subplots(1,2)
        p[0].set_xlim(xlo,xhi)
        p[0].plot(x, y1,linewidth=2.0,color="k")
        p[0].grid(color='k', linestyle='-', linewidth=1)
        p[0].set_xlabel(xlabel1)
        p[0].set_title(title1, fontdict = font1)

        p[1].set_xlim(xlo,xhi)
        p[1].plot(x, y2,linewidth=2.0,color="k")
        p[1].grid(color='k', linestyle='-', linewidth=1)
        p[1].set_xlabel(xlabel2)
        p[1].set_title(title2, fontdict = font1)
        plt.show()
#--------------------------------------------------------------------------------
def plotfft(reqs,mag,phase,fl0,fhi):
        font1 = {'family': 'arial',
        'color':  'black',
        'weight': 'heavy',
        'size': 10,
        }
        [fig, p] = plt.subplots(1,2)

        p[0].set_xlim(flo,fhi)
        p[0].plot(freqs, mag,linewidth=2.0,color="k")
        p[0].grid(color='k', linestyle='-', linewidth=1)
        p[0].set_xlabel("Frenquency Hertz", fontdict = font1)
        p[0].set_title("Magnitude", fontdict = font1)
         
        p[1].set_xlim(flo,fhi)
        p[1].set_ylim(-90,90)
        p[1].plot(freqs, phase,linewidth=2.0,color="k")
        p[1].grid(color='k', linestyle='-', linewidth=1)
        p[1].set_ylabel("Degrees", fontdict = font1)
        p[1].set_xlabel("Frenquency Hertz", fontdict = font1)
        p[1].set_title("Phase", fontdict = font1)
        plt.show()
#----------------------------------------------------------------------------------
def make_sin(f1,f2,a1,a2,os1,os2,phi1,phi2 ,t):
   #sigbal plus noise
    n = len(t)
    y = np.ndarray(shape=(n),dtype=np.double)
    noise = np.ndarray(shape=(n),dtype=np.double)
    for i in range(n):
            y1 = a1*np.sin(2*np.pi*f1*t[i]+phi1)+os1
            y2 = a2*np.sin(2*np.pi*f2*t[i]+phi2)+os2
            noise[i] = np.random.random_integers(-2,1) + np.random.random()
            y[i] = y1 + y2 + noise[i]
    return y,noise
#-----------------------------------------------------------------------------------
def make_time_freq(fs,fr,t):
    n = len(t)
    dt = 1/fs
    df = fs/float(n)
    freq = 0.
    for i in range(n):
        t[i] = i*dt
        fr[i] = freq
        freq += df
#---------------------------------------------------------------------------------
def mag_phase(y):
        #manitude and phase from complex y
        n = len(y)
        m = np.ndarray(shape=(n),dtype=np.double)
        p = np.ndarray(shape=(n),dtype=np.double)
        for i in range(n):
                if(i==0):m[i] = y[i].real  # first bin is average value
                else:m[i] = abs(y[i])*2
                p[i] = np.arctan(y[i].imag/y[i].real)*180/np.pi
        return m,p
#-------------------------------------------------------------------------------
def noise_filter(y,threshold):
        n = len(y)
        for i in range(n):
                if(abs(y[i]) < threshold):y[i] = 0. +0.j
#------------------------------------------------------------------------------
n = pow(2,16)
fs = 2000.
f1 = 60
f2 = 200
a1 = 1
a2 = 0
os1 = 0
os2 = 0
phi1 = 0
phi2 = 0
dt = 1./fs
tmax = n*dt
t = np.linspace(0,tmax,n)
ys,noise = make_sin(f1,f2,a1,a2,os1,os2,phi1,phi2 ,t)
freqs = np.fft.fftfreq(n, d=dt)
ytf = np.fft.fft(ys)/n
yinv = np.fft.ifft(ytf*n)   # frequncy domain to time domain                                                                        
mag,phase = mag_phase(ytf)
yinv = np.fft.ifft(ytf*n)
flo = 0  #Hertz
fhi =500
tlo = 0. # seconds
thi=  .05
plotxy(0,.05,t,noise,"Noise","Seconds")
plotfft(freqs,mag,phase,flo,fhi)
plot2x(tlo,thi,t,ys,yinv,"Original Sgnal","Recovered Signal","secomds","seconds")

noise_filter(ytf,.3)
               
yinv = np.fft.ifft(ytf*n)
ytf = np.fft.fft(yinv)/n  
plotfft(freqs,abs(ytf),phase,0,500)          
plot2x(tlo,thi,t,ys,yinv,"Original Sgnal","Filtered Signal","secomds","seconds

 

[steve_bank]

Responding to Jaryn's question on transform of a step function.

The true unit step is theoretical, nothing changes in zero time. It is important in electronics. A fast rise time digital signal generates radiated electromagnetic noise, and the frequency spectrum of the noise is the transform of the step. Put an antenna near the signal attached to a spectrum analyzer and you will see the spectrum.

import numpy as np
import matplotlib.pyplot as plt

#------------------------------------------------------------------------------------------------
def fft_scale(fs,t,f):
        n = len(t)
        df = fs/n
        dt = 1./fs
        freq = 0
        for i in range(n):
                f[i] = freq
                t[i] = i * dt
                freq += df
#-----------------------------------------------------------------------
def make_pulse(tlo,thi,vlo,vhi,t,y):
    n = len(y)
    for i in range(n):
        if t[i] >= tlo and t[i] <= thi:
            y[i] = vhi
        else:
            y[i] = vlo
#------------------------------------------------------------
def plotxy(xlo,xhi,ylo,yhi,x,y,title,xlabel,ylabel):
        # single y plot
        font1 = {'family': 'arial',
        'color':  'black',
        'weight': 'heavy',
        'size': 15,
        }
        [fig, p1] = plt.subplots(1)
        p1.plot(x,y,linewidth=2.0,color="k")
        p1.grid(color='k', linestyle='-', linewidth=1)
        p1.grid(which='major', color='k',linestyle='-', linewidth=0.8)
        p1.grid(which='minor', color='k', linestyle='-', linewidth=0.3)
        p1.set_xlim(xlo,xhi)
        p1.set_ylim(ylo,yhi)          
        p1.set_title(title, fontdict = font1)
        p1.set_xlabel(xlabel, fontdict = font1)
        p1.set_ylabel(ylabel, fontdict = font1)
        p1.minorticks_on()
        plt.show()
#-------------------------------------------------------------------------------  
n = pow(2,16)
y = np.ndarray(shape=(n),dtype=np.double)
f = np.ndarray(shape=(n),dtype=np.double)
yfft = np.ndarray(shape=(n),dtype=np.cdouble)
t = np.ndarray(shape=(n),dtype=np.double)
mag = np.ndarray(shape=(n),dtype=np.double)

fs = 1000 #sampling frequency
fft_scale(fs,t,f)   # time and frequency x axis                  
make_pulse(1,2,0,10,t,y)
yfft = np.fft.fft(y)/n
mag = abs(yfft)  #yfft complex, take abs for magnitude
plotxy(0,2,-1,11,t,y,"Signal","time","")
plotxy(0,40,0,max(mag),f,mag,"FFT","Frequency","")

.

 

[steve_bank]

And the transform of step function is sin(x)/x

n = 1000
y = np.ndarray(shape=(n),dtype = np.double)
t = np.linspace(0,20*np.pi,n)  #radians
for i in range(n):y[i] = np.sin(t[i])/t[i]
y = abs(y)
plotxy(0,max(t),-1,1,t,y,"","","")