Lab 8 - ManuBarcellos/computational_physics GitHub Wiki

Goal:

The goal of this lab is to explore the Fourier transform and understand how complex waveforms can be broken down into simpler sine and cosine components. We’ll sum multiple cosine waves and examine their combination, analyze real-world data like sunspot counts using FFT (Fast Fourier Transform), reconstruct complex signals (like square waves) using sine components, and compare musical instrument sounds through spectral analysis.

Overview:

Part 1: Summing cosine functions

We begin by plotting three individual cosine functions and their sum. Each cosine has a different amplitude, frequency, and phase shift.

# use linspace to set up time array called t
t = np.linspace(0, 4, 100)
omega = [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
f1 = A[0] * np.cos(omega[0] * t + phi[0])
f2 = A[1] * np.cos(omega[1] * t + phi[1])
f3 = A[2] * np.cos(omega[2] * t + phi[2])

fsum = f1 + f2 + f3

Plot the individual and summed cosine waves:

# create your graph here
plt.subplots_adjust(hspace = 0.6, wspace = 0.4)
plt.subplot(2, 2, 1)
plt.plot(t, f1)
plt.title('f1 Graph')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

plt.subplot(2, 2, 2)
plt.plot(t, f2)
plt.title('f2 Graph')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

plt.subplot(2, 2, 3)
plt.plot(t, f3)
plt.title('f3 Graph')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

plt.subplot(2, 2, 4)
plt.plot(t, fsum)
plt.title('f1 + f2 + f3 Graph')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

Part 2: Computing the Fourier Transform of the Summed Wave

We analyzed the frequency content of the summed waveform created in Part 1 using the Fast Fourier Transform (FFT). This tells us which frequencies contribute most to the shape of the waveform.

# 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(fsum)

# set power equal to the absolute value of Y
power = np.abs(Y)/50 # 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
# label the axes
# add a grid using plt.grid()
plt.plot(freq, power)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power (W)')
plt.title('Power x Frequency')
plt.grid()

This graph shows which frequencies dominate in the summed waveform.

To confirm that the Fourier transform contains all the information needed to reconstruct the original wave, we apply the inverse FFT and compare it with the original sum.

# 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 = 'Inverse Fourier vs t')

# overplot sumwave vs t using a dashed line
plt.plot(t, fsum, '--', label = 'f1 + f2 + f3 vs t')
plt.title('Inverse Fourier Transform')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()

The two lines should perfectly overlap, demonstrating that the Fourier transform preserves all waveform information.

0.5 Hz, 1.0 Hz, and 2.0 Hz are the frequencies that have a significant amplitude
The significant frequencies in the Fourier transform exactly match the frequencies of the original three cosine waves.
The amplitude of the Fourier transform corresponds to the amplitude of the original cosine waves (the sum of the 3 cosine waves).

Part 3: Analyzing the Periodicity in Sunspots

In this part, we’ll explore the periodic behavior of sunspots using real-world data and Fourier analysis. We'll identify the sunspot cycle's period and compare it to known scientific values.
We begin by downloading and reading the file sunspots.txt, which contains monthly sunspot counts starting from January 1749.

# wget command here
!wget http://www-personal.umich.edu/~mejn/cp/data/sunspots.txt

# code to read in file
# save the month as month and the number of sunspots as nspots
month, nspots = np.loadtxt('sunspots.txt',unpack=True)

Plotting Number of Sunspots vs Months

# make your graph here
plt.figure(figsize=(20, 5))
plt.plot(month, nspots)
plt.xlim(0, 500)
plt.title('Number of Sunspots vs Months')
plt.xlabel('Months')
plt.ylabel('# of Sunspots')
plt.grid()

The estimate of the period, by looking at the graph, is about 130 months.

We now apply a Fourier transform to analyze the periodicity mathematically.


# run np.fft.rfft on nspots
X = np.fft.rfft(nspots)

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

# calculate the corresponding frequencies using np.fft.rfftfreq
n = len(month)
dt = month[1] - month[0]
frequency1 = np.fft.rfftfreq(n, d = dt)

# calculate the period corresponding to each frequency
period = 1/frequency1

# plot the power versus period, rather than versus frequency
plt.plot(period, power1)
plt.title('Power vs Period')
plt.xlabel('Power (W)')
plt.ylabel('Period (Months)')
#plt.xlim(0, 500)

print(f'The max period is: {period[np.argmax(power1[1:])]:.3f} months')

The max period is: 136.652 months.

The period for the Power vs Period graph is around 136 months, which is very close to the period that I estimated on the Number of Sunspots vs Time graph, which was 130 months.
The known period is around 11 years, 132 months. Therefore, the estimators were really close to the actual value.

Part 4: Building a square wave from cosine functions

-In this part, we approximate a square wave by summing multiple sine waves of increasing frequencies and decreasing amplitudes

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

t = np.linspace(0, 4, 1000)
sum = np.zeros(len(t))
# try adding 50 sine functions together
# let the frequency increase each time, and
# let the amplitude decrease as the frequency increases

for i in range(1,100,2):
  # compute a new sine function
  amp = a0 / i
  freq = w0 * i
  sine = amp*np.sin(freq*t)

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

  # add the new function to the sum
  sum += sine

# plot the sumwave vs time
plt.plot(t, sum, 'black')
plt.title('Square Wave')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

Part 5: Analyzing Sound of Musical Instruments

We analyzed how the same musical note differs when played by two different instruments, a piano and a trumpet, by examining their waveforms and frequency content using Fourier Transforms.
We used wget to download two .txt files containing amplitude data sampled, and plotted both of the waveforms on the same figure with one column and two rows.

# wget files
!wget http://www-personal.umich.edu/~mejn/cp/data/piano.txt
!wget http://www-personal.umich.edu/~mejn/cp/data/trumpet.txt


# load files
piano = np.loadtxt('piano.txt',unpack=True)
trumpet = np.loadtxt('trumpet.txt',unpack=True)

# plot the waveforms
plt.subplots_adjust(hspace = 0.6, wspace = 0.4)
plt.subplot(2, 1, 1)
plt.plot(piano)
plt.title('Piano')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

plt.subplot(2, 1, 2)
plt.plot(trumpet)
plt.title('Trumpet')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

The piano waveform appears smoother and more sinusoidal, while the trumpet waveform shows more complex oscillations, indicating a richer harmonic structure.
We applied the Real Fast Fourier Transform (rFFT) to convert the time-domain signals into the frequency domain, and plotted the Power vs Frequency Graph.

# 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 = np.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 = np.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 = (15,5))
plt.subplot(1, 2, 1)
plt.plot(freq, Ppiano)
plt.title('Piano')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power (W)')
#plt.xlim(400, 600)

# plot power vs frequency for trumpet
plt.subplot(1, 2, 2)
plt.plot(freq, Ptrumpet)
plt.title('Trumpet')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power (W)')
#plt.xlim(400, 600)

The piano does not have as prominent overtones as the trumpet. This can be seen on the graph because while the piano has one distinct peak at the frequency of the intended note, the trumpet has its biggest peak, surrounded by other sizeable peaks that indicate overtones of the desired note.
From the Fourier transform results, both instruments are playing at C5, as the peak for the piano and the trumpet were around 525 Hz.
That note is audible, as humans can hear from 20 Hz to 20,000 Hz.