steve_bank
Diabetic retinopathy and poor eyesight. Typos ...
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.
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.
Code:
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")