Fourier Transforms

• Goal In this lab we will be using Fourier transform methods to model functions that we will create using sum waves. The Fourier transform method is a mathematical tool that can decompose any function into its frequency, this allows us to transform it into a time vs frequency representation.
• How we will test it making a wave that is the sum of 3 cosine functions. We will then take Fourier transform of the summed wave and interpret the results.

analyzing the number of sunspots versus time using a Fourier transform.

determining the frequencies and harmonics that give different instruments their unique tone frequency is the number of cycles or waves per second. Period is the time needed to complete one cycle of the wave. -Part 1 We will be summing cosine functions (f(t)=Acos(ωt+ϕ)), 0<t<4 with 100 steps,: ω=[π,2π,4π],A=[0.5,1,0.75],ϕ=[π/2,π,3π/2]

t = np.linspace(0,4,100)
# set up the three cosine curves
cos1 = 0.5*np.cos(np.pi*t+np.pi/2)
cos2 = 1*np.cos(2*np.pi*t+np.pi)
cos3 = 0.75*np.cos(4*np.pi*t+3*np.pi/2)
sum = cos1 + cos2 + cos3
A = [0.5, 1, 0.75]
omega = [np.pi, 2*np.pi, 4*np.pi]
phi = [np.pi/2, np.pi, 3*np.pi/2]

plt.plot(t, cos1, label='cos1')
plt.plot(t, cos2, label='cos2')
plt.plot(t, cos3, label='cos3')
plt.plot(t, sum, label='sum')
plt.legend()

• Part 2 computing the Fourier transform of the sum wave above We will be using a function not used in pervious labs, np.fft.rfft which will returns the amplitudes of the fourier transform components of the input data, giving component frequencies that make up the input function.The code below will show the application of this code

# take the FFT of the summed wave by passing sumwave into np.fft.rfft
fft = np.fft.rfft(sum)
# save the output into the variable Y
Y =np.fft.rfft(sum)

# set power equal to the absolute value of Y
power = np.abs(Y)

# set n equal to the length of the time array
n=len(t)

# calculate the time step between data points

dt = t[1]-t[0]

# use rfftfreq to get the array of frequencies
freq = np.fft.rfftfreq(n, d = dt)
plt.plot(freq, power)

# create a figure showing power versus frequency
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
# label the axes
# add a grid using plt.grid()
plt.grid()
plt.show()

Now we will create the inverse using irfft code will be below

 feed Y into irfft and set the result equal to Yinv
Yinv = np.fft.irfft(Y)

# create a plot that shows Yinv vs t
plt.plot(t, Yinv, label='Yinv')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid()
plt.show()


# overplot sumwave vs t using a dashed line
plt.plot(t, sum, '--', label='sumwave')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid()
plt.legend()
plt.show()

• Part 3: Analyzing the Periodicity in Sunspots

• In this part we will be analyzing sunspot data starting from 1749(wow) the data we will import contains months followed by the number of sunspots The goal is to plot and estimate the sunspot cycle. The code will be below

# code to read in file
month, nspots = np.loadtxt('sunspots.txt',unpack=True)
# save the month as month and the number of sunspots as nspots
dt = month[1]-month[0]
plt.plot(month, nspots)
plt.xlabel('Month')
plt.ylabel('Number of Sunspots')
plt.grid()
plt.show()

• Now with Fourier transform

# copy code from part 2 to convert the output of np.fft.rfft to power
power = np.abs(nspots_fft)


# calculate the corresponding frequencies using np.fft.rfftfreq
freq = np.fft.rfftfreq(len(nspots), d = dt)

# calculate the period corresponding to each frequency
period = 1/freq
plt.plot(period, power)
plt.xlabel('Period (months)')
plt.ylabel('Power')
plt.grid()
plt.show()

• Part 4 creating a square wave In this part we will show that any function can be modeled using a series of waves but to create this square wave we will have to constantly adjust the amplitude and the frequency of the wave.The wave must also have the frequency and the amplitude of successive waves must decrease in a specific pattern to create the square.

T = 1# Period is 1 sec
a0 = 1 # no DC offset
w0 = 2*np.pi/T

# try adding 50 sine functions together
sumwave = a0*np.sin(w0*T)
# let the frequency increase each time, and
frequency = np.arange(1,100,2) * w0
# let the amplitude decrease as the frequency increases
amp = 1/np.arange(1,100,2)

for i in range(len(frequency)):
  # compute a new sine function
  newsine = amp[i]*np.sin(frequency[i]*t)

  # plot the sine versus time
  plt.plot(t, newsine)

  # add the new function to the sum
  sumwave += newsine

# plot the sumwave vs time
plt.plot(t, sumwave, color='k', lw=3)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Square Wave')
plt.show()

• part 5 applying it to a simple concept such as music In this part we will be applying all the concepts and theories used above on piano and trumpet data that we will be importing first by plotting the wave function then applying fourier transform method onto that wave to get a better idea of how the frequency of notes are when played.

# plot the waveforms
piano = np.loadtxt('piano.txt')
trumpet = np.loadtxt('trumpet.txt')
plt.subplots(2,1)
plt.subplot(2,1,1)
plt.plot(piano)
plt.title('Piano')
plt.subplot(2,1,2)
plt.plot(trumpet)
plt.title('Trumpet')
plt.tight_layout()
plt.show()

• Now applying Fourier code

!wget http://www-personal.umich.edu/~mejn/cp/data/piano.txt
piano = np.loadtxt('piano.txt')
!wget http://www-personal.umich.edu/~mejn/cp/data/trumpet.txt
trumpet = np.loadtxt('trumpet.txt')

# compute the Fourier transform of the piano data using rfft
fortran = np.fft.rfft(piano)
# save the results as Ypiano
Ypiano = np.fft.rfft(piano)
# set Ppiano to the absolute value of Ypiano
Ppiano = np.abs(Ypiano)
# compute the Fourier transform of the trumpet data using rfft
fortran2 = np.fft.rfft(trumpet)
# save the results as Ytrumpet
Ytrumpet = np.fft.rfft(trumpet)
# set Ptrumpet to the absolute value of Ytrumpet
Ptrumpet = np.abs(Ytrumpet)
# calculate the time array
timearr = np.arange(len(piano))/44100

# calculate dt
dt = timearr[1] - timearr[0]

# calculate the frequency array
freq = np.fft.rfftfreq(len(piano), d = dt)

# convert the output from rfft to a frequency spectrum
spectrum_piano = np.abs(Ypiano)
spectrum_trumpet = np.abs(Ytrumpet)

# plot power vs frequency for piano
plt.subplots(1,2)
plt.subplot(1,2,1)
plt.plot(freq, spectrum_piano)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.title('Piano')
# plot power vs frequency for trumpet
plt.subplot(1,2,2)
plt.plot(freq, spectrum_trumpet)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.title('Trumpet')
plt.tight_layout()
plt.show()

In the end you can see how the power varies as we change the frequency of each note played in the piano and in the trumpet