<jupyter_start><jupyter_text>Assignment 7 Problem 1One thing we didn't do in our lectures, is to plot the Fourier Transform with the appropriate x-axis. The appropriate x-axis should be in frequencies. So how should we get the frequencies? A simple guess is, just take the inverse of time, $\frac{1}{T}$, and let it be our x-axis. (Later on, when I looked at the Jupyter notebook "fft_test" in the Data Analysis folder, the first example on the sine function does mention this too! The frequency is $f=\frac{k}{N}$.)Then we ran into another issue. the dates are in years or months for our CO2 data. Take the example of the data from "co2_mm_mlo.txt", we have a data point each month. Take the number of month, and let $1$ be divided by it, certainly makes no sense. So how about we treat time as if it is counted in months, instead of seconds, and we let the start time be 0? How do we implement this?Here I will follow the example in "fft_test", use the data in there and modify the example to illustrate our point.<jupyter_code># This is to enable linting, making sure we write our python code in style.
%load_ext pycodestyle_magic
%pycodestyle_on
import matplotlib.pyplot as plt

from fft import fft as fft
import numpy as np
import math
# Make the plots a bit bigger to see
# NOTE: Must be done in a separate cell
plt.rcParams['figure.dpi'] = 100
# A python code called read_co2 is used to help us read the data from files
from read_co2 import read_co2

# read_co2.py extracts two columns, correponding to our x and y, x is number of month,
# y is the average CO2 level in that month. Inside the bracket is the name of our data file.
x, y = read_co2('co2_mm_mlo.txt')
y_valid = y >= 0.
y = y[y_valid]
x = x[y_valid]

# To use fast Fourier transform based on the Cooley–Tukey algorithm, the number of
# our data points is in powers of 2, hence we need to cut the number of data.
# We could also choose 128 or 512, they are also powers of 2.
y = y[0:256]
x = x[0:256]

# Implementing the Fourier transform on our data
Y = fft(y)
Y_abs = abs(Y)

# Plot our data first
f1 = plt.figure(1)
plt.plot(x, y)
plt.ylabel('CO${2}$,ppm')
plt.xlabel('Date')
plt.show()<jupyter_output><empty_output><jupyter_text>Obviously our comments are a bit too long. Anyway, let's carry on, this is the plot of CO2 level VS month, notice we haven't plotted the Fourier transformed function yet. Because we need to get our x-axis redefined. I will first generate an array of numbers from 1 to 256. They will correspond to our actual time, sort of like 1 to 256 seconds, except here the time is counted in months, rather than seconds.<jupyter_code>N = 256
T = np.arange(1, N+1)

# If we just have np.range(N), this would give us an array starting with 0.
# But we will divide 1 with this array of numbers, so 0 is no good.
# We must start with 1
T<jupyter_output><empty_output><jupyter_text>Next, I think our frequency is just $\frac{1}{T}$<jupyter_code>w = 1/T<jupyter_output><empty_output><jupyter_text>Now we plot the fourier transform VS our frequency w<jupyter_code>f3 = plt.figure(3)
plt.plot(w, Y_abs)

plt.ylabel('Fourier component')
plt.yscale('log')
plt.xlabel('Frequency (month^(-1))')

plt.show()<jupyter_output><empty_output><jupyter_text>Now this looks more like a Fourier transform plot I think. We should zoom into the peaks at the start, I think that's where the genuine frequencies are. So we limit the range of our x-axis. I find that the most interesting peaks are between 0.0 to 0.1. The first really large peak is 1/256. This peak arises because we can think of the whole data set as one oscillation, that is not an interesting feature. We should probably ignore that too.<jupyter_code>f4 = plt.figure(4)
plt.plot(w, Y_abs)

plt.ylabel('Fourier component')
plt.yscale('log')
plt.xlabel('Frequency (month^(-1))')
plt.xlim([0.02, 0.10])

