The wave divisor function consists of a pulse outline modulated with a high frequency component. The real solution of the wave divisor function is:
$$ \Re(\sigma_{0})=\sum_{\mathbb{X}=2}^{\infty}\cos^{N} \left( \frac{\pi}{\mathbb{X}}x \right) \cos \left( \frac{N\pi}{\mathbb{X}}x \right) $$$N$ is determined by the pulse width of $cos^{N}$ and calculated with ($L$ pulseheight at position $\Delta x$). For every $\mathbb{X}$ a $N$ is calculated, this way all waves in the summation have similar pulsewidths. N should be an positive even integer to obtain positive pulses only.
More information: Wave Divisor Function, Wiki Fourier Transform, Wolfram Alpha
%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt1
import matplotlib.pyplot as plt2
import ipywidgets as widgets
from operator import add
from IPython.display import Audio
from IPython.display import display
divisors = 125
samplerate = 3000
samples = int(40001 * samplerate / 3000)
x1 = np.linspace(0, divisors, samples)
fig, ax1 = plt1.subplots(1, figsize=(9, 4))
plt1.suptitle('$\sigma_{0}$ Wave Divisor Function')
def update_plot(dx, L, wave):
ax1.clear()
#Set zero list
y=[0]*samples
#Calc Re divisor solution for all selected divisor waves
for w in wave:
N=-2*(w**2)*np.log(L)/((np.pi**2)*(dx**2))
N=2*round(0.5*N,0)
yw = ((np.cos(x1*np.pi/w))**N)*(np.cos(np.pi*N*x1/w))
y=list(map(add, y, yw) )
#Determine scaling for y axis (x=0 is excluded)
countMax=max(y[int(samples*(2)/100):samples])
countMin=min(y[int(samples*(2)/100):samples])
units = '$\Delta x$ = {}, $L$ = {}'
#update graph
ax1.plot(x1, y, label=units.format(dx, L))
ax1.axis([0, divisors, countMin-5,countMax+5])
ax1.legend(loc=1)
ax1.set_xlabel('$x$')
ax1.set_ylabel('$\sigma_{0}$')
ax1.grid(which='major', color='#666666', linestyle='-')
plt1.show()
print("Download WAV for playback")
display(Audio(y, rate = samplerate))
print("Options:")
print("L is pulse height at dx")
print("Select number of divisors waves (use: <Ctrl> or <Alt> to select multiples.)")
print("note: initial click on t=0, summation of all all waves occur there.)")
dx = widgets.FloatSlider(min=0.15, max=0.99, value=0.25, step=0.01, description='$\Delta x$:')
L = widgets.FloatSlider(min=0.15, max=0.99, value=0.5, step=0.01, description='$L$:')
wave = widgets.SelectMultiple(options=list(range(2,101)), value=list(range(2,101)), description="$\mathbb{X}$:")
widgets.interactive(update_plot, dx=dx, L=L, wave=wave)
Options: L is pulse height at dx Select number of divisors waves (use: <Ctrl> or <Alt> to select multiples.)
interactive(children=(FloatSlider(value=0.25, description='$\\Delta x$:', max=0.99, min=0.15, step=0.01), Floa…
The wave divisor function consists of a pulse outline modulated with a high frequency component. The real solution of the wave divisor function is:
$$ \Re(\sigma_{0})=\sum_{\mathbb{X}=2}^{\infty}\cos^{N} \left( \frac{\pi}{\mathbb{X}}x \right) \cos \left( \frac{N\pi}{\mathbb{X}}x \right) $$$N$ is determined by the pulse width of $cos^{N}$ and calculated with ($L$ pulseheight at position $\Delta x$). For every $\mathbb{X}$ a $N$ is calculated, this way all waves in the summation have similar pulsewidths. N should be an positive even integer to obtain positive pulses only:
$$ N(\mathbb{X}) = \frac{\log(L)}{\log \left( \cos \left( \frac {\pi}{\mathbb{X} } \Delta x \right) \right)} \approx - \frac{2 \mathbb{X}^2 \log(L)}{\pi^2 \Delta x^2} + \frac{\log(L)}{3}+ \mathcal{O} \left( \frac{1}{\mathbb{X}^2} \right)$$The first term $cos^N$ can also be simplified, this is the pulse outline. The pulse outline forms a bell shaped distribution arround the origin for $\mathbb{X} \rightarrow \infty$:
$$ O(x)=\lim_{\mathbb{X} \rightarrow \infty}\cos^{N} \left( \frac{\pi}{\mathbb{X}}x \right)= e^{a x^{2}}$$$$ a=\frac{\log(L) \space}{\Delta x^{2}}=constant$$The high frequency component $HF(\mathbb{X})$ scales linear with $\mathbb{X}$ (see link for more information) for: $\mathbb{X} \rightarrow \infty$.
$$ HF(\mathbb{X})= \cos \left( \frac{N\pi}{\mathbb{X}} x \right) \approx \cos (b x)$$$$ b(\mathbb{X}) = \frac{N}{\mathbb{X}}\pi \approx - \frac{2 \space \log(L)}{\pi \space \Delta x^{2}} \mathbb{X} = constant \cdot \mathbb{X}$$So for $\mathbb{X} \rightarrow \infty$ the wave divisor function becomes:
$$ \Re(\sigma_{0})\rightarrow \sum_{\mathbb{X}=2}^{\infty}e^{a x^{2}} \cos (b x) $$The wave divisor at infinity can be Fourier transformed in the frequency domain. The following Fourier transform definitation was used:
$$ \hat{f}(\xi)=\int_{-\infty}^{\infty}f(x) \mspace{3mu} e^{-2 \pi ix \xi} \mspace{3mu} dx$$With help of Wolfram Alpha the Fourier transform is determined (see link below). The frequency spectra of an individual divisor wave will consist of a bell shape mirrored in the y-axis.
$$ \hat{\sigma}_{0}(\xi)= \frac{\sqrt{\pi}}{2 \sqrt{-a}} \left( e^{(b-2 \pi \xi)^{2} /4a} + e^{(b+2 \pi \xi)^{2} /4a} \right) $$Every number will have at least on divisor wave. Because of the linearity properties of the Fourier transform we can sum the spectra to obtain the complete spectra of a number. The simulation below shows the time domain wave and the frequency spectra. Also the wave has been transposed to an audible signal.
More information: Wave Divisor Function, Wiki Fourier Transform, Wolfram Alpha
#Create Plot grid
fig= plt2.figure(figsize=(9, 4), constrained_layout=True)
widths = [4.5,4.5]
heights = [4]
gs=fig.add_gridspec(1,2,width_ratios=widths, height_ratios=heights, wspace=0.05)
ax1=fig.add_subplot(gs[0,0])
ax2=fig.add_subplot(gs[0,1])
#fig, ax2= plt2.subplots(1,2, figsize=(9, 4))
def update_plot(dx2, L2, sx):
xf=np.linspace(sx-0.5,sx+0.5,5000,endpoint=True)
ax1.clear()
ax2.clear()
reD=[0]*5000
#imD=[0]*2000
#amplification
amp=[10]*5000
#Create list with waves X=2 to X=100
wave2=list(range(2,101))
#Calculate Solution Wave Devisor Function
for w2 in wave2:
N2=(np.log(L2))/(np.log(np.cos(np.pi*dx2/w2)))
N2=2*round(0.5*N2,0)
reDw = ((np.cos(xf*np.pi/w2))**N2)*(np.cos(np.pi*N2*xf/w2))
#imDw = (-(np.cos(xf*np.pi/w2))**N2)*(np.sin(np.pi*N2*xf/w2))
reD=list(map(add, reD, reDw))
#imD=list(map(add, imD, imDw))
#Determine maximum Divisor Count
countD=max(reD)
#Plot Divisor Function
units2 = '$\Delta x$={}, $L$={}, $x$={}'
ax1.plot(xf, reD,color='#1f77b4', label=units2.format(dx2, L2, sx))
ax1.plot([sx],[countD], color='orange', marker='o')
ax1.legend(loc=2)
ax1.set_title('Divisor Wave $\sigma_{0}(x)$')
ax1.set_xlabel('$x$')
ax1.set_ylabel('$\sigma_{0}$')
ax1.axis([(sx-0.5), (sx+0.5), None,(countD+countD/3)])
ax1.grid(which='major', color='#666666', linestyle='-')
#Calculate Fourier Transform set amplitude summation 0.
ampliS=[0]*10000
#Maximum Frequency Range
N2=-2*(sx**2)*np.log(L2)/((np.pi**2)*(dx2**2))
N2=2*round(0.5*N2,0)
fmax=N2/sx
frange=np.linspace(-fmax,fmax,10000)
#Fourier Transform Calculated. Create graph label.
lab='Divisors of ' + str(sx) +':'
for w2 in wave2:
#Determine coeficients: a, b calculate Fourier Transform.
N2=(np.log(L2))/(np.log(np.cos(np.pi*dx2/w2)))
N2=2*round(0.5*N2,0)
a=np.log(L2)/(dx2**2)
b=np.pi*N2/w2
#b=-w2*(2/np.pi)*np.log(L2)/(dx2**2)
#Only add waves from divisors of x (modules).
if (sx%w2)==0:
Spec = (np.sqrt(np.pi))/(2*np.sqrt(-a))*(np.exp(((b-2*np.pi*frange)**2)/(4*a)) + np.exp(((b+2*np.pi*frange)**2)/(4*a)))
ampliS=list(map(add, ampliS, Spec))
lab=lab+'\n $\mathbb{X}$='+str(w2) + '$, f$='+str(np.round(0.5*N2/w2,1))
#Plot individual divisor frequencies
ax2.fill_between(frange,Spec, color='orange')
#Plot summation frequencies.
ax2.set_title('Divisor Spectrum $\hat{\sigma}_{0}$')
ax2.annotate(lab, xy=(fmax-fmax/3.1,0.01))
ax2.plot(frange, ampliS,color='#1f77b4')
ax2.set_xlabel('$Frequency$')
ax2.set_ylabel('$\hat{\sigma}_{0}$')
ax2.axis([0,fmax, 0,None])
ax2.grid(which='major', color='#666666', linestyle='-')
plt2.show();
#print(np.sqrt(-a*2))
#Create Audiofile
display(Audio(reD, rate=20000))
#display(Audio(imD, rate=20000))
print('Options: L is pulse height at dx, Audio is from Real Part.')
print('Orange dot in Wave graph indicates divisor count of: x')
print('Blue line in spectrum indicates total spectrum')
dx2 = widgets.Dropdown(options=[0.05, 0.10,0.15, 0.20,0.25,0.30,0.35,0.40,0.45,0.5], value=0.2, description='$\Delta x$:')
L2 = widgets.Dropdown(options=[0.10,0.15, 0.20,0.25,0.30,0.35,0.40,0.45,0.5], value=0.5, description='$L$:')
sx = widgets.Dropdown(options=list(range(2,101)), description='$x$:',value=30)
widgets.interactive(update_plot, dx2=dx2, L2=L2, sx=sx)
Options: L is pulse height at dx, Audio is from Real Part. Orange dot in Wave graph indicates divisor count of: x Blue line in spectrum indicates total spectrum
interactive(children=(Dropdown(description='$\\Delta x$:', index=3, options=(0.05, 0.1, 0.15, 0.2, 0.25, 0.3, …
An movie has been created where the divisor function and audio are synchronized. The movie has been created in the range till x=1000. The displayed movie below is done with the following settings: dx=0.15, L=0.5, Rate=3000, BPM=600. For other settings see my youtube channel.
from IPython.display import YouTubeVideo
#Wave divisor function audio
YouTubeVideo('8qtZJ6yp5D0')
When the pulsewidth in the time domain gets smaller the frequency spectrum of the divisors tend to be identified more clear. From Fourier transform properties one would expect the frequency bandwidth to become wider as the time domain pulse gets narrow: Uncertainty principle.
The spectrum of the wave divisor function seems to behave opposite to the uncertainty principle. Below the z-score of the wave divisor spectra is calculated. The z-score describes the behaviour and the uncertainty principle is maintained.
Time domain $f(x)$:
$$ \Re(\sigma_{0})\rightarrow \sum_{\mathbb{X}=2}^{\infty}e^{a x^{2}} \cos (b x) $$The pulsewidth in the time domain is determined by: $L$ pulseheight at position $\Delta x$. In the equations described later we will vary the pulsewidth in the time domain. Onward we set $L=0.5$ as an constant and the time domain pulsewidth is varied by reducing $\Delta x \rightarrow 0$.
$$ a=\frac{\log(L) }{\Delta x^{2}}=constant$$$$ b(\mathbb{X}) = \frac{N}{\mathbb{X}}\pi \approx - \frac{2 \log(L)}{\pi \Delta x^{2}} \mathbb{X} = constant \cdot \mathbb{X}$$Frequency domain $\hat{f} (\xi)$:
$$ \hat{\sigma}_{0}(\xi)= \frac{\sqrt{\pi}}{2 \sqrt{-a}} \left( e^{(b-2 \pi \xi)^{2} /4a} + e^{(b+2 \pi \xi)^{2} /4a} \right) $$The frequency pulses can be seen as normal distributions . The standard deviation of a pulse in the frequency domain is proportional to:
$$ Stdev(\hat{\sigma}_{0}(\xi)) \propto \sqrt{-a}$$The minimal frequency distance between two neigbour pulses is:
$$ \Delta \xi = b(\mathbb{X}+1)-b(\mathbb{X})=b(1)$$The z-score between to neighbour frequency pulses then is:
$$ Z \propto \frac{b(1)}{\sqrt{-a}} \propto \frac{1}{\Delta x}$$When the time domain pulse gets narrow $\Delta x \rightarrow 0$ the $z-score$ in the frequency domain gets bigger. Thus the individual pulses in the frequency domain become better identified. One can say that the pulsewidth in frequency domain $\sqrt{-a}$ grows more slowly then the frequency difference between two neighbour divisors $b$.
More information: Uncertainty principle, z-score