1

No.86

Nov. 2021

ISSN 0917-2211

Variable Star Bulletin
A code for two-dimensional frequency analysis using
the Least Absolute Shrinkage and Selection Operator (Lasso)
for multidisciplinary use
Taichi Kato1
1

Department of Astronomy, Kyoto University, Sakyo-ku, Kyoto 606-8502, Japantkato@kusastro.kyoto-u.ac.jp

Received 2021 Nov. 19

Abstract
In Kato and Uemura (2012), we introduced the Least Absolute Shrinkage and Selection Operator (Lasso)
method, a kind of sparse modeling, to study frequency structures of variable stars. A very high frequency
resolution was achieved compared to traditional Fourier-type frequency analysis. This method has been
extended to two-dimensional frequency analysis to obtain dynamic spectra. This two-dimensional Lassofrequency analysis yielded a wide range of results including separation of the orbital, superhump and negative
superhump signals in Kepler data of SU UMa stars. In this paper, I brieﬂy reviewed the progress andapplications of this method. I present a full R code with examples of its usage. This code has been conﬁrmed
to detect the appearance of the orbital signal and the variation of the spin period after the eruption of the nova
V1674 Her. This code also can be used in multidisciplinary purposes and I provide applications to analysis ofavian vocalizations. I found ﬁne structures in the call of the Eurasian wren (Troglodytes troglodytes), which
is likely to be used species identiﬁcation. This code would be a new tool in studying avian vocalizations withhigh temporal and frequency resolutions. Interpretations of the power spectra of avian vocalizations will also
be helpful in interpreting the power spectra of variable stars.

1

Introduction

Frequency analysis is essential in many ﬁelds of science. In astronomy (I particularly deal with period/frequency
analysis of variable stars), many types methods have been used. Fourier-type methods include discrete Fouriertransform (Deeming 1975), least-squares ﬁt of sinusoids for unevenly spaced data (Lomb-Scargle Periodogram,Scargle 1982; Horne and Baliunas 1986; Zechmeister and Kürster 2009) and CLEAN algorithm (Roberts et al.
1987). They are related to a variety of Cohen’s class transformations for evenly sampled data, which include
short-time Fourier transform (STFT), wavelet transform and Wigner-Ville distribution (Hlawatsch and Auger
2008). These methods give generally the same frequency resolution given by the Heisenberg-Gabor limit.
There are also methods by evaluating dispersions in phase-sorted or phase-binned data. This class of approach
includes the string-length method (Dworetsky 1983), Laﬂer-Kinman’s method (Laﬂer and Kinman 1965), Phase
Dispersion Minimization (PDM, Stellingwerf 1978) and Analysis of Variance (AoV, Schwarzenberg-Czerny 1989).
These methods give an estimate with a smaller error than Cohen’s class transformations for the frequency/period
when a single signal is present (Fernie 1989; Kato et al. 2010). These methods are advantageous over Cohen’sclass transformations for a non-sinusoidal signal, which is usual for variable stars.
In recent years, a completely new approach in Fourier-type analysis has developed with an assumption that
only small number of frequencies are present at the same time (condition usually called “sparse”). This is an2application of a rapidly developing ﬁeld of compressed sensing (e.g. Donoho 2006). Under this assumption,
we can, in principle, break the Heisenberg-Gabor limit [A simple explanation why the Lasso-type analysis canovercome the Heisenberg-Gabor limit in a diﬀerent ﬁeld of science can be found in Herman and Strohmer (2009)].
In Kato and Uemura (2012), we applied the least absolute shrinkage and selection operator (Lasso), a form ofcompressed sensing to frequency analysis of astronomical time-series observations (also called “sparse Fourier”
analysis). This method has a resolution superior to Cohen’s class transformations and Laﬂer-Kinman’s classmethods, and is particularly useful in separating closely spaced multiple frequency components, such as theorbital signals and superhumps signals in SU UMa-type dwarf novae.
For studying frequency variations, two-dimensional frequency analysis (dynamic spectrum) is used. For
evenly sampled data, STFT is usually used (e.g. Still et al. 2010). Although Kato and Uemura (2012) presented
a code for single-epoch analysis, two-dimensional frequency analysis using Lasso has been employed and various
results have been published (Kato and Maehara 2013; Osaki and Kato 2013; Kato and Osaki 2013a,b; Osaki
and Kato 2014; Kato et al. 2014b; Ohshima et al. 2014; Kato et al. 2014a; Pavlenko et al. 2014; Kato et al.
2015b; Kato et al. 2015a; Kato et al. 2016; Niijima et al. 2021). In Osaki and Kato (2013, 2014), Lasso analysis
played an important role in demonstrating the conclusion that the thermal-tidal instability model (Osaki 1989)
is the only viable model to explain superoutbursts in SU UMa stars. In these series of papers, we have shownthat Lasso analysis is advantageous in several respects: (1) the resultant signal is very sharp, (2) the frequency
resolution is very high, and (3) the result is less aﬀected by uneven sampling than in Fourier analysis. The lastpoint is particularly important in real astronomical data, in which observations can usually be obtained unevenly
in time. Comparisons of spectrograms between STFT and Lasso-based Fourier analysis are shown in Osaki andKato (2013); Kato and Osaki (2013a); Kato and Maehara (2013). The advantage of Lasso or sparse modelingin Fourier-type application has been extended to super-resolution imaging with radio interferometry by Honma
et al. (2014).
Considering that two-dimensional frequency analysis using Lasso is useful not only in astronomy but also inother ﬁelds of science, I present a multi-purpose full code to calculate and draw two-dimensional Lasso spectral
analysis. The code was originally written for Kato and Maehara (2013), and has been modiﬁed here for moregeneral purposes.

2

Methods

