ABSTRACT (original article): We demonstrate the ability of a simple algorithm based on the venerable method of images (MOIs), to accurately model the detailed frequency response of a multidimensional, rectangular, lossy resonant cavity. The convergence properties of the model’s infinite series solution are shown to be determined by the cavity’s quality factor Q. A 1D example demonstrates that the MOI series converges to the exact solution. Next, a comparison to precisely measure 2D cavity data confirms that a straightforward extension of the 1D algorithm to multiple dimensions provides accurate results. The algorithm is short, easily understandable by undergraduate students and relatively undemanding to code. An example using ®Mathematica is provided. CITATION (original article): Frank Rice, Teresa Riedel, Isaiah Curtis; Modeling the frequency response of an acoustic cavity using the method of images. Am. J. Phys. 1 April 2023; 91 (4): 288–297. https://doi.org/10.1119/5.0124415
This is a heavily annotated version of a Method of Images (MOI) algorithm used to calculate the driven, steady-state response of a rectangular, lossy resonant cavity. The model is described in: “Modeling the frequency response of an acoustic cavity using the method of images”, Am. J. Phys. 91, 288–297 (2023), https://doi.org/10.1119/5.0124415. A manuscript of the article is available: http://sophphx.caltech.edu/MOI/MOImanuscript.pdf.
This Mathematica code is densely commented and not very compactly written to help make it more understandable. The default values used for various parameters match those used to generate the results presented in the article.
Perform example MOI function calculations:
Perform example MOI function calculations:
(set to False to skip those potentially time-consuming examples)
In[]:=
ExamplesQ=True;
Definitions of user interface symbols and their default values
Default values are what were used for paper’s model
Definitions of user interface symbols and their default values
Default values are what were used for paper’s model
Default values are what were used for paper’s model
Rectangular cavity characteristics (2 - D)
Rectangular cavity characteristics (2 - D)
Lx - Long dimension of cavity (will be in x direction) (meters)
Ly - Short dimension (y direction, meters)
Ly - Short dimension (y direction, meters)
In[]:=
Lx=0.152247(*meters*);Ly=0.114417(*meters*);
Speed of Sound for freq f ↔ wave number k conversions (m/sec):
Speed of Sound for freq f ↔ wave number k conversions (m/sec):
In[]:=
Cs=344.794(*meter/sec*);
Cavity quality factor Q as a function of frequency f (kHz):
Cavity quality factor Q as a function of frequency f (kHz):
This is the simple, linear Q(f) model used in the paper.
Change to Qfunc[f_]:=Q to use a constant value Q independent of frequency.
Change to Qfunc[f_]:=Q to use a constant value Q independent of frequency.
In[]:=
Qfunc[f_(*Hertz*)]:=EvaluateSimplify(f-1100.)+80.
220.-80.(*ΔQ*)
4500.-1100.(*Δf*)
Source and Receiver position offsets away from their respective, diagonally opposite corners
Source and Receiver position offsets away from their respective, diagonally opposite corners
All offsets should be nonnegative:
Sx - Source offset away from its corner in x direction (meters)
Sy - in y direction (meters)
Sy - in y direction (meters)
In[]:=
Sx=0.(*meters*);Sy=0.(*meters*);
x - Detector offset away from its corner in x direction (meters)
y - in y direction (meters)
Note: [x, y] values can be supplied as function arguments to Rtable (below). These definitions provide default values.
y - in y direction (meters)
Note: [x, y] values can be supplied as function arguments to Rtable (below). These definitions provide default values.
In[]:=
x=0.(*meters*);y=0.(*meters*);
Number of source image ‘clusters’ to use in the MOI calculation (see above graphic)
Number of source image ‘clusters’ to use in the MOI calculation (see above graphic)
Nmax - number of source image ‘clusters’ away from the source’s corner along the long cavity dimension Lx.
The image cluster array will extend ± 2NmaxLx away from source corner in both x and y directions, forming a
square-shaped area containing the source and its images.
The image cluster array will extend ± 2NmaxLx away from source corner in both x and y directions, forming a
square-shaped area containing the source and its images.
In[]:=
Nmax=300;
Functions for the MOI Algorithm
Functions for the MOI Algorithm
Rtable[]
Generate a list of distances from the array of image source locations to detector’s(x,y)
Rtable[]
Generate a list of distances from the array of image source locations to detector’s(x,y)
Generate a list of distances from the array of image source locations to detector’s(x,y)
Rtable[x,y] uses global values set for Lx, Ly, Sx, Sy, and Nmax.
Detector position offsets x, y are supplied as function arguments.
Returns a sorted list (largest first) of {count, distance} pairs to the source and its images.
‘count’ is the total number of images calculated to be at ‘distance’ from the detector.
Detector position offsets x, y are supplied as function arguments.
Returns a sorted list (largest first) of {count, distance} pairs to the source and its images.
‘count’ is the total number of images calculated to be at ‘distance’ from the detector.
Clear[Rtable];Rtable[x_,y_(*detectoroffsetpositionfromitscorner*)]:=Block{Rfunc,X=1.x,Y=1.y,Ny=Ceiling[(Lx/Ly)Nmax]},(*******************************************************)(*Rfunc[i,j]isalocalfunctionwhichreturnsthedistancesfrom{x,y}tothesourceimagesclusteredaroundanindexed[i,j]latticepoint.ByusingBlock[],itisgivenalocalredefinitionwitheachcalltoRtable.*)(*lookatsourceoffsetpositioninthecavity*)IfSx0.&&Sy0.,(*IfSx=Sy=0,thenthesourceisinitscorner,andonlyoneimage(ratherthan4)isfoundateachlatticepoint*)Rfunc=Compile{{i,_Integer},{j,_Integer}},(*compileforspeed*)Evaluate√(((2i-1)Lx+X)^2+((2j-1)Ly+Y)^2),RuntimeOptions"Speed",(*elsewecalculate4imagesaroundeachlatticepoint,±Sx,±Sy*)Rfunc=Compile{{i,_Integer},{j,_Integer}},Evaluate√(((2i-1)Lx-Sx+X)^2+((2j-1)Ly-Sy+Y)^2),√(((2i-1)Lx-Sx+X)^2+((2j-1)Ly+Sy+Y)^2),√(((2i-1)Lx+Sx+X)^2+((2j-1)Ly-Sy+Y)^2),√(((2i-1)Lx+Sx+X)^2+((2j-1)Ly+Sy+Y)^2),RuntimeOptions"Speed";(*endthelocaldefinitionofRfunc[]*)(*******************************************************)(*Here'swhatthenext,quitepowerful,expressiondoes:1.Table[]:iteratesoverthelatticeofimages,calculatingthedistancestotheimagesusingRfunc[i,j];2.//Flatten:getsridofallinnerbraces,leavingasimplelistofdistances;3.//Gather[]&:organizesthelistintosublistsofequal-valueddistances;4.//Sort[]&:sortsthesublistsfromlongesttoshortestdistancestoimproveprecision;5.//Map[]:replaceseachsublistwitha{count,distance}pair(bothwillbeoftypeReal).*)Table[Rfunc[i,j],{i,-Nmax,Nmax},{j,-Ny,Ny}]//Flatten//Gather[#,Equal]&//Sort[#,First[#1]>First[#2]&]&//Map[{1.Length[#],First[#]}&](*notethatcountwillbetypeReal*)(*endofBlock[]andthedefinitionofRtable*)
Rtable[] : use the default detector position offsets x, y assigned above.
In[]:=
Rtable[]:=Rtable[x,y]
Term[]: calculate one term of the sum of complex amplitudes from the sources
Term[]: calculate one term of the sum of complex amplitudes from the sources
The wave number argument k must have units of cycles/meter, not radians/meter. This makes it
quicker to get the phase of k×d, which will be 2π times its fractional part, rather that use modulo 2π.
Thus: k=f/Cs. (frequency f in Hertz)
quicker to get the phase of k×d, which will be 2π times its fractional part, rather that use modulo 2π.
Thus: k=f/Cs. (frequency f in Hertz)
In[]:=
Term=Compile {k,_Real},(*wavenumberincycles/meter(f/Cs)*) {Q,_Real},(*Qatthisk(orfreq)*) {m,_Real},(*countfromRtableentry*) {d,_Real}(*distancefromRtableentry*) , EvaluateExp[I2Pi*FractionalPart[k*d]]Exp[-Pi*k*d/Q]
m
d
Out[]=
CompiledFunction
Response[]: Total Term[] calls threaded over the list generated by Rtable[]
Response[]: Total Term[] calls threaded over the list generated by Rtable[]
Response takes a list of image source distances returned by Rtable as its first argument. If a single frequency (in Hertz) is its
second argument, Response returns the sum of Term[] calls over the elements in Rtable, generating
the complex-valued response phasor at the supplied frequency.
second argument, Response returns the sum of Term[] calls over the elements in Rtable, generating
the complex-valued response phasor at the supplied frequency.
In[]:=
Response RtableList:{{_Real,_Real}..},(*listofresultsfromcalltoRtable*) f_(*frequency(Hertz)tobeevaluated*) := Block {k=f/Cs,q=Qfunc[1.f]}, Total[Term[k,q,First[#],Last[#]]&/@RtableList]
If a list of frequencies (in Hertz) is input as Response’s second argument, then it returns a list of
{frequency, complex response} pairs suitable for generating frequency vs. gain or phase plots, etc. (examples below)
{frequency, complex response} pairs suitable for generating frequency vs. gain or phase plots, etc. (examples below)
Examples using the MOI Functions
Examples using the MOI Functions
Examples of Rtable[]
using various parameter values
Examples of Rtable[]
using various parameter values
using various parameter values
◼
Call with source and detector locations in opposite corners. The suffix //Timing
will return a list containing the CPU execution time (in seconds) required along with the result of the expression.
will return a list containing the CPU execution time (in seconds) required along with the result of the expression.
◼
Now call with the detector located in the center of the cavity. The list will be longer because
of the reduced symmetry of the images around the detector, taking a few more seconds to generate.
of the reduced symmetry of the images around the detector, taking a few more seconds to generate.
◼
With both the source and detector positioned away from their respective corners the symmetry is broken, and
the list of distances is approximately 16 times that for both in their corners. So is the time it takes to
generate the results. Using Block[] gives temporary, local values to Sx and Sy.
the list of distances is approximately 16 times that for both in their corners. So is the time it takes to
generate the results. Using Block[] gives temporary, local values to Sx and Sy.
Response[] Function Examples
Uses the list distances calculated by the Rtable[0,0] example above.
Response[] Function Examples
Uses the list distances calculated by the Rtable[0,0] example above.
Uses the list distances calculated by the Rtable[0,0] example above.
◼
First, we’ll use the response magnitude at 3kHz as the magnitude (response gain) reference.
◼
Now calculate the MOI model response at a few hundred frequencies and return the CPU execution time.
Scale the modeled response by the GainRef value calculated above.
Scale the modeled response by the GainRef value calculated above.
◼
Plot the resulting MOI model response gain and phase vs. frequency.
Going Deeper
Unwrap the phase plot and compare with the gain plot.
Going Deeper
Unwrap the phase plot and compare with the gain plot.
Unwrap the phase plot and compare with the gain plot.
◼
Unwrap[]: a function to remove phase data wrapping at ±180 degrees.
Use the function to unwrap the model phase data vs. frequency saved in modelphase.
Save the unwrapped phase data in unwrappedphase.
Use the function to unwrap the model phase data vs. frequency saved in modelphase.
Save the unwrapped phase data in unwrappedphase.
◼
Every response resonance is indicated by a phase change of 180 degrees. Thus, we next scale the unwrapped phase data to indicate an accumulated count of resonances passed as frequency goes up through the data range.
Note that every resonance peak is associated with an additional phase shift of 180 degrees: the response at frequencies above the resonance are phase shifted by 180 degrees relative to the response below the resonance. Conversely, an “anti-resonance,” a very sharp dip sometimes found between a pair of resonances, corresponds to a 180 degree phase shift with the opposite sign.
Thus, over the 0.9 - 5.0 kHz range modeled, the total phase shift is 8 × 180°: (14 resonances - 6 anti-resonances)×180°. We’ve counted the small peak and dip that nearly coincide at around 4.52 kHz as a resonance - anti-resonance pair, although they are too close in frequency to exhibit their full ±180° shifts.
Thus, over the 0.9 - 5.0 kHz range modeled, the total phase shift is 8 × 180°: (14 resonances - 6 anti-resonances)×180°. We’ve counted the small peak and dip that nearly coincide at around 4.52 kHz as a resonance - anti-resonance pair, although they are too close in frequency to exhibit their full ±180° shifts.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Modeling the frequency response of an acoustic cavity using the method of images
by Frank Rice, Teresa Riedel, Isaiah Curtis
Wolfram Community, STAFF PICKS, August 26, 2024
https://community.wolfram.com/groups/-/m/t/3256594
by Frank Rice, Teresa Riedel, Isaiah Curtis
Wolfram Community, STAFF PICKS, August 26, 2024
https://community.wolfram.com/groups/-/m/t/3256594