Lab 08 - Curlss/Computational_physics GitHub Wiki
The goal of this lab was to learn more about the Fourier transform. We also wanted to get better at using np.fft.rfft, np.fft.rfftfreq, and np.fft.irfft to compute, understand, and reconstruct signals so that we could use them to look at data and audio files and make graphs and charts.
Since this lab was based on the Fourier transform, we wanted to make some Fourier transforms by adding up the cosine and sin functions. We did this to get the Fourier transform of the total wave, which we then used to figure out what the results meant. We learned about the Fourier transform in this lab by adding up three cosine functions and then using FFT and inverse FFT to figure out its frequency components. We were able to show how different frequency components and phase relationships work together in the time domain in this lab. We were also able to show how to decompose a signal into its constituent parts, find its peaks, and rebuild it without any loss. We could also see how fourier spectra can be used to find pitch and sound, as well as the name of an instrument. Lastly, we used what we had learned to make scripts that could be used again and again to load data, make plots, do spectral analysis, and filter out noise. These scripts were based on what we had done in the lab before.
Lab
We had to add up the cosine functions in the first part of this lab so that we could use them in the fourier transform. To do this, we had to make a 100-point time array t that ran from 0 to 4 s and then set up three cosines with amplitudes A = [0.5, 1.0, 0.75], angular frequencies ω = [π, 2π, 4π], and phase shifts ϕ = [π/2, π, 3π/2].
# use linspace to set up time array called t
t = np.linspace(0,4,100)
w = [np.pi,2*np.pi,4*np.pi]
A = [0.5, 1, 0.75]
phi = [np.pi/2,np.pi,3*np.pi/2]
# set up the three cosine curves
cos1 = A[0]*np.cos(w[0]*t + phi[0])
cos2 = A[1]*np.cos(w[1]*t + phi[1])
cos3 = A[2]*np.cos(w[2]*t + phi[2])
# create your graph here
plt.figure(figsize=(10,6))
plt.subplots_adjust(hspace = .9, )
plt.subplot(4,1,1)
plt.plot(t,cos1,label='cos1')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Cosine 1 function')
plt.grid()
plt.subplot(4,1,2)
plt.plot(t,cos2,label='cos2')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Cosine 2 function')
plt.grid()
plt.subplot(4,1,3)
plt.plot(t,cos3,label='cos3')
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Cosine 3 function')
plt.grid()
sumwave = cos1 + cos2 + cos3
plt.subplot(4,1,4)
plt.plot(t,sumwave,label='sumwave')
plt.legend()
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Three cosine functions Sum')
plt.grid()
plt.show()
To find the Fourier transform of a set of summed waves, we used np.fft.rfft, which gives us the amplitudes of the fourier transform components of the input data along with the frequencies that make up the input function; np.fft.rfftfreq, which figures out the frequency array that goes with the fourier transform amplitudes; and np.fft.irfft, which finds the inverse fourier transform, which turns frequency data back into time data. This helped us figure out our power, and then we could plot power vs. frequency, which showed us which bands were most important in a signal.
# take the FFT of the summed wave by passing sumwave into np.fft.rfft
# save the output into the variable Y
Y = np.fft.rfft(sumwave)
# set power equal to the absolute value of Y
power = abs(Y) # the power
# 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)
# create a figure showing power versus frequency
plt.figure()
plt.plot(freq,power,label='power')
plt.legend()
# label the axes
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.title('Power vs Frequency')
# Add a grid using plt.grid()
plt.grid()
We looked at our sunspots data by getting sunspots.txt and putting columns into month and nspots. Once that was done, we could plot nspots against months to get a rough idea of a 10–12 year cycle. By using FFT on nspots, which calculated power vs. time, we were able to find a main peak that happened around 11 years. This proved that our calculations for the Fourier transform were correct.
# run np.fft.rfft on nspots
Yspots = np.fft.rfft(nspots)
# copy code from part 2 to convert the output of np.fft.rfft to power
# take the FFT of the summed wave by passing sumwave into np.fft.rfft
# save the output into the variable Y
# set power equal to the absolute value of Y
power = abs(Yspots) # the power
n = len(month)
dt = month[1] - month[0]
# calculate the corresponding frequencies using np.fft.rfftfreq
freq = np.fft.rfftfreq(n, d = dt)
# calculate the period corresponding to each frequency
period = 1/freq
# plot the power versus period rather than versus frequency
plt.figure()
plt.plot(period,power,label='power')
plt.legend()
plt.xlabel('Period')
plt.ylabel('Power')
plt.title('Power vs Period')
plt.grid()
<ipython-input-11-f6ffff7d6c1a>:14: RuntimeWarning: divide by zero encountered in divide
period = 1/freq
By adding up sine waves with different frequencies and amplitudes, we made a square wave.
T = 1# Period is 1 sec
a0 = 1 # no DC offset
w0 = 2*np.pi/T
# try adding 50 sine functions together
# let the frequency increase each time, and
# let the amplitude decrease as the frequency increases
sumwave = np.zeros_like(t)
for i in range(1,100,2):
# compute a new sine function
sine = (1/i)*np.sin(i*w0*t)
# plot the sine versus time
plt.plot(t,sine)
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Sine function')
plt.grid()
# add the new function to the sum
sumwave += sine
# plot the sumwave vs time
plt.plot(t,sumwave,'k',label='sumwave')
plt.legend()
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Square Wave')
plt.grid()
Finally it was all about figuring out how different musical instruments sound. Our code could download and read instrument files, as well as load them into arrays. This let us do 44,100 Hz sampling for FFTs, power spectra, and frequency arrays. We were able to plot power vs. frequency next to each other and see that both instruments have the same basic frequency peak. The trumpet has clear, strong, evenly spaced harmonics, while the piano has a richer set of overtones that decay more quickly.
# compute the Fourier transform of the piano data using rfft
# save the results as Ypiano
Ypiano = np.fft.rfft(piano)
# set Ppiano to the absolute value of Ypiano
Ppiano = abs(Ypiano)
# compute the Fourier transform of the trumpet data using rfft
# save the results as Ytrumpet
Ytrumpet = np.fft.rfft(trumpet)
# set Ptrumpet to the absolute value of Ytrumpet
Ptrumpet = abs(Ytrumpet)
# calculate the time array
timearr = np.arange(len(piano))/44100
# calculate dt
dt = timearr[1] - timearr[0]
# convert the output from rfft to a frequency spectrum
freq = np.fft.rfftfreq(len(timearr), d = dt)
# plot power vs frequency for piano
plt.figure(figsize=(10,6))
plt.subplots_adjust(hspace = .9, )
plt.subplot(2,1,1)
plt.plot(freq,Ppiano,label='Ppiano')
plt.legend()
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.title('Piano Frequency Spectrum')
plt.grid()
# plot power vs frequency for trumpet
plt.subplot(2,1,2)
plt.plot(freq,Ptrumpet,label='Ptrumpet')
plt.legend()
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.title('Trumpet Frequency Spectrum')
plt.grid()
We used what we had learned in the first part of the lab to help write our own code that would perform the Fourier transform and read and plot the results for the four different devices' input signals and fft spectrums.