The simpliﬁed mathematical basis of our Lasso analysis is given here. For a more complete treatment andcomparison with other frequency determination algorithms, please refer to Kato and Uemura (2012). This
mathematical basis is not always necessary to use the code, and readers can skip this section to understandactual applications and to learn how to use the code.
The mathematical formulation is based on Tanaka (2010). The data (in amplitude) have Y (ti ) with timesof observations at ti . The mean of Y is set to be zero. The observation can be expressed as a sum of signal (Ys )and random errors (n) :(1)
Yi = Y (ti ) = Ys (ti ) + n(ti ).
The signal is assumed to be composed of a sum of strictly periodic functions. Using sine and cosine Fouriercomponents, Ys can be expressed as :∑∑Ys (ti ) =aj cos(ωj ti ) +bj sin(ωj ti ),
(2)
j

j

where ω are frequencies and a and b are amplitudes. This equation can be rewritten as
y = (Y1 , · · · , YN )T(3)

x = (a1 , · · · , aM , b1 , · · · , bM )T ,

(4)

and
where N and M are number of observations and number of diﬀerent ω, respectively. A set of equations 1 and 2
can be rewritten as a form of
y = Ax + n(5)using a 2M × N observation matrix A composed of{cos(ωi tj ),Ai,j =sin(ωi−M tj ),(i ≤ M )(i > M ).

(6)

3
The vector x0 is what to be estimated.
∑N
In Lasso, “1-norm” is used, which is ||x||1 ≡ i=1 |xi |, and x̂LAR :x̂LAR = arg min(x1||y − Ax||2 + λ||x||1 ),2N(7)

where ||y − Ax||2 is the ordinary squared sum of residuals, is chosen, and λ||x||1 is the ℓ1 -norm penalty function
with a parameter λ having a value of λ ≥ 0. arg minx f (x) means the value of x which minimizes f (x). Thisestimate becomes identical with a least-squares estimation at λ = 0. It has been known that this ℓ1 regularization
provides a sparse solution (small number of non-zero elements in x0 ) in most cases (cf. Donoho 2006; Candes
and Tao 2006) and there is a fast algorithm, known as least angle regression (LAR) by Efron et al. (2004) tosolve this problem.
In this Lasso analysis, I used LAR implementation lars package by Efron et al. (2004) on R software1
combined with glmnet (generalized linear model via penalized maximum likelihood; Friedman et al. 2010) as awrapper.

3

Application to frequency analysis of variable stars

3.1

Usage of the code

Since previous applications in research of variable stars are already given in section 1, I only provide an example
to run a code. Assuming that readers have the data of KIC 8751494 (Kato and Maehara 2013) as a form of aplain text with BJD−2400000 and magnitudes ready to read with R, the analysis can be done with:source("lassobase.R")
dat <- read.table("datafile")x <- dat$V1 - mean(dat$V1)
le <- lowess(x,dat$V2,f=0.05)
res <- dat$V2 - le$ytab <- data.frame(V1=dat$V1,V2=res-mean(res))pgm <- getpergrmlasso2(tab,56106,56205,1/0.130,1/0.105,200,10,1)
The second part is for prewhitening the data using the locally-weighted polynomial regression (LOWESS: Cleveland 1979). This part needs to be modiﬁed depending on the complexity of the light curve. This is particularly
important for analyzing light curves of dwarf novae. In papers such as Osaki and Kato (2013), we used a verysmall smoother span (f in function lowess) and a small delta in the LOWESS regression to remove the global
patterns of outbursts in dwarf novae and processed the entire data as a whole. This prescription led to favorable
results for many dwarf novae despite that the amplitudes of the target signals (superhumps, orbital variations andothers) are greatly reduced. To conserve the amplitudes of the target signals, it is better to divide the entire light
curve to shorter segments so that the LOWESS regression requires a larger f (typically larger than 0.1) to obtain
suﬃciently prewhitened segments but with well-conserved target signals. These short segments are combined toobtain the prewhitened entire data for the two-dimensional Lasso analysis.
The second and third arguments for getpergrmlasso2 gives the time range, and the fourth and ﬁfth for thefrequency range (c/d in this case). The sixth argument 200 is the number of frequency bins, the seventh 10 is
the window length and the ﬁnal 1 is the shift value for the windows. The resultant pgm can be saved using savefunction in R for further analysis or graphical representation. The result pgm can be readily seen on graphicalwindow by
drawpgmlasso(pgm,-2.7,-13,-6,3)
The second argument −2.7 is log10 λ and the third and fourth −13 and −6 are the display range of the powerstrength (in logarithmic scale) and the ﬁnal 3 is time bins to smear (use 0 when no smearing is needed). The
R graphical window would look like ﬁgure 1. The λ value can be optimized by looking at the resultant twodimensional spectrogram by changing the gray scale (by changing the display range of power strength), examining
1 The R Foundation for Statistical Computing: <http://cran.r-project.org/>.

0.00
−0.15

−0.10

−0.05

amplitude

0.05

0.10

0.15

4

56140

56120

56140

56160

56180

56200

56160

56180

56200

8.5
8.0

Frequency (c/d)

9.0

9.5

56120

Time (d)Figure 1: Example of graphical representation of two-dimensional Lasso analysis on R using the code presented
in this paper. Kepler data of KIC 8751494 were used (see Kato and Maehara 2013). The signal at 8.74 c/d isthe orbital signal. The varying signal between 8.0 and 8.2 c/d is the superhump signal. See text how to produce
this ﬁgure.

5whether target frequencies are well visible above the background noise. The list of diﬀerent values of λ by lars isstored in the vector pgm$lst[[i]]$lassopow$lambda, where i is the bin number (1 ≤ i ≤ length(pgm$lst)),and one can ﬁnd an adequate range of log10 λ by using these values. One can also draw a spectrogram of a
speciﬁed element (n) within a series of diﬀerent λ given by lars by setting the argument loglambda to FALSE:
drawpgmlasso(pgm,n,-13,-6,3,loglambda=FALSE)(n must be an integer and 1 ≤ n ≤ length(pgm$lst[[i]]$lassopow$lambda)). By changing n, one can easilyﬁnd the optimal λ.
Although one can use the popular cross-validation technique to obtain λ expressing the most regularized
model (see Kato and Uemura 2012), this is not always optimal for actual applications. This is probably because
the model using sine and cosine functions is too simple for describing the real variation.

3.2

Application to the super-fast nova V1674 Her

I give another example of a Lasso two-dimensional power spectrum of the fastest known nova V1674 Her. This nova
was discovered by Seiji Ueda on 2021 June 12,2 and is the fastest ever recorded classical nova with t2 ∼ 1.2 d(Quimby et al. 2021). This nova was reported to show coherent pulses with a period of 501.428 s in Zwicky
Transient Facility (ZTF; Masci et al. 2019) data before the eruption (Mroz et al. 2021) and this period is
considered to reﬂect the spin period of white dwarf. In addition to this period, orbital modulations with aperiod of 0.15302(2) d was identiﬁed after the nova eruption (Shugarov and Afonina 2021a; Patterson et al. 2021;
Shugarov and Afonina 2021b). Patterson et al. (2021) reported a spin period of 501.52(2) s after the eruption.
This value is 0.09(2) s longer than the value before the eruption and is consider to be a consequence of the nova
eruption, either by an increase in the moment of inertia or by a loss of mass that carried away angular momentum(Drake et al. 2021).
I examined whether these signals can be seen in Lasso two-dimensional power spectra. I used the datareported to VSNET Collaboration (Kato et al. 2004) (observers: Tonny Vanmunster, Michael Richmond, Stephen
Brincat, Charles Galdies and Franz-Josef Hambsch). The analysis presented here is preliminary and is planned
to appear in VSNET Collaboration in prep. with more details. In the middle panel of ﬁgure 2, the signal of theorbital period appeared around BJD 2459420 (∼40 d after the eruption). There was some hint of modulations
with a frequency slightly lower than that of the orbital signal before BJD 2459410. This signal has a period
longer than the orbital period and may be some sort of superhumps, although the data were not suﬃcient todraw a conclusion. In the lower panel of ﬁgure 2, a spin frequency was detected below the one observed before theeruption. This has conﬁrmed the results by Drake et al. (2021). Although period variations of the order of 0.1–
0.3 s suggested from X-ray observations (Drake et al. 2021) were not apparent, the Lasso analysis suggests thatthe spin frequency had a tendency to return to the frequency before the nova eruption. This may be supportiveof an idea that the decrease in the spin frequency was the result of the increase in the moment of inertia due tothe nova explosion.

4

Application to frequency analysis of vocalization of birds

4.1

General introduction and usage of the code

In analyzing vocalizations of birds, two-dimensional Fourier power spectra (sonograms) have been widely used
(e.g. Marler and Slabbekoorn 2004). STFT is usually used to construct two-dimensional spectra. Both opensource, cross-platform audio software audacity3 , or software specially designed for analysis of vocalizations ofbirds Raven Lite4 can produce research-grade sonograms. On R, dynspec function in seewave package5 canproduce a two-dimensional STFT power spectrum both as a form of matrix and graphical representation. They
are suﬃcient to study sonograms of vocalizations of birds in many situations.
Although this method is well-established in bioacoustic studies of birds, time resolutions or frequency resolutions limited by the Heisenberg-Gabor limit (section 1; see also Beecher 1988) could become a problem especiallywhen the durations of sounds are short. In analysis of song features, frequency tracking or methods on focusing
2 <http://www.cbat.eps.harvard.edu/unconf/followups/J18573095+1653396.html>.
3 <https://www.audacityteam.org/>
4 <https://ravensoundsoftware.com/software/raven-lite/>.
5 <https://cran.r-project.org/web/packages/seewave/index.html>.

6

8

10

12

14

16

18
7.0

6.8

6.6

Porb

6.4

6.2

6.0

172.4

Pspin before eruption
172.3

172.2

172.1

59400

59420

59440

59460

59480

59500

59520

Figure 2: Lasso two-dimensional power spectrum of the nova V1674 Her. The data were from observations
reported to VSNET Collaboration (Kato et al. 2004). (Upper) Light curve. The data were binned to 0.1 d.
(Middle) Lasso power spectrum around the orbital period The orbital signal (6.54 c/d, corresponding to a periodof 0.1530 d) has been present since BJD 2459420. (Lower) Lasso power spectrum around the spin period. The
window length and shift value for computing these spectra are 20 d and 1 d, respectively.

7on the main frequency have been applied to overcome this diﬃculty (Gall et al. 2012; Stowell and Plumbley 2012,2014). Such methods, however, can have diﬃculties in handling sounds with multiple frequency components.
Furthermore, it has been demonstrated that ultrafast (over 200 Hz) movements of vocal muscles control birdsong
in certain species (Elemans et al. 2008). In such vocalizations, the duration of each note is about 1 msec. Even
the optimal selection of parameters can give a frequency resolution of 1 kHz by the traditional STFT method,and it may not be adequate for tracing the variation of frequencies within individual notes. Elemans et al. (2008)indeed employed the sound strength as a tracer of vocalization, but if we could trace the frequency variation, it
would provide more information about the characteristic of notes and mechanisms to produce and control thesesounds.
In order to analyze vocalizations of birds by the present code, the sound data either in linear pulse code
modulated audio (LPCM) format or conventional compressed format such as MP3 need to be read just liketime-series data for variable stars. For reading LPCM wav format ﬁles, readWave function in tuneR package6
is useful. The actual code would be like:library(tuneR)
source("lassobase.R")s <- readWave("wren.wav",from=0,to=2,units="seconds")
y <- s@left
rate <- s@samp.rate
y <- y - mean(y)x <- (1:length(y))/rate * 1000d <- data.frame(V1=x,V2=y)pgm <- getpergrmlasso2(d,1390,1425,2,16,100,2,0.4)
In this example, I used a call of an Eurasian Wren (Troglodytes troglodytes) recorded locally in Japan using a
SONY recorder ICD-SX1000 in the LPCM mode (96 kHz/24-bit). The usage of getpergrmlasso2 would be
self-evident when compared to section 3. In sound analysis, the times are given in msec and frequencies are
given in kHz. In order to obtain a graphical representation drawpgmlasso is used as in analysis of variable stars.
In sound analysis, the amplitudes are usually much larger than in variable stars (of course, depending on thedeﬁnition of the amplitudes) and λ and the display range have larger values. In this case,drawpgmlasso(pgm,2.7,7,12,3,xlab="Time (msec)",ylab="Frequency (kHz)")yielded a sonogram of optimal quality. The parameters xlab and ylab are necessary to override the default values
designed for study of variable stars.
For MP3 ﬁles, readMP3 in tuneR package can be used instead:library(tuneR)
source("lassobase.R")s <- readMP3("sound.mp3")y <- s@left
rate <- s@samp.rate
y <- y - mean(y)x <- (1:length(y))/rate * 1000d <- data.frame(V1=x,V2=y)and the usage of drawpgmlasso is the same as the case as with LPCM recordings. I must note, however,the quality of MP3 ﬁles are highly variable depending on recording devices and the MP3 compression quality.
Although the use of LPCM recordings is highly recommended, MP3 ﬁles recorded by suitable recorders withhigh-quality compressions (such as at 320 kbps) can be used in frequency analysis usually without a problem.
In the following subsections, I describe selected results of analysis of vocalizations of birds by this code. They
may be useful not only to biologists but also to astronomers to interpret the features in Lasso spectrograms. The
6 <https://cran.r-project.org/web/packages/tuneR/index.html>.

8data used for this analysis were either recorded by myself in the LPCM mode or from xeno-canto public archive7(recording numbers are given after XC). The current archive of xeno-canto consists of MP3 recordings and I used
relatively high-quality ones.

4.2

Birdsong with frequency modulations

I used a song of a Common Chiﬀchaﬀ (Phylloscopus collybita) (XC214244) for a comparison between STFT
and Lasso spectrograms (ﬁgure 3). This song contains only one fundamental tone and shows strong frequency
modulations. Since birdsong usually has a broader range of frequency modulations than in many astronomicalsignals, the advantage of Lasso analysis in high frequency resolution is less apparent than in a similar comparisonin Kato and Maehara (2013). The frequencies, however, were detected as much sharper signals than in the STFT
analysis. This result is as sharp as those by frequency extraction methods reported in Stowell and Plumbley(2014). Although the second harmonic (ﬁrst overtone) was detected in both analyses, the signal was more sharply
detected in the Lasso analysis. The strength of the second harmonic tends to increase when the frequency of thefundamental tone decreases. This feature is less apparent in the STFT analysis. The comparison suggests thatboth STFT and Lasso methods can be equally used in analysis of song with relatively slow frequency modulation.
The diﬀerence may be more striking when more rapid frequency modulations are involved.
This result for the Lasso analysis can be reproduced by the following code:
library(tuneR)
source("lassobase.R")s <- readMP3("XC214244.mp3")y <- s@left
rate <- s@samp.rate
y <- y - mean(y)x <- (1:length(y))/rate * 1000d <- data.frame(V1=x,V2=y)pgm <- getpergrmlasso2(d,2650,2830,2,12,100,4,0.8)
drawpgmlasso(pgm,0,0,8,5,xlab="Time (msec)",ylab="Frequency (kHz)")

4.3

Analysis of a call of a Crested Honey Buzzard

In this case, I deal with a call with structures. The Crested Honey Buzzard (Pernis ptilorhynchus) is a diurnal
raptor species. The call is described as “high, four-note whistle wee hey wee hey or pii-yoo pii-ee, or a whistledscream kleeeur during breeding season” (Brazil 2009).8 Kabaya and Matsuda (2001) described: “The voice isaround 3 kHz. There is a short note pyu initially, then a prolonged iiii- voice followed by a slightly lower voice.
The entire voice has a falling intonation.” (translated from Japanese by myself).
Since the sonogram in Kabaya and Matsuda (2001) indicated the presence of a prolonged voice with aconstant frequency lasting 1 s overlapped on the rest of the notes, I examined how a Lasso two-dimensional
power spectrum can resolve these features. The voice was recorded by myself as in the case of the Eurasian
Wren. The result is shown in ﬁgure 4. The prolonged voice at 2.5 kHz is also seen in this spectrogram up to
1000 msec. Since this feature is seen in many other recordings of this species, it would be interesting to consider
the possibilities: (1) the bird continuously issued this frequency while issuing other frequencies, or (2) it is anechoed signal of the initial note or (3) the vibration somehow persisted in the body of the bird. Although echoing
is a promising explanation, this recording was obtained in an open space and it seemed to me diﬃcult to consider
a structure causing echoes within 150 m (the distance required from the duration of the signal). If the possibility(1) is the case, this phenomenon is called polyphony and it suggests that the bird controls a pair of muscles ofthe syrinx (vocal organ of birds) diﬀerently. Such a case is rarely documented in non-passerine birds (songbirds)
and would be worth studying in more detail.
7 <http://www.xeno-canto.org/>.
8 One might see how diﬃcult it is to describe vocalizations of a bird by human language!

9

1.0

0.5

0.0

−0.5

12

10

8

6

4

2
12

10

8

6

4

2
2660

2680

2700

2720

2740

2760

2780

2800

2820

Figure 3: Comparison of Lasso and STFT two-dimensional power spectra of a Common Chiﬀchaﬀ Phylloscopus
collybita (XC214244 in xeno-canto). The time is in msec from the start of the audio ﬁle. (Upper) Amplitudeof the signal. (Middle) STFT two-dimensional power spectrum. 100 frequency bins and an 8 msec window size
shifted by 0.8 msec were used. The window size was optimized to obtain the best resolution. (Lower) Lasso
two-dimensional power spectrum. 100 frequency bins and a 4 msec window size shifted by 0.8 msec were used.
The signal is resolved more sharply than in the STFT analysis.

10

1.0

0.5

0.0

−0.5

8

6

4

2

200

300

400

500

600

700

800

900

1000

Figure 4: Lasso two-dimensional power spectrum of a call of a Crested Honey Buzzard (Pernis ptilorhynchus).
(Upper) Amplitude of the signal. (Lower) Lasso two-dimensional power spectrum. 200 frequency bins and a
20 msec window size shifted by 3 msec were used. The complex structure of the vocalization is well depicted
together with the prolonged component around 2.5 kHz.

11
The entire sonogram is characterized by the variable fundamental frequency and the second (sometimes
the third) harmonic. Between 370 and 430 msec, the voice strength became smaller but with a higher (around3.8 kHz) frequency with multiple frequency components. The entire sequence appears to match the descriptionof a four-note whistle by Brazil (2009), although many other notes of the same species are not as complex as this
one.

4.4

Analysis of call of an Eurasian Wren

In this case, I deal with a short call with complex structures. A comparison of Lasso and STFT spectra of acall of an Eurasian Wren is shown in ﬁgure 5. This call is typical for this species, and I call this type of call“usual call”. This call is short and rich in rapidly varying signals. Although the STFT analysis, which has been
optimized for resolutions, showed the signal as a barely resolved band with rapid frequency modulations, theLasso analysis clearly showed that this signal is composed of rapidly varying closely separated frequencies. Thesefrequencies varied in parallel.
Such a pattern of frequencies can be interpreted in two ways: (a) They are adjacent higher harmonics. (b)They are the result of amplitude variation of a single frequency. Consider the following equation in which f andfm are the original frequency and the frequency of amplitude modulation, respectively (usually fm ≪ f ). The
relative amplitude of the modulation is given as a.
a(8)signal = {1 + a sin(2πfm t)} sin(2πf t) = sin(2πf t) + [cos{2π(f − fm )t} − cos{2π(f + fm )t}]2
In the resultant power spectrum, two close frequencies f ± fm appear in addition to the central frequency f .
I consider the interpretation (b) more likely since the amplitudes showed periodic modulations directly visible
to the eyes (upper panel of ﬁgure 5), and both the central and pararell frequencies are present and the strengthsof the pararell frequencies are almost the same. By numerical experiments with artiﬁcially modulated data, a
90% modulation in the amplitude is required to reproduce the pattern in the spectrogram of this recording. This
is in agreement with the actual amplitude variation.
In variable stars, two closely separated frequencies can cause a beat phenomenon. In such cases, however,the strengths of the signals are diﬀerent and the resultant amplitude modulations due to the beat phenomenon
are smaller than in the present case. Frequency analysis of variable stars with amplitude modulations requires
attention in interpreting the spectrogram as in the present case.

4.5

Diﬀerent types of calls in the Eurasian Wren

Since many avian species have rich vocabularies in calls, the characteristics of a call of the Eurasian Wren
described in subsection 4.4 may not be applicable to diﬀerent types of calls of the same species. In ﬁgure 6, I
show two samples (XC192853, XC192854) of the alarm call and one sample (XC194033) of the trill-type call (callcomposed of quickly repeated short notes) of the same species in Japan. The same duration (35 msec) is used
for analysis. In the alarm call, the initial part of the note is characterized by a single rising frequency in contrast
to the usual call. In the trill-type call, the structure is simple without closely separated multiple frequencies.
The frequency of the trill-type call more quickly decreases within a single note than in the usual call. We cansee that closely separated multiple frequencies (probably reﬂecting amplitude modulation) are not always seen inthis species.

4.6

Comparison between subspecies of the Eurasian Wren

I have made a comparative analysis of subspecies using xeno-canto database to see whether the spectral characteristics I found (subsection 4.4) are unique to the Japanese subspecies fumigatus or are useful in broadly identifying
the species. Out of 22 subspecies listed in xeno-canto archive, I analyzed indigenus, zetlandicus, szetschuanus
and taivanus.
The results are shown in ﬁgure 7. Although one of them (XC23363) is labeled as the alarm call, all thesounds have characteristics similar to the usual call recorded in Japan in that they show closely separated
multiple components with decreasing frequencies. Since these subspecies are widespread around the Palearctic,
it is highly likely this frequency structure is common to other subspecies of the Eurasian Wren, and can be used
for species identiﬁcation. The result of XC23363 (szetschuanus) even more clearly illustrates the presence ofequally spaced multiple frequencies than in Japanese recordings. The overall analysis suggests that the frequency
structures of calls in the Eurasian Wren are highly conserved between diﬀerent subspecies.

12

1.0

0.5

0.0

−0.5

15

10

5

15

10

5

1390

1395

1400

1405

1410

1415

1420

1425

Figure 5: Comparison of Lasso and STFT two-dimensional power spectrum of an Eurasian Wren. The time is in
msec from the start of the audio ﬁle. (Upper) Amplitude of the signal. (Middle) STFT two-dimensional power
spectrum. 100 frequency bins and a 4 msec window size shifted by 0.4 msec were used. The window size was
optimized to obtain the best resolution. (Lower) Lasso two-dimensional power spectrum. 100 frequency bins anda 2 msec window size shifted by 0.4 msec were used. Around the strongest part of the sound, the signal is split
into three frequencies in Lasso analysis.

13(comp)XC192853

1.0

1.0

0.5

0.5

0.0

0.0

−0.5

−0.5

15

15

10

10

5

5

1390 1395 1400 1405 1410 1415 1420 1425

610

615

620

XC192854

630

635

640

645

310

315

320

XC194033

1.0

1.0

0.5

0.5

0.0

0.0

−0.5

−0.5

15

15

10

10

5

5

1255 1260 1265 1270 1275 1280 1285 1290

625

285

290

295

300

305

Figure 6: Lasso spectral analysis of diﬀerent types of calls of the Eurasian Wren in Japan. The initial one (comp)is the same as the one we used in subsection 4.4.

14

indigenus (XC66279)

zetlandicus (XC48210)

1.0

1.0

0.5

0.5

0.0

0.0

−0.5

−0.5

15

15

10

10

5

5

4040 4045 4050 4055 4060 4065 4070

3200 3205 3210 3215 3220 3225 3230

szetschuanus (XC23363)

taivanus (XC59754)

1.0

1.0

0.5

0.5

0.0

0.0

−0.5

−0.5

15

15

10

10

5

5

8275 8280 8285 8290 8295 8300 8305 8310

4700 4705 4710 4715 4720 4725 4730 4735

Figure 7: Lasso spectral analysis of diﬀerent subspecies of the Eurasian Wren.

15

5

Code (lassobase.R)

This code is also available at <http://www.kusastro.kyoto-u.ac.jp/∼tkato/pub/lassobase.R>.
# The first part is the same as in Kato and Uemura (2012)# Publications of the Astronomical Society of Japan, 64, 122library(lars)
library(glmnet)
seqfreq <- function(a,b,...) {return(1/seq(1/b,1/a,...))}
makematlasso <- function(d,p,ndiv) {nd <- length(p)m <- matrix(0,nrow(d),nd*ndiv)for (i in 1:nd) {ph <- ((d$V1/p[i]) %% 1)*pi*2for (j in 0:(ndiv-1)) {m[,i+nd*j] <- sin(ph+pi*j/ndiv)}}
return(m)}perlasso <- function(d,p,ndiv=2,alpha=1,
cv=FALSE) {nd <- length(p)mat <- makematlasso(d,p,ndiv)
y <- d$V2 - mean(d$V2)
m <- glmnet(mat,y,alpha=alpha)ndim <- m$dim[2]pow <- matrix(0,nd,ndim)
for (i in 1:ndim) {v <- m$beta[,i]for (j in 0:(ndiv-1)) {pow[,i] <- pow[,i] +v[(nd*j+1):(nd*(j+1))]^2}}
nmin <- NULL
gcv <- NULL
if (cv) {gcv <- cv.glmnet(mat,d$V2,alpha=alpha)minl <- gcv$lambda.min
nmin <- which.min(abs(m$lambda-gcv$lambda.min))}r <- list(pow=pow,p=p,lambda=m$lambda,
nmin=nmin,m=m,gcv=gcv,mat=mat)class(r) <- c("lassopow",class(r))
return(r)}plot.lassopow <- function(pow,n,...) {16
p <- pow$p
pw <- pow$pow
plot(p,pw[,n],typ="l",xlab="Period",ylab="Relative Power",...)}
minper <- function(pow,n,num=1) {p <- pow$p
pw <- pow$pow
pp <- numeric()
while (num > 0) {i <- which.max(pw[,n])pp <- c(pp,p[i])num <- num-1
pw[i,n] <- 0}
return(pp)}lassofit <- function(pow,x,n,ndiv=2) {d <- data.frame(V1=x,V2=numeric(length(x)))
m <- pow$mat
return(m %*% pow$m$beta[,n])}lassoexpect <- function(pow,x,n,ndiv=2) {d <- data.frame(V1=x,V2=numeric(length(x)))
m <- makematlasso(d,pow$p,ndiv)
return(m %*% pow$m$beta[,n])}# 2-D version
pergrmlasso <- function(d,fmin,fmax,div) {d$V1 <- d$V1 - mean(d$V1)
d$V2 <- d$V2 - mean(d$V2)
p <- seq(fmin,fmax,length=div+1)
pl <- perlasso(d,1/p)
las <- list(pow=pl$pow,p=pl$p,lambda=pl$lambda)
r <- list(period=p,lassopow=las)class(r) <- c("PergrmLasso",class(r))
return(r)}pergrm2dlasso2 <- function(d,fmin,fmax,div,len,step) {nd <- floor((diff(range(d$V1))-len)/step+1)
lst <- list()st <- min(d$V1)cat("bins="); cat(nd); cat("\n")for (i in 1:nd) {t <- subset(d,V1 >= st+step*(i-1) & V1 < st+len+step*(i-1))
if (nrow(t) > 10) {p <- pergrmlasso(t,fmin,fmax,div)} else {p <- NA}17lst[[i]] <- p}return(lst)}getpergrmlasso2 <- function(d,start,end,fmin,fmax,div,len,step) {# range should be expressed in frequencyorgdata <- dd <- data.frame(V1=d$V1,V2=d$V2)
f <- subset(d,V1 > start & V1 < end)r <- list(g=f,f=orgdata)g <- r$g
f <- r$f
dt <- 1lst <- pergrm2dlasso2(g,fmin,fmax,div,len,step)yy <- seq(fmin,fmax,length=div+1)
nd <- floor((diff(range(g$V1))-len)/step+1)
xx <- min(g$V1)+(0:(nd-1))*step+len/2
r <- list(lst=lst,d=orgdata,f=f,g=g,xx=xx,yy=yy,len=len,dt=dt)
class(r) <- c("Pergrm2lasso",class(r))
return(r)}drawpgmsub <- function(pgm,pp,powmin,powmax,xlab,ylab) {d <- pgm$d
xx <- pgm$xx
yy <- pgm$yy
dt <- pgm$dt
len <- pgm$len
freq <- pgm$freq
def.par <- par(no.readonly = TRUE)
mat <- matrix(c(1,1,2,2,2), 5, 1, byrow = TRUE)
layout(mat)par(mar=c(2,5,0,4))par(oma=c(1,0,4,0))par(cex.lab=1.4)
par(cex.axis=1.2)
xlim <- c(min(xx)-dt*len/2,max(xx)+dt*len/2)
plot(d$V1,d$V2,xlim=xlim,ylim=range(d$V2),typ="l",xaxs="i",pch=19,cex=0.3,xlab="",ylab="amplitude")par(mar=c(4,5,0,4))
image(x=xx,y=yy,z=pp,xlim=xlim,
zlim=c(powmin,powmax),col=gray.colors(100,start=1,end=0,gamma=2),
xlab=xlab,ylab=ylab)par(def.par)}drawpgmlasso <- function(pgm,n,powmin,powmax,dup=0,
xlab="Time (d)",ylab="Frequency (c/d)",loglambda=TRUE) {18
if (!inherits(pgm,"Pergrm2lasso")) {stop("only Pergrm2lasso class can be drawn")}
nn <- length(pgm$lst)
pp <- matrix(0,nn,length(pgm$yy))
zerovec <- rep(0,length(pgm$yy))
if (dup == 0) {dwin <- 1} else {dwin <- 10^(-abs(((-dup):dup)/dup*10))}for (i in 1:nn) {if (is.na(pgm$lst[[i]][1])) {pw <- zerovec} else {if (loglambda == FALSE) {pw <- pgm$lst[[i]]$lassopow$pow[,n]} else {ni <- which.min(abs(log10(pgm$lst[[i]]$lassopow$lambda)-n))pw <- pgm$lst[[i]]$lassopow$pow[,ni]}}for (j in (-dup):dup) {if (i+j > 0 && i+j <= nn) {pp[i+j,] <- pp[i+j,] + dwin[j+dup+1]*pw}}}
pp <- log10(pp)pp <- ifelse(pp > powmax,powmax,pp)
drawpgmsub(pgm,pp,powmin,powmax,xlab,ylab)}

Acknowledgements
This work was supported by JSPS KAKENHI Grant Number 21K03616. This work was originally initiatedunder the grant “Initiative for High-Dimensional Data-Driven Science through Deepening of Sparse Modeling”(2013–2018) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, and I
am grateful to Toshiyuki Tanaka, Shiro Ikeda and Masato Okada for comments on compressed sensing and fortheir encouragement in application of compressed sensing to bioacoustics. I am grateful to the observers whoreported observations of V1674 Her to VSNET. I am grateful to the xeno-canto team and recordists for providing
the database and recordings available to the public.

References
Beecher, M. D. (1988) Spectrographic analysis of animal vocalizations: implications of the “uncertainty principle”.
Bioacoustics 1, 187
Brazil, M. (2009) Birds of East Asia: China, Taiwan, Korea, Japan, and Russia (Princeton: Princeton University
Press)
Candes, E. J., & Tao, T. (2006) Near-optimal signal recovery from random projections: Universal encoding
strategies?. IEEE Trans. Inf. Theory 52, 5406

19
Cleveland, W. S. (1979) Robust locally weighted regression and smoothing scatterplots. J. Amer. Statist. Assoc.
74, 829
Deeming, T. J. (1975) Fourier analysis with unequally-spaced data. Ap&SS 36, 137
Donoho, D. L. (2006) Compressed sensing. IEEE Trans. Inf. Theory 52, 1289
Drake, J. J. et al. (2021) The remarkable spin-down and ultrafast outﬂows of the highly pulsed supersoft sourceof Nova Herculis 2021. ApJ 922, L42
Dworetsky, M. M. (1983) A period-ﬁnding method for sparse randomly spaced observations of ‘How long is apiece of string?. MNRAS 203, 917
Efron, B., Hastie, T., Johnstone, I., & Tibshirani, R. (2004) Least angle regression. Ann. Statist. 32, 407
Elemans, C. P. H., Mead, A. F., Rome, L. C., & Goller, F. (2008) Superfast vocal muscles control song production
in songbirds. PLoS ONE 3, e2581
Fernie, J. D. (1989) Uncertainties in period determinations. PASP 101, 225
Friedman, J., Hastie, T., & Tibshirani, R. (2010) Regularization paths for generalized linear models via coordinate
descent. J. Statistical Software 33, 1
Gall, M. D., Brierley, L. E., & Lucas, J. R. (2012) The sender-receiver matching hypothesis: Support from theperipheral coding of acoustic features in songbirds. Journal of Experimental Biology 215, 3742
Herman, M. A., & Strohmer, T. (2009) High-resolution radar via compressed sensing. IEEE Transactions on
Signal Processing 57, 2275
Hlawatsch, A., & Auger, F. (2008) Time-Frequency Analysis: Concepts and Methods (Hoboken: Wiley)
Honma, M., Akiyama, K., Uemura, M., & Ikeda, S. (2014) Super-resolution imaging with radio interferometry
using sparse modeling. PASJ 66, 95
Horne, J. H., & Baliunas, S. L. (1986) A prescription for period analysis of unevenly sampled time series. ApJ
302, 757
Kabaya, T., & Matsuda, M. (2001) The Songs & Calls of 420 Birds in Japan (Tokyo: Shogakukan)
Kato, T. et al. (2014a) Survey of period variations of superhumps in SU UMa-type dwarf novae. VI: The ﬁfth
year (2013–2014). PASJ 66, 90
Kato, T. et al. (2015a) Survey of period variations of superhumps in SU UMa-type dwarf novae. VII: The sixthyear (2014–2015). PASJ 67, 105
Kato, T. et al. (2014b) Survey of period variations of superhumps in SU UMa-type dwarf novae. V: The ﬁfth
year (2012–2013). PASJ 66, 30
Kato, T., Hambsch, F.-J., Oksanen, A., Starr, P., & Henden, A. (2015b) CC Sculptoris: Eclipsing SU UMa-type
intermediate polar. PASJ 67, 3
Kato, T. et al. (2016) RZ Leonis Minoris bridging between ER Ursae Majoris-type dwarf nova and nova-likesystem. PASJ 68, 107
Kato, T., & Maehara, H. (2013) Analysis of Kepler light curve of the novalike cataclysmic variable KIC 8751494.
PASJ 65, 76
Kato, T. et al. (2010) Survey of Period Variations of Superhumps in SU UMa-Type Dwarf Novae. II. The Second
Year (2009-2010). PASJ 62, 1525
Kato, T., & Osaki, Y. (2013a) Analysis of three SU UMa-type dwarf novae in the Kepler ﬁeld. PASJ 65, 97
Kato, T., & Osaki, Y. (2013b) KIC 7524178 – an SU UMa-type dwarf nova showing predominantly negative
superhumps throughout supercycle. PASJ 65, L13

20
Kato, T., & Uemura, M. (2012) Period analysis using the Least Absolute Shrinkage and Selection Operator(Lasso). PASJ 64, 122
Kato, T., Uemura, M., Ishioka, R., Nogami, D., Kunjaya, C., Baba, H., & Yamaoka, H. (2004) Variable Star
Network: World center for transient object astronomy and variable stars. PASJ 56, S1
Laﬂer, J., & Kinman, T. D. (1965) An RR Lyrae star survey with the Lick 20-inch Astrograph II. the calculationof RR Lyrae periods by electronic computer. ApJS 11, 216
Marler, P., & Slabbekoorn, H. (2004) Nature’s Music. The Science of Birdsong (Waltham: Academic Press,
Elsevier)
Masci, F.-J. et al. (2019) The Zwicky Transient Facility: Data processing, products, and archive. PASP 131,018003Mroz, P., Burdge, K., van Roestel, J., Prince, T., Kong, A. K. H., & Li, K.-L. (2021) An 8.4 min period in thearchival ZTF light curve of Nova Herculis 2021. Astron. Telegram 14720
Niijima, K. et al. (2021) Optical variability correlated with X-ray spectral transition in the black-hole transient
ASASSN-18ey = MAXI J1820+070. VSOLJ Variable Star Bull. 74, (arXiv:2107.03681)
Ohshima, T. et al. (2014) Study of negative and positive superhumps in ER Ursae Majoris. PASJ 66, 67
Osaki, Y. (1989) A model for the superoutburst phenomenon of SU Ursae Majoris stars. PASJ 41, 1005
Osaki, Y., & Kato, T. (2013) Study of superoutbursts and superhumps in SU UMa stars by the Kepler light
curves of V344 Lyrae and V1504 Cygni. PASJ 65, 95
Osaki, Y., & Kato, T. (2014) A further study of superoutbursts and superhumps in SU UMa stars by the Kepler
light curves of V1504 Cygni and V344 Lyrae. PASJ 66, 15
Patterson, J., Epstein-Martin, M., Vanmunster, T., & Kemp, J. (2021) The orbital and pulse periods of V1674
Herculis. Astron. Telegram 14856
Pavlenko, E., Kato, T., Sosnovskij, A. A., Andreev, M. V., Ohshima, T., Sklyanov, A. S., Bikmaev, I. F., &Galeev, A. I. (2014) Dwarf nova EZ Lyncis second visit to instability strip. PASJ 66, 113
Quimby, R. M., Shafter, A. W., & Corbett, H. (2021) The detailed light-curve evolution of V1674 Her (Nova Her2021). Research Notes of the American Astronom. Soc. 5, 160
Roberts, D. H., Lehar, J., & Dreher, J. W. (1987) Time series analysis with Clean – part one – derivation of aspectrum. AJ 93, 968
Scargle, J. D. (1982) Studies in astronomical time series analysis. II – statistical aspects of spectral analysis of
unevenly spaced data. ApJ 263, 835
Schwarzenberg-Czerny, A. (1989) On the advantage of using analysis of variance for period search. MNRAS 241,
153
Shugarov, S., & Afonina, M. (2021a) Detection of a possibly orbital period of the classical Nova V1674 Her.
Astron. Telegram 14835
Shugarov, S., & Afonina, M. (2021b) Photometric study of classical Nova V1674 Her. Perem. Zvezdy 41, 4Stellingwerf, R. F. (1978) Period determination using phase dispersion minimization. ApJ 224, 953
Still, M., Howell, S. B., Wood, M. A., Cannizzo, J. K., & Smale, A. P. (2010) Quiescent superhumps detected in
the dwarf nova V344 Lyrae by Kepler. ApJ 717, L113
Stowell, D., & Plumbley, M. D. (2012) Framewise heterodyne chirp analysis of birdsong. Proceedings of theEuropean Signal Processing Conference (EUSIPCO) pp 2694–2698
Stowell, D., & Plumbley, M. D. (2014) Large-scale analysis of frequency modulation in birdsong databases.
Methods in Ecology and Evolution 5, 901

21Tanaka, T. (2010) Mathematics of Compressed Sensing. IEICE Fundamentals Review 4, 39
Zechmeister, M., & Kürster, M. (2009) The generalised Lomb-Scargle periodogram. A new formalism for theﬂoating-mean and Keplerian periodograms. A&A 496, 577

VSOLJ
c/o Keiichi Saijo National Science Museum, Ueno-Park, Tokyo Japan
Editor Seiichiro Kiyota
e-mail: skiyotax@gmail.com