plt.show()<jupyter_output><empty_output><jupyter_text>The peak around the middle in the last graph, I think it actually corresponds to something! It is the yearly cycle I think, $12/256 = 0.046875$. Note that we have a logarithmic scale for the Fourier component. So the first peak is not that strong compared to the peak in the middle. So I would say that the most prominent frequency we have here is $0.046875$ per month, which correponds to the yearly oscillation. And we can see that in the first plot where we plotted the CO2 level VS real time, each bump in the curve correponds to a year.However, there is a little uneasiness in me. If we are talking about the yearly oscillation, shouldn't the freqeuncy be 1/12 per month? But that is only true if we count one year's data I think?On a second thought, if we take the inverse of $12/256$ and divide it by $256$, we get the frequency $1/12$. So maybe the frequency is actually given by $T/N$, where $N = 256$ is the number of data points. But the unit for $T/N$ is month, which is not the appropriate unit for frequency.---  Fixing our graphA problem with our 1/T approach is that, as we can see that our data points for w are not equally spaced!<jupyter_code>w<jupyter_output><empty_output><jupyter_text>There is also the problem of having to center our x-axis around 0 frequency. These two intricacies might have made our graph weird. But I think it still has the right spirit in it... How do we fix our graph then?Simply create an equally spaced array running from 0 to 0.5 (minimum and maximum freqeuncies in Fourier transform) in steps of $1/256$.<jupyter_code>f = np.arange(0, .5, 1/256)
len(f)<jupyter_output><empty_output><jupyter_text>Next we plot our graph versus f. Note there are only 128 data points in this array, (we ignore the negative frequencies) so we have to change the code slightly.<jupyter_code>f2 = plt.figure(2)
plt.plot(f, Y_abs[:128])

plt.ylabel('Fourier component')
plt.yscale('log')
plt.xlabel('Frequency (month^(-1))')


plt.show()<jupyter_output><empty_output><jupyter_text>Let's limit our x range slightly<jupyter_code>f2 = plt.figure(2)
plt.plot(f, Y_abs[:128])

plt.ylabel('Fourier component')
plt.yscale('log')
plt.xlabel('Frequency (month^(-1))')
plt.xlim([0.0, 0.3])

plt.show()<jupyter_output><empty_output><jupyter_text>Now we see a peak around $1/12=0.083$, yes! There is also a peak around $2/12=0.167$, this is called the second harmonic. Essentially, if we have one oscillation every year, there is also a repeat of the pattern every two years. This looks complicated, so the easiest way is to just use the numpy function fftfreq. It automatically generates the frequencies for us.<jupyter_code># input is (number of data points, step size)
freq = np.fft.fftfreq(256, 1)
f2 = plt.figure(2)
plt.plot(freq, Y_abs)

plt.ylabel('Fourier component')
plt.yscale('log')
plt.xlabel('Frequency (month^(-1))')

plt.show()<jupyter_output><empty_output><jupyter_text>---Next, we will do some filtering. I think I will first try filter out the higher frequencies. Here I borrowed the code from "fft_filtering.ipynb".<jupyter_code>from fft import ifft

# Get the FFT
Y = fft(y)
# Smooth the data in the Fourier domain.
# Adjust this to change the frequencies to delete.
maxfreq = 20
Y[maxfreq:len(Y)-maxfreq] = 0.0

# Now go back to the time domain.
# Compare the data before and after filtering.
yfiltered = ifft(Y)
yfiltered_abs = abs(yfiltered)

f6 = plt.figure(6)
plt.plot(x, y, label='Original')
plt.plot(x, yfiltered_abs, label='Filtered')
plt.ylabel('CO${2}$,ppm')
plt.xlabel('Date')
plt.show()<jupyter_output><empty_output><jupyter_text>This works. Can I modify it to filter out the lower frequencies?<jupyter_code># Get the FFT
Y = fft(y)
# Smooth the data in the Fourier domain.
# Adjust this to change the frequencies to delete.
minfreq = 20
Y[np.abs(Y) < minfreq] = 0.0

# Now go back to the time domain.
# Compare the data before and after filtering.
yfiltered = ifft(Y)
yfiltered_abs = abs(yfiltered)

f7 = plt.figure(7)
plt.plot(x, y, label='Original')
plt.plot(x, yfiltered_abs, label='Filtered')
plt.ylabel('CO${2}$,ppm')
plt.xlabel('Date')
plt.show()<jupyter_output><empty_output><jupyter_text>This also seems to work. Let's zoom in a little.<jupyter_code># Get the FFT
Y = fft(y)
# Smooth the data in the Fourier domain.
# Adjust this to change the frequencies to delete.
minfreq = 20
Y[np.abs(Y) < minfreq] = 0.0

# Now go back to the time domain.
# Compare the data before and after filtering.
yfiltered = ifft(Y)
yfiltered_abs = abs(yfiltered)

f8 = plt.figure(8)
plt.plot(x, y, label='Original')
plt.plot(x, yfiltered_abs, label='Filtered')
plt.ylabel('CO${2}$,ppm')
plt.xlabel('Date')
plt.xlim([1960, 1965])
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09578677839851024, 'fraction_numerical': 0.023044692737430168, 'mean_word_length': 4.751673360107095, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '48560', 'license': 'Apache-2.0', 'n_tokens_mistral': 2981, 'n_tokens_neox': 2749, 'n_words': 1213, 'path': 'HW7.1_fixed.ipynb', 'repo_name': 'piaohanx/CP1_GreenhouseGases', 'size': '8592'}