<jupyter_start><jupyter_text>ThinkDSPThis notebook contains code examples from Chapter 9: Differentiation and IntegrationCopyright 2015 <NAME>License: [Creative Commons Attribution 4.0 International](http://creativecommons.org/licenses/by/4.0/)<jupyter_code>from __future__ import print_function, division

import thinkdsp
import thinkplot
import math
import numpy as np
import pandas
import scipy.signal

PI2 = 2 * math.pi

%matplotlib inline<jupyter_output><empty_output><jupyter_text>As the first example, I'll look at daily closing stock prices for Facebook, from its IPO on 17 May 2012 to 8 December 2015 (note: the dataset includes only trading days )<jupyter_code>df = pandas.read_csv('fb.csv', parse_dates=[0])
ys = df.Close.values[::-1]
len(ys)<jupyter_output><empty_output><jupyter_text>Here's what the time series looks like:<jupyter_code>close = thinkdsp.Wave(ys, framerate=1)
close.plot()<jupyter_output><empty_output><jupyter_text>And here's the spectrum on a log-log scale.<jupyter_code>close_spectrum = close.make_spectrum()
close_spectrum.plot()
thinkplot.config(xscale='log', yscale='log', legend=False)<jupyter_output><empty_output><jupyter_text>The slope is -1.86, which is similar to red noise (which should have a slope of -2).<jupyter_code>close_spectrum.estimate_slope()<jupyter_output><empty_output><jupyter_text>We can use `np.diff` to compute the difference between successive elements, which is the daily change.<jupyter_code>change = thinkdsp.Wave(np.diff(ys), framerate=1)
change.plot()<jupyter_output><empty_output><jupyter_text>And here's the spectrum of the daily changes:<jupyter_code>change_spectrum = change.make_spectrum()
change_spectrum.plot()<jupyter_output><empty_output><jupyter_text>And here's the spectrum on a log-log scale.<jupyter_code>change_spectrum.plot()
thinkplot.config(xscale='log', yscale='log', legend=False)<jupyter_output><empty_output><jupyter_text>The estimated slope is close to zero, which is consistent with white noise.<jupyter_code>change_spectrum.estimate_slope()<jupyter_output><empty_output><jupyter_text>We can think the diff operation as convolution with a difference window, [1, -1].  Convolution with this window corresponds to multiplication by a filter.The following function computes the filter that corresponds to the window.<jupyter_code>def make_filter(window, wave):
    """Computes the filter that corresponds to a window.
    
    window: NumPy array
    wave: wave used to choose the length and framerate
    
    returns: new Spectrum
    """
    padded = thinkdsp.zero_pad(window, len(wave))
    wave = thinkdsp.Wave(padded, framerate=wave.framerate)
    philter = wave.make_spectrum()
    return philter<jupyter_output><empty_output><jupyter_text>And here's what the corresponding filter looks like:<jupyter_code>diff_window = np.array([1.0, -1.0])
diff_filter = make_filter(diff_window, close)
diff_filter.plot()<jupyter_output><empty_output><jupyter_text>And the angles:<jupyter_code>thinkplot.plot(diff_filter.angles())<jupyter_output><empty_output><jupyter_text>So we could also compute the daily changes by multiplying the spectrum of closing prices by the diff filter:<jupyter_code>change_spectrum2 = close_spectrum * diff_filter
change_spectrum2.plot()<jupyter_output><empty_output><jupyter_text>And then converting the spectrum to a wave.In the result, the first element is wrong, so I remove it.<jupyter_code>change2 = change_spectrum2.make_wave()
change2.ys = change2.ys[1:]
change2.plot()<jupyter_output><empty_output><jupyter_text>The we can confirm that we get the same result both ways (within floating point error)<jupyter_code>max(abs(change.ys - change2.ys))<jupyter_output><empty_output><jupyter_text>This diff operation is an approximation of differentiation, and we can compute the filter for differentiation analytically: each complex component is multiplied by $2 \pi i f$.<jupyter_code>deriv_filter = close.make_spectrum()
i = complex(0, 1)
deriv_filter.hs = PI2 * i * deriv_filter.fs
deriv_filter.plot()<jupyter_output><empty_output><jupyter_text>Now we can multiply the closing prices by the derivative filter:<jupyter_code>deriv_spectrum = close_spectrum * deriv_filter
deriv_spectrum.plot()<jupyter_output><empty_output><jupyter_text>The results are similar, with some differences due to (1) the difference window is only a coarse approximation of the derivative, especially at higher frequencies, and (2) the spectral derivative is based on the assumption that the signal is periodic, so the behavior at the beginning and end is different.<jupyter_code>deriv = deriv_spectrum.make_wave()
len(deriv), len(change)
deriv = deriv_spectrum.make_wave()
change.plot(color='gray', alpha=0.5)
deriv.plot(color='red', alpha=0.3)<jupyter_output><empty_output><jupyter_text>We can see the differences more clearly by zooming in on a section in the middle:<jupyter_code>low, high = 400, 430
thinkplot.plot(change.ys[low:high], color='gray', alpha=0.5)
thinkplot.plot(deriv.ys[low:high], color='red', alpha=0.3)<jupyter_output><empty_output><jupyter_text>The differences and derivative are similar in many places, but sometimes substantially different. Here's the difference between the derivative filter and the difference filter:<jupyter_code>deriv_filter.plot()
diff_filter.plot()<jupyter_output><empty_output><jupyter_text>The difference filter does not amplify the highest frequencies as much, which is why it is smoother than the derivative.  IntegrationNow let's think about integration.  We can compute the filter for integration analytically: each frequency component gets divided by $2 \pi i f$.I plot the result on a log-y scale so we can see it more clearly.<jupyter_code>integ_filter = deriv_filter.copy()
i = complex(0, 1)
integ_filter.hs = 1 / (PI2 * i * integ_filter.fs)
integ_filter.hs[0] = 0
integ_filter.plot()
thinkplot.config(yscale='log')
len(deriv_spectrum), len(integ_filter)<jupyter_output><empty_output><jupyter_text>We can confirm that the integration filter is correct by applying it to the spectrum of the derivative we just computed:<jupyter_code>integ_spectrum = deriv_spectrum * integ_filter
integ_spectrum.plot()
thinkplot.config(yscale='log', legend=False)<jupyter_output><empty_output><jupyter_text>And then converting back to a wave.  The result is identical to the daily closing prices we started with, but shifted so the mean is 0.  The problem is that the derivative clobbers the first element of the spectrum, which is the bias.  Once the bias information is lost, integration can't restore it.  So the result has an unspecified constant of integration.<jupyter_code>thinkplot.preplot(2)
integ_wave = integ_spectrum.make_wave()
close.plot(label='closing prices')
integ_wave.plot(label='integrated derivative')
shift = np.mean(close.ys) - np.mean(integ_wave.ys)
max(abs(integ_wave.ys-close.ys+shift))<jupyter_output><empty_output><jupyter_text>In the same way that the diff operator approximates differentiation, the cumulative sum approximates integration.  The following figure shows the cumulative sum of the daily changes:<jupyter_code>cumsum_wave = change.copy()
cumsum_wave.ys = np.cumsum(cumsum_wave.ys)
cumsum_wave.plot()<jupyter_output><empty_output><jupyter_text>Unlike the difference operator, it's not obvious how to write the cumulative sum in terms of convolution with a window.  But since cumsum is the inverse operation of diff, it's filter must be the inverse of `diff_filter`:<jupyter_code>cumsum_filter = make_filter(diff_window, change).invert()
cumsum_filter.plot()
thinkplot.config(yscale='log', legend=False)<jupyter_output><empty_output><jupyter_text>We can confirm that by computing the ratio of the cumsum spectrum to the change spectrum:<jupyter_code>cumsum_spectrum = cumsum_wave.make_spectrum()
ratio_spectrum = cumsum_spectrum.ratio(change_spectrum)
ratio_spectrum.plot(style='.', color='red', markersize=5)
thinkplot.config(yscale='log')<jupyter_output><empty_output><jupyter_text>The computed ratios have the same general shape as the filter, but they are noisy and might be offet by a constant factor.  I'm not sure what's going on there.<jupyter_code>thinkplot.preplot(2)
ratio_spectrum.plot(style='.', color='red', markersize=5)
cumsum_filter.plot(label='cumsum')
integ_filter.plot(label='integration')
thinkplot.config(yscale='log', legend=True)<jupyter_output><empty_output><jupyter_text>An example that's actually periodicWe can do the same thing using a periodic signal like a sawtooth:<jupyter_code>in_wave = thinkdsp.SawtoothSignal(freq=50).make_wave(duration=0.1, framerate=44100)
in_wave.plot()<jupyter_output><empty_output><jupyter_text>Here's the spectrum before the cumulative sum:<jupyter_code>in_spectrum = in_wave.make_spectrum()
in_spectrum.plot()<jupyter_output><empty_output><jupyter_text>The output wave is the cumulative sum of the input<jupyter_code>out_wave = in_wave.cumsum()
out_wave.unbias()
out_wave.plot()<jupyter_output><empty_output><jupyter_text>And here's its spectrum<jupyter_code>out_spectrum = out_wave.make_spectrum()
out_spectrum.plot()<jupyter_output><empty_output><jupyter_text>Now we compute the ratio of the output to the input:<jupyter_code>sum(in_spectrum.amps < 1), len(in_spectrum)<jupyter_output><empty_output><jupyter_text>In between the harmonics, the input componenents are small, so I set those ratios to NaN.<jupyter_code>ratio_spectrum = out_spectrum.ratio(in_spectrum, thresh=1e-11, val=np.nan)
ratio_spectrum.plot(style='.')

thinkplot.config(yscale='log', legend=False)
sum(np.isnan(ratio_spectrum.hs)), len(ratio_spectrum)<jupyter_output><empty_output><jupyter_text>To get the cumsum filter, I compute the diff filter again and invert it.<jupyter_code>diff_window = np.array([1.0, -1.0])
diff_filter2 = make_filter(diff_window, in_wave)
cumsum_filter2 = diff_filter2.invert()
cumsum_filter2.hs[0] = 0

thinkplot.preplot(1)
cumsum_filter2.plot()
thinkplot.config(yscale='log', legend=False)<jupyter_output><empty_output><jupyter_text>Finally, we can compare the computed ratios to the filter.  They match, confirming that the cumsum filter is the inverse of the diff filter.<jupyter_code>cumsum_filter2.plot(label='filter', color='gray', alpha=0.5, linewidth=6)
ratio_spectrum.plot(label='ratio', style='.', color='red', alpha=0.7, markersize=1)
thinkplot.config(yscale='log', legend=True)<jupyter_output><empty_output><jupyter_text>Now we can compute the output wave using the convolution theorem, and compare the results:<jupyter_code>len(in_spectrum), len(cumsum_filter2)
thinkplot.preplot(2)

out_wave.plot(label='cumsum')

out_wave2 = (in_spectrum * cumsum_filter2).make_wave()
out_wave2.plot(label='filter')

thinkplot.config(legend=True, loc='lower right')<jupyter_output><empty_output><jupyter_text>They are the same, within floating point error.<jupyter_code>max(abs(out_wave.ys - out_wave2.ys))<jupyter_output><empty_output><jupyter_text>Since we're figured out the filter that corresponds to cumulative sum, we can compute the window that corresponds to the filter.<jupyter_code>cumsum_window2 = cumsum_filter2.make_wave()
cumsum_window2.ys += 0
cumsum_window2.plot()<jupyter_output><empty_output><jupyter_text>Convolving with this window should be the same as applying the cumsum function.<jupyter_code>len(cumsum_window2), len(in_wave)
thinkplot.preplot(2)

out_wave.plot(label='cumsum')

out_wave3 = in_wave.convolve(cumsum_window2)
out_wave3.plot(label='convolution')

thinkplot.config(legend=True, loc='lower right')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0998433147632312, 'fraction_numerical': 0.01088091922005571, 'mean_word_length': 7.824116743471582, 'pattern_counts': {'":': 0, '<': 159, '<?xml version=': 0, '>': 158, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '39', 'dir': 'jupyter-structured-clean-dedup', 'id': '293117', 'license': 'MIT', 'n_tokens_mistral': 3865, 'n_tokens_neox': 3649, 'n_words': 969, 'path': 'ThinkDSP-master/code/diff_int.ipynb', 'repo_name': 'beckermy/popcorn-analyzer', 'size': '11494'}