dos33fsprogs/graphics/gr/plasmagoria/sines/sine_fft.py
2023-09-03 23:32:26 -04:00

52 lines
3.9 KiB
Python

# https://stackoverflow.com/questions/59725933/plot-fft-as-a-set-of-sine-waves-in-python
import matplotlib.pyplot as plt
import numpy as np
import cmath
#sin1
#sin3 = (46,48,50,52,53,54,56,58,60,60,62,64,65,66,68,69,71,71,73,74,75,76,77,78,79,80,81,82,83,83,84,84,85,85,86,87,87,88,88,87,88,88,88,88,88,88,88,88,88,87,87,87,86,86,85,84,85,84,83,82,82,81,80,79,78,78,77,76,75,75,74,73,72,71,70,69,69,68,66,66,65,65,63,63,61,61,60,59,59,57,57,57,56,56,55,54,54,53,53,52,52,51,50,50,50,49,49,49,48,49,48,48,48,48,47,47,48,47,47,47,47,47,47,47,46,47,47,47,46,47,47,47,47,46,47,47,47,46,47,47,46,46,47,46,46,45,46,45,45,45,44,44,44,43,43,43,42,42,41,40,40,39,39,38,38,37,37,35,35,34,33,33,32,31,31,29,29,28,27,26,25,25,23,22,22,21,20,19,19,18,17,16,15,15,14,13,12,12,11,10,9,9,8,8,8,7,6,7,6,6,6,6,5,6,5,5,6,5,6,6,7,7,8,8,9,9,10,11,11,12,12,13,15,15,16,18,18,20,21,22,23,25,26,27,29,30,32,33,34,36,38,39,40,42,44,46)
#sin2
#sin3 = (46,51,56,60,64,67,71,75,78,81,84,86,89,90,92,93,93,94,94,93,92,90,89,87,85,83,79,76,73,70,66,62,58,54,50,46,42,38,35,31,28,24,21,18,16,14,12,10,9,8,7,7,7,7,9,10,11,13,15,17,19,22,25,28,31,34,38,41,44,47,50,54,56,59,62,63,66,68,70,71,72,73,75,75,75,74,74,73,73,72,70,68,67,65,63,60,58,56,53,51,48,46,44,42,40,38,36,34,33,32,31,31,30,30,29,29,30,30,31,32,33,34,36,37,39,41,43,45,46,48,51,53,55,56,58,60,61,62,63,63,64,64,65,64,64,63,63,62,61,59,58,56,54,52,49,47,45,43,41,37,35,33,31,29,27,25,24,22,21,20,20,19,19,19,19,20,22,23,24,26,28,29,32,35,38,40,44,46,50,53,56,59,62,65,69,72,75,76,79,81,83,84,85,85,87,87,87,86,85,83,82,80,78,75,73,69,66,63,59,55,52,48,44,39,35,31,28,24,20,17,14,11,9,7,5,3,2,1,0,0,1,1,2,3,5,7,10,13,16,19,23,26,30,34,38,42)
# sin3
#sin3 = (38,44,49,53,57,61,64,66,68,69,69,70,69,67,66,64,60,58,56,54,51,49,48,47,47,46,47,47,48,51,51,54,55,58,60,60,62,62,61,61,59,57,54,52,48,43,40,35,29,25,20,17,12,9,7,4,3,3,3,3,4,7,9,12,15,19,22,24,27,30,32,34,34,35,36,36,35,34,33,32,29,28,27,26,25,25,25,26,28,30,32,35,39,43,47,51,55,61,64,68,71,74,76,77,78,78,77,76,74,71,69,65,60,57,53,50,46,43,40,38,37,35,35,34,34,36,36,37,38,41,42,42,43,44,43,43,41,40,37,35,32,28,25,21,16,13,9,7,4,2,1,0,0,0,2,3,6,10,13,17,21,27,31,35,39,43,45,48,50,51,52,53,53,51,51,50,48,46,45,44,43,42,42,42,43,44,46,48,50,54,56,59,62,66,69,71,73,75,75,75,74,73,71,69,66,61,58,53,48,43,38,34,30,26,23,20,19,17,16,16,16,18,18,20,21,24,26,27,29,30,31,31,31,31,30,29,27,24,22,20,16,14,12,11,9,8,8,9,10,12,14,17,20,25,29,34)
# sin4
sin3 = (64,65,67,68,70,71,73,74,76,77,79,80,82,83,85,86,88,89,91,92,93,95,96,97,99,100,101,103,104,105,106,107,108,109,111,112,113,114,115,115,116,117,118,119,120,120,121,122,122,123,123,124,124,125,125,125,126,126,126,127,127,127,127,127,127,127,127,127,127,127,126,126,126,125,125,125,124,124,123,123,122,122,121,120,120,119,118,117,116,115,115,114,113,112,111,109,108,107,106,105,104,103,101,100,99,97,96,95,93,92,91,89,88,86,85,83,82,80,79,77,76,74,73,71,70,68,67,65,63,62,60,59,57,56,54,53,51,50,48,47,45,44,42,41,39,38,36,35,34,32,31,30,28,27,26,24,23,22,21,20,19,18,16,15,14,13,12,12,11,10,9,8,7,7,6,5,5,4,4,3,3,2,2,2,1,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,2,2,2,3,3,4,4,5,5,6,7,7,8,9,10,11,12,12,13,14,15,16,18,19,20,21,22,23,24,26,27,28,30,31,32,34,35,36,38,39,41,42,44,45,47,48,50,51,53,54,56,57,59,60,62)
#x = np.arange(0,256,1)
# discrete works better if periodicity matches range (?)
x = np.linspace(0,4*np.pi, 256)
#fft3=np.fft.fft(sin3)
#
#plt.plot(x, sin3)
#plt.show()
fft3 = np.fft.fft(sin3)
freqs = np.fft.fftfreq(len(x),.01)
threshold = 0.0
recomb = np.zeros((len(x),))
middle = len(x)//2 + 1
for i in range(middle):
if abs(fft3[i])/(len(x)) > threshold:
if i == 0:
coeff = 2
else:
coeff = 1
sinusoid = 1/(len(x)*coeff/2)*(abs(fft3[i])*np.cos(freqs[i]*2*np.pi*x+cmath.phase(fft3[i])))
recomb += sinusoid
print(abs(fft3[i]))
plt.plot(x,sinusoid)
#plt.show()
plt.plot(x,recomb,x,sin3)
plt.show()
#plt.plot(fft3.imag)
#plt.show()
#for i in range(1,256):
# plt.plot(x, fft3.imag[i] * np.sin(i*x)/100)
#plt.show()