playing with fourier transforms and discrete informationby bob jantzen[execute this worksheet to restore output]1. Encoding the telephone numberwith(plots):Enter your 10 digit telephone number as a MAPLE list:N:=10:b:=[6,1,0,5,1,9,7,3,3,5];
and create the following wave function signal on the unit interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkobWZlbmNlZEdGJDYmLUYjNiUtSSNtbkdGJDYkUSIwRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiLEYnRjQvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjE2JFEiMUYnRjRGNC8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJw== of the space variable LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw==:s:=(x,t)->sum(b['n']*sin('n'*Pi*x)*cos('n'*Pi*t),'n'=1..N);
with initial value function at LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2I1EhRictRiM2JS1GLDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEiPUYnL0Y4USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGQi8lKXN0cmV0Y2h5R0ZCLyUqc3ltbWV0cmljR0ZCLyUobGFyZ2VvcEdGQi8lLm1vdmFibGVsaW1pdHNHRkIvJSdhY2NlbnRHRkIvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZRLUkjbW5HRiQ2JFEiMEYnRj5GKw== :f:=unapply(s(x,0),x): f(x);
which looks like:plot(f(x),x=0..1,color=green); signal0:=%:
2. Graphical representation of construction of initial signal in spaceWe take the first 10 harmonics on the unit interval (unit amplitude sine waves that vanish at the endpoints, labeled by the number of half oscillations from 1 to 10):plot([seq(sin(n*Pi*x),n=1..N)],x=0..1);
and weight them by the ordered coefficients which are your 10 digit telephone number:plot([seq(b[n]*sin(n*Pi*x),n=1..N)],x=0..1);
and then combine them by addition to form a waveform for a signal:plot(s(x,0),x=0..1,color=green); signal0:=%:
3. Comparison of different oscillation frequencies: animation of individual modesWe create a movie showing the first mode (one long bump) and the last mode (10 short bumps) in motion: the last mode moves at a frequency 10 times as fast as the first. You can use the VCR controls to slow it down in the MAPLE worksheet.animate(b[1]*sin(1*Pi*x)*cos(1*Pi*t),
x=0..1,t=0..2,color=red,
frames=50):
movie1:=%:animate(b[10]*sin(10*Pi*x)*cos(10*Pi*t),
x=0..1,t=0..2,color=blue,
frames=50,numpoints=80):
movie2:=%:display(movie1,movie2);
4. Animation of whole signalNext we view how this signal composed of all 10 modes behaves in time:animate( # make a movie of
s(x,t), # space function to be plotted
x=0..1, # space interval to be plotted
t=0..2, # time interval (one period)
frames=30, # number of frames in movie
numpoints=100, # number of points in each frame plot
color=green); # color of curves
5. Extracting the telephone number exactlyThrough the miracle of Fourier analysis, one can easily recover the coefficients of the different modes from the whole signal by an exact integration process. Setting the value of the digit LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== between 1 and 10 in the following integral process recovers that digit from the telephone number. Try a couple digits by changing its value 2 in the MAPLE worksheet.n:=2:
Int('2*f(x)*sin(n*Pi*x)',x=0..1)
=2*int(f(x)*sin(n*Pi*x),x=0..1);
n:='n':
6. Extracting the telephone number numericallyFirst we set the value of LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEiTUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== and determine the width of the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2I1EhRictRiM2JS1JI21uR0YkNiRRIjJGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUTEmSW52aXNpYmxlVGltZXM7RidGNS8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGPi8lKXN0cmV0Y2h5R0Y+LyUqc3ltbWV0cmljR0Y+LyUobGFyZ2VvcEdGPi8lLm1vdmFibGVsaW1pdHNHRj4vJSdhY2NlbnRHRj4vJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZNLUYsNiVRIk1GJy8lJ2l0YWxpY0dRJXRydWVGJy9GNlEnaXRhbGljRidGKw== intervals into which we divide the unit interval LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkobWZlbmNlZEdGJDYmLUYjNiUtSSNtbkdGJDYkUSIwRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVEiLEYnRjQvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHUSV0cnVlRicvJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRjE2JFEiMUYnRjRGNC8lJW9wZW5HUSJbRicvJSZjbG9zZUdRIl1GJw==:M:=11:Dx:=1.0/(2*M); # width of each of 2 M sampling subintervals
Now we create a MAPLE sequence of the numerical values of the subinterval endpoints LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiWEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYjLUYvNiVRImlGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJw== (labeled by integer subscripts from 1 to LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiUtSSNtbkdGJDYkUSIyRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUkjbW9HRiQ2LVExJkludmlzaWJsZVRpbWVzO0YnRjcvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkAvJSlzdHJldGNoeUdGQC8lKnN5bW1ldHJpY0dGQC8lKGxhcmdlb3BHRkAvJS5tb3ZhYmxlbGltaXRzR0ZALyUnYWNjZW50R0ZALyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTy1GLDYlUSJNRicvJSdpdGFsaWNHUSV0cnVlRicvRjhRJ2l0YWxpY0YnLUY7Ni1RIitGJ0Y3Rj5GQUZDRkVGR0ZJRksvRk5RLDAuMjIyMjIyMmVtRicvRlFGaG4tRjQ2JFEiMUYnRjdGKw==, one extra for the leftmost endpoint) and check the first and last such endpoints as a check:X:=[seq(0+1.0*i*Dx,i=0..2*M)]:
X[1],X[2*M+1];
Next we numerically sample (floating point evaluate) the inital signal function at these input values (we only show 4 digits to shorten the output)Y:=map(evalf,map(f,X)): evalf(%,4);
and zip the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== values together to form the list of coordinate pairs for these sampled points so that we can plot the data points together with the original initial signal function in the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw==-LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== plane:data:=zip((x,y)->[x,y],X,Y);
plot(data,style=point): dataplt:=%:
display(signal0,dataplt);
Now we need the numerical integration procedure to extract the coefficients of the modes from the numerical signal using the following Simpson's formula:'Dx'*Sum('( F(X[2*'i'-2+1])+4*F(X[2*'i'-1+1])+F(X[2*'i'+1]) )/3'
,'i'=1..M);
We create a 10 digit array to contain the numerical values of the coefficients whose rounded values should produce the 10 digits of our telephone number if we are successful.B:=array(1..N): convert(B,list): B[1];
Now we do a loop to produce each of the 10 mode coefficients.i:='i':for i from 1 to 10 do
F:=unapply(s(x,0)*2*sin(i*Pi*x),x):
B[i] :=
Dx*sum( evalf(
(F(X[2*'i'-2+1])+4*F(X[2*'i'-1+1])+F(X[2*'i'+1]))/3
),'i'=1..M);
od:
evalm(B);
We don't need all those decimals so we shorten up the numbers to 4 digits:map(z->evalf(z,4),B);
and now round them off to the nearest integer:map(round,B);
7. Finishing up: commentary and HTML export