New and Improved Data

Length[CircularDrumData]

CircularDrumData={310.,527,737,1018,1292,1495,1768,2006,2277,2552,2780,3040,3379,3632,3892,4158,4396,4725,5022,5295,5947};
ScaledCirc=(CircularDrumData/(Min[CircularDrumData]))

Clear[α]
ScaledTheoreticalFreq[α_]:=Module[
{PredictedFreq},
PredictedFreq= Flatten[
Table[
N[(Sqrt[((N[BesselJZero[n,k]])^2)-(α^2)])],{n,0, 5},{k,1,5}]
];
PredictedFreq=Sort[Select[PredictedFreq,Im[#]==0.0&]];
PredictedFreq/Min[PredictedFreq]
]
ScaledCirc
α=11.53;
ScaledTheoreticalFreq[α]

Manipulate[
ListPlot[ 
ScaledTheoreticalFreq[α],
GridLines->{Automatic,ScaledCirc},
PlotRange->All
],
{α,8,12,0.1}]

ScaledTheoreticalFreq[α]

Wave Eq: 
Dirichlet Boundary Conditions

So far we have solved eigen problems for the -Laplacian 
	Δϕ=-λ ϕ
on several domains.

Heat Eq

We solved the heat equation
	(Subscript[u, t])=k Δu
on these domains by seperating the time behavior
	u=T[t] ϕ[space]
and plugging in 
	T'[t]ϕ[space]=k T[t]Δϕ[space]
and separating to get 
	((T'[t])/(k T[t]))=((Δϕ[space])/(ϕ[space]))=-λ
and solving the T equation to get 
	T[t]=(e^(-λ k t))
where λ was one of the eigenvalues of -Δ.

Wave Eq

We solve the wave equation
	(Subscript[u, tt])=k Δu
on these domains in the same way. We seperate the time behavior
	u=T[t] ϕ[space]
and plugging in
	T''[t]ϕ[space]=k T[t]Δϕ[space]
and separating to get 
	((T''[t])/(k T[t]))=((Δϕ[space])/(ϕ[space]))=-λ
where λ is any of the eigenvalues of -Δ.   The time equation is now
	T''[t]+λ T[t]==0
which has solutions 
	T[t]=(Subscript[C, 1])Cos[(Sqrt[λ])t]+(Subscript[C, 2])Sin[(Sqrt[λ])t]
since all the eigenvalues λ are positive.
If we have the eigenfunctions etc it is not much harder to solve the wave equation than the heat equation.

Damped Wave Eq

The wave equation
	(Subscript[u, tt])+α (Subscript[u, t])=k Δu
works the same. Seperate the time behavior
	u=T[t] ϕ[space]
and plugging in
	T''[t]ϕ[space]+α T'[t] ϕ[space]=k T[t]Δϕ[space]
and separating to get 
	((T''[t] + α T'[t])/(k T[t]))=((Δϕ[space])/(ϕ[space]))=-λ
where λ is any of the eigenvalues of -Δ.   The time equation is now
	T''[t]+α T'[t]+λ T[t]==0
which is again something we know how to solve from our ODE course.

Our Geometry List now

1D Interval 0<X<L

Dirichlet: ϕ[0]=ϕ[L]=0

Mode Shape (Subscript[ϕ, n])=Sin[n π x]

Wave Equation Time Dependence
(Subscript[a, n])Cos[(Subscript[ω, n])t]+(Subscript[b, n])Sin[(Subscript[ω, n])t]

(Subscript[ω, n])=(Sqrt[(((nπ/L))^2)])=n (π/L)

2D Disk 0<r<(Subscript[r, max]) and 0<θ<2π

Dirichlet: ϕ[(Subscript[r, max]),θ]=0

Mode Shapes
(BesselJ[n,  (BesselJZero[n, k]) r/(Subscript[r, max])])Cos[n θ] and
(BesselJ[n,  (BesselJZero[n, k])   r/(Subscript[r, max])])Sin[n θ]

Wave Equation Time Dependence
(Subscript[a, n,k])Cos[(Subscript[ω, n,k])t]+(Subscript[b, n,k])Sin[(Subscript[ω, n,k])t]

(Subscript[ω, n])=(Sqrt[((BesselJZero[n,k])^2)/((Subscript[r, max])^2)])=BesselJZero[n,k]

2D Rectangle 0<x<a and 0<y<b

Dirichlet: ϕ[0,y]=ϕ[a,y]=ϕ[x,0]=ϕ[x,b]=0

Mode Shapes
Sin[n π x/a] Sin[m π y/b] and

Wave Equation Time Dependence
(Subscript[a, n,m])Cos[(Subscript[ω, n,m])t]+(Subscript[b, n,m])Sin[(Subscript[ω, n,m])t]

(Subscript[ω, n])=(Sqrt[(π^2)(((m^2)/(b^2))+ ((n^2)/(a^2))) ])=(π/a)(Sqrt[(n^2)+(((a^2)(m^2))/(b^2)) ])=(Sqrt[(n^2)+(ϵ^2)(m^2) ])
where ϵ =a/b.

Predicted Scaled Spectra(with damping)

If we include damping into the mix then we should have the PDE
	(Subscript[u, tt])+2α (Subscript[u, t])=Δu
and as noted before the eigen problem -Δϕ=λ ϕ is the same as before while the temporal behavior for the mode with eigenvalue λ is 
	((T''+2α T')/T)=-λ
or equivalently
	 T''+2α T'+λ T=0
which has solutions 
	T=(e^(-α t))((Subscript[C, 1])cos(β t)+(Subscript[C, 2])sin(β t))
where β=(Sqrt[λ-(α^2)]) provided λ>(α^2).  As you can tell the 2 is folded out of the damping constant purely for neatness.

In other words the damping causes the mode amplitude to decay exponentially as well as shifting the frequency slightly downwards.  The big question is how can we measure the damping in each of our drums.  Two ways come to mind.  We could fit the frequencies and the damping from the peaks we have in the data.  We could fit the damping from the decay of one of the frequencies.

Clear[PianoFreq]
PianoFreq[α_]:= Table[ ((Sqrt[(((n π))^2) - (α^2)])/(Sqrt[(π^2)-(α^2)])),{n,1,9}]
Plot[PianoFreq[α],{α,0,2}]

Clear[CircularDrumFreq]
CircularDrumFreq[α_]=Sort[
Flatten[
Table[
((Sqrt[((N[BesselJZero[n,k]])^2)-(α^2)])/(Sqrt[((N[BesselJZero[0,1]])^2)-(α^2)])),{n,0,3},{k,1,3}]
]
];
Plot[
Evaluate[CircularDrumFreq[α]],
{α,0,2}]

RectangularDrumFreq[α_,ϵ_]:=Sort[Flatten[
N[Table[
((Sqrt[((n^2)+(ϵ^2)(m^2))(π^2) -(α^2)])/(Sqrt[(1+(ϵ^2))(π^2)-(α^2)])),{n,1,3},{m,1,3}]
]
]
]
Plot[
Evaluate[RectangularDrumFreq[α,1.4]],
{α,0,4}]

Estimating α from the Mode Decay

Lets see if we can estimate the α for the Piano example sound.

snd=ExampleData[{"Sound","Piano"}]
Spectrogram[snd]

First we need the data. Here I am

Map[Head,snd[[1]]]
snd[[1,2]]
Dimensions[snd[[1,1]]]

Here are both tracks.

{Track1, Track2}=snd[[1,1]];
ListPlot[ {Track1, Track2},
PlotLegends->{"1","2"},
PlotRange->All
]

ListLogPlot[ Abs[Track1]]
ListLogPlot[ Abs[Track2]]

The first chunk of Track1 has the initial decay of the string vibrations in it.

CutTrack1=Track1[[3000;;14000]];
ListPlot[CutTrack1]
ListLogPlot[Abs[CutTrack1]]
Length[CutTrack1]

We need to look at the Frequency data.  The tool everyone woudl use is the Fast Fourier Trabsform (FFT).  The command in Mathematica is Fourier. The FFT data is comple.  For real input the Real part of the FFT is symmetric about the middle while the Imaginary part of the FFT is anti-symmetric about the middle. For a clean sound the large entries in the FFT are near the edges.

FreqData=Fourier[CutTrack1];
ListPlot[Im[FreqData]];
ListLogPlot[Abs[FreqData],
PlotRange->All]

We are going to focus in on the left hand side.  Something very cool is going on right here because you can see the frequency shift that damping should put into the harmonics.

Fund = 263;
ListLogPlot[Abs[FreqData],
PlotRange->{{1, 2000},All},
Frame->True,
GridLines->{ Fund* Range[1, 7],Automatic}
]

ListPlot[ Log[10,Abs[FreqData]],
PlotRange->{{1, 2000},All},
Frame->True,
GridLines->{ Fund* Range[1, 7],Automatic}
]

But first we need to understand what the FFT means. First location: the peak at 263 means that there is a lot of 263 wiggle in the sample which is of length 11,001;  the peak at 263*2 means there is quite a lot of 526 wiggle in the sample; etc.  Now height, the size of the peak says how much of the various wiggles are in the sample:  the peak at 263 is about (10^1) tall; the peak at 526 is about (10^(-0.4)) tall; the peak at 1578 is about (10^(-0.7)) tall.  This tells you that the weights for the various wiggles are roughly {10,(10^(-0.4)),(10^(-0.7))}={10,0.4,0.2}.  The last thing you need to know is that the real part is a cosine and the imaginary part is a sine - since you can not hear the difference it does not matter and we are going to use Cosine

Fund =1;
Plot[ 10 Cos[ Fund ((2 π t)/11001)]+ 0.4 Cos[ 2 Fund ((2 π t)/11001)]+0.2 Cos[Fund 3 ((2 π t)/11001)],{t,0,11001}]

We can resynthesize an approximate noise to reassure ourselves that we are not doing anything to crazy.

Fund =263;
sampleRate=22050;
Play[ 
 10 Cos[ Fund (sampleRate/11001) 2 π t]+ 0.4 Cos[ 2 Fund (sampleRate/11001) 2 π t]+0.2 Cos[Fund 3 (sampleRate/11001)2 π t],

{t, 0,0.5}]

It sounds lousy becaue we do not have the decay in there it also sounds lousy cause we do not have all the overtones in.  Here is the original and a spectrogram of our first snippet.

ListPlay[CutTrack1,SampleRate->22050]
Spectrogram[CutTrack1]

What we want to do is see how the coefficients decay with time.  What we have right now is the average.  We need to chop the data up into pieces and see how the intensity changes as time passes. We need to pick a chunk size to analyze: It should be big enough to have a fair number of wiggles but small enough that we get a fair number of samples out of our 0.5 seconds of data.

Fund=263;
SingleWiggle=Ceiling[(11001/Fund)];
WiggleCount=23;
SampleCutData=CutTrack1[[1;;WiggleCount*SingleWiggle]];
FreqSampleCutData=Fourier[SampleCutData];
ListPlot[SampleCutData]

CloseUpSize=50;
ListLogPlot[Abs[FreqSampleCutData][[1;;CloseUpSize]],
PlotRange->All,
GridLines->{{WiggleCount},Automatic},
Frame->True,
DataRange->{0,CloseUpSize-1}
]

5

Now all we need to do is sweep through the data and see how the coefficient for the fundamental mode changes as we go down.

WiggleCount=3;
FData=Table[
Abs[ Fourier[CutTrack1[[n;;n+WiggleCount*SingleWiggle]]][[1;;2WiggleCount]]],
{n, 1, Length[CutTrack1]-WiggleCount*SingleWiggle}
];
ListLogPlot[Transpose[FData],
PlotRange->All,
Frame->True,
PlotLabel->StringForm["WiggleCount = ``",WiggleCount],
PlotLegends->SwatchLegend[Automatic,LegendMarkerSize->8,LegendFunction->Frame]
]

PeakData=FData[[All,WiggleCount+1]];
LogPeakData=Log[PeakData];
ListPlot[LogPeakData]
logFit[x_]=Fit[LogPeakData,{1,x},x]

Lets check how the thing worked.

funFit[x_]=(E^(logFit[x]))
Show[
LogPlot[ funFit[x],{x,1, (10^4)}],
ListLogPlot[PeakData]
]