Sound analysis and synthesis

Description

seewave provides functions for analysing, manipulating, displaying, editing and synthesizing time waves (particularly sound). This package processes in particular time analysis (oscillograms and envelopes), spectral content, resonance quality factor, entropy, cross correlation and autocorrelation, zero-crossing, frequency coherence, dominant frequency, analytic signal, 2D and 3D spectrograms.

Details

Author(s)

Jerome Sueur <sueur@mnhn.fr>
Thierry Aubin
Caroline Simonis
Maintainer: Jerome Sueur <sueur@mnhn.fr>

Acoustic Complexity Index

Description

This function computes the Acoustic Complexity Index (ACI) as described in Pieretti et al. (2011)

Usage

ACI(wave, f, channel = 1, wl = 512, ovlp = 0,  wn = "hamming", flim = NULL, nbwindows = 1)

Arguments

Details

The function computes first a short-term Fourier transform and then the ACI index.
The function returns only the ACI total, ACI tot in Pieretti et al. (2010).
See the references for details on computation.

Value

A vector of length 1 returning the ACI total.

Note

Values returned were checked with the results provided by the add-on Soundscapemeter for the software Wavesurfer.

Author(s)

Laurent Lellouch, improved by Amandine Gasc and Morgane Papin

References

Pieretti N, Farina A, Morri FD (2011) A new methodology to infer the singing activity of an avian community: the Acoustic Complexity Index (ACI). Ecological Indicators, 11, 868-873.
Farina A, Pieretti N, Piccioli L (2011) The soundscape methodology for long-term bird monitoring: a Mediterranean Europe case-study. Ecological Informatics, 6, 354-363.

See Also

spectro, specflux

Examples

data(tico)
ACI(tico)
## dividing the sound sample into 4 windows of equal duration
ACI(tico, nbwindows=4)
## selection of a frequency band
ACI(tico, flim=c(2,6))

Acoustic Richness index

Description

This function computes the Acoustic Richness index based on M and Ht indices

Usage

AR(..., datatype = "objects", envt = "hil",
msmooth = NULL, ksmooth = NULL, ssmooth = NULL,
pattern = "[wav]$|[WAV]$|[mp3]$")

Arguments

Details

AR is ranked index based on the rank of the M and Ht indices obtained with the functions M and th respectively following:

AR = \frac{rank(M) \times rank(H_{t})}{n^2}

with

0 \leq AR \leq 1

Value

A data.frame with three columns (M, Ht, AR) and n columns, with n the number of objects (respectively files) used as input.

Note

As a ranked index, the results returned by AR strongly depends with the set of objects (respectively files) used as input. Comparaison between different data sets may be spurious. Computing AR on a set of a single object does not make any sense but is allowed.

Author(s)

Jerome Sueur and Marion Depraetere

References

Depraetere M, Pavoine S, Jiguet F, Gasc A, Duvail S, Sueur J (2012) Monitoring animal diversity using acoustic indices: implementation in a temperate woodland. Ecological Indicators, 13, 46-54.

See Also

M, th, env

Examples

## input as R objects
data(orni)
data(tico)
AR(orni, tico)
## give names to objects if you wish to have them as row names of the returned data.frame
AR(orni=orni, tico=tico)
## input as files stored in the working directory
## file names will be used as row names of the returned data.frame
## Not run: 
require(tuneR)
AR(getwd(), datatype="files")

## End(Not run)

Total entropy

Description

This function estimates the total entropy of a time wave.

Usage

H(wave, f, channel = 1, wl = 512, envt="hil", msmooth = NULL, ksmooth = NULL)

Arguments

Details

This function computes the product between the values obtained with sh and th functions.
This then gives a global (time and frequency) estimation of signal entropy.
The frequency mean spectrum and the amplitude envelope needed for computing respectively sh and th are automatically generated. They can be controlled through wl and smooth arguments respectively. See examples below and examples in sh and th for implications on the results.

Value

A single value varying between 0 and 1 is returned. The value has no unit.

Note

The entropy of a noisy signal will tend towards 1 whereas the entropy of a pure tone signal will tend towards 0.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Sueur, J., Pavoine, S., Hamerlynck, O. & Duvail, S. (2008) - Rapid acoustic survey for biodiversity appraisal. PLoS ONE, 3(12): e4065.

See Also

sh, th, csh

Examples

data(orni)
H(orni,f=22050)
# changing the spectral parameter (wl)
H(orni,f=22050,wl=1024)
# changing the temporal parameter (msmooth)
H(orni,f=22050,msmooth=c(20,0))

Median of the amplitude envelope

Description

This function computes an acoustic index based on the median of the amplitude envelope.

Usage

M(wave, f, channel = 1, envt = "hil", plot = FALSE, ...)

Arguments

Details

This amplitude index M is computed according to:

M = \bar{A}(t) \times 2^{1-depth}

with

0 \leq M \leq 1

where A(t) is the amplitude envelope and depth is the signal digitization depth in number of bits.

Value

A numeric vector of length 1 between 0 and 1, without unit.

Author(s)

Jerome Sueur and Marion Depraetere

References

Depraetere M, Pavoine S, Jiguet F, Gasc A, Duvail S, Sueur J (2012) Monitoring animal diversity using acoustic indices: implementation in a temperate woodland. Ecological Indicators, 13, 46-54.

See Also

env, AR

Examples

data(tico)
M(tico)
# smoothing the amplitude may change slightly the result
M(tico, msmooth=c(500,50), plot=TRUE)

Normalized Difference Soundscape Index

Description

This function computes the Normalized Difference Soundscape Index as described by Kasten et al. (2012).

Usage

NDSI(x, anthropophony = 1, biophony = 2:8, max = FALSE)

Arguments

Details

NDSI aims at estimating the level of anthropogenic disturbance on the soundscape by computing the ratio of human-generated (anthropophony) to biological (biophony) acoustic components found in field collected sound samples. In terms of frequency, the anthropophony is defined as the [1-2[ kHz frequency bin and the biophony as the [2-8[ kHz frequency bins of a soundscape frequency spectrum (see soundscapespec).

NDSI is computed according to:

NDSI = \frac{(biophony - anthropophony)}{(biophony + anthropophony)}

NDSI varies between -1 and +1, where +1 indicates a signal containing no anthropophony.

Value

A numeric vector of length 1 giving the NDSI value.

Author(s)

Jerome Sueur

References

Kasten, E.P., Gage, S.H., Fox, J. & Joo, W. (2012). The remote environmental assessment laboratory's acoustic library: an archive for studying soundscape ecology. Ecological Informatics, 12, 50-67.

See Also

soundscapespec, SAX, NDSI

Examples

## Note that 'tico' is not a soundscape recording...
data(tico)
spec <- soundscapespec(tico, plot=FALSE)
NDSI(spec)
NDSI(spec, max=TRUE)

Resonance quality factor of a frequency spectrum

Description

This function estimates the frequency pureness of a time wave by returning the resonant quality factor Q at a specific dB level.

Usage

Q(spec, f = NULL, level = -3, mel = FALSE, plot = TRUE, colval = "red",
cexval = 1, fontval = 1, flab = NULL,
alab = "Relative amplitude (dB)", type = "l", ...)

Arguments

Details

A high Q value indicates a highly resonant system.

Value

A list is returned with the following four items:

Note

This function is based on an linear interpolation of the spectrum so that the result should be considered as an estimation, not an exact measure.

Author(s)

Jerome Sueur, improved by Laurent Lellouch

See Also

spec, meanspec, corspec, fft.

Examples

# bird song
data(tico)
t<-spec(tico,f=22050,at=1.1,plot=FALSE,dB="max0")
op<-par(mfrow=c(2,1),las=1)
Q(t,type="l")
Q(t,type="l",xlim=c(3.8,4.2),ylim=c(-60,0))
title("zoom in")
par(op)
# cricket, changing the dB level
data(pellucens)
p<-spec(pellucens,f=11025,at=0.5,plot=FALSE,dB="max0")
op<-par(mfrow=c(3,1))
Q(p,type="l",xlim=c(1.8,2.6),ylim=c(-70,0))
title("level = - 3 (default value)",col.main="red")
Q(p,type="l",level=-6,
    xlim=c(1.8,2.6),ylim=c(-70,0),colval="blue")
title("level = - 6",col.main="blue")
Q(p,type="l",level=-9,
    xlim=c(1.8,2.6),ylim=c(-70,0),colval="green")
title("level = - 9",col.main="green")
par(op)

Symbolic Aggregate approXimation

Description

This function converts a numeric times seris into a series of letters with a specific length and alphabet.

Usage

SAX(x, alphabet_size, PAA_number,
breakpoints = "gaussian", collapse = NULL)

Arguments

Details

The SAX method has been developed to reduce the dimensionality of a numerical series into a short chain of characters. SAX follows a two-step process: (1) Piecewise Aggregate Approximation (PAA) and (2) conversion a PAA sequence into a series of letters.

PAA consists in a Z-normalisation, a segmentation of the series of length n into w segments, and the computation of each segment average.

The conversion of the PAA into a series of letters is achieved by attributing with equiprobability each value of the PAA to a letter in reference to a Gaussian distribution. This process therefore assumes that the distribution of the numeric series x follows a Gaussian distribution. To relax the constraints of normality we here added the possibility to directly work on the quantiles of the original data distribution or to specify particular breakpoints along the distribution of x. See the examples.

Value

A character vector of length (when collapse is NULL) or number of character (when collapse is not NULL) corresponding to PAA_number argument.

Note

SAX has been used recently to search similar times series in a soundcape data base (Kasten et al., 2012).

Author(s)

Laurent Lellouch. An improvement added by Pavel Senin.

References

Kasten, E.P., Gage, S.H., Fox, J. & Joo, W. (2012). The remote environmental assessment laboratory's acoustic library: an archive for studying soundscape ecology. Ecological Informatics, 12, 50 - 67.

Lin, J., Keogh, E., Lonardi, S., Chiu, B., June (2003). A symbolic representation of time series with implications for streaming algorithms. Proceedings of the 8th ACM SIGMOD Workshop on Research Issues in Data Mining and Knowledge Discovery. San Diego, California, USA.

See Also

discrets, symba, soundscapespec

Examples

data(tico)
spec <- soundscapespec(tico, plot=FALSE)[,2]
SAX(spec, alphabet = 5, PAA = 10)

# change breakpoints
SAX(spec,  alphabet = 5, PAA = 10, breakpoints="quantiles")
SAX(spec,  alphabet = 5, PAA = 10, breakpoints=c(0, 0.5, 0.75, 1))
SAX(spec,  alphabet = 5, PAA = 10, breakpoints=c(0, 0.33, 0.66, 1))

# different output formats
SAX(spec,  alphabet = 5, PAA = 10, collapse="")
SAX(spec,  alphabet = 5, PAA = 10, collapse="-")

normalized Time and Frequency Second Derivative

Description

This function computes the normalized Time and Frequency Second Derivative as described by Aumond et al. (2017).

Usage

TFSD(wave, f, channel = 1, ovlp = 0,  wn = "hamming", flim = c(2,6), nbwindows = 1)

Arguments

Details

The TFSD aims at estimating the time of presence of avian or human vocalizations within a sound environment. It calculates the variation in time and frequency of a signal around frequencies of interest, normalized by the spectral time variation of a signal as a whole.

Warning, this index was initially developed to work from a third octave spectrogram with a time sampling of 125 ms.

TFSD is computed according to formulation in reference.

The higher the TFSD varies between 0 and 1, the greater the temporal presence of avian or human vocalizations. With the default configuration, a TFSD > 0.3 indicates a very important presence time of the vocalizations in the signal. The TFSD is always greater than 0.

Value

A numeric vector of length nbwindows giving the TFSD values.

Author(s)

Pierre Aumond, Guillaume Corbeau

References

Aumond, P., Can, A., De Coensel, B., Botteldooren, D., Ribeiro, C., & Lavandier, C. (2017). Modeling soundscape pleasantness using perceptual assessments and acoustic measurements along paths in urban context. Acta Acustica united with Acustica, 12, 50-67.

Gontier, F., Lavandier, C., Aumond, P., Lagrange, M., & Petiot, J. F. (2019). Estimation of the perceived time of presence of sources in urban acoustic environments using deep learning techniques. Acta Acustica united with Acustica, 105(6), 1053-1066.

See Also

ACI, NDSI

Examples

## Note that 'tico' is not a soundscape recording...
data(tico)
TFSD(tico)
## dividing the sound sample into 4 windows of equal duration
TFSD(tico, nbwindows=4)
## selection of a frequency band
TFSD(tico, flim=c(2,6))

Teager-Kaiser energy tracking operator

Description

This function computes the Teager-Kaiser energy operator.

Usage

TKEO(wave, f, channel = 1, m = 1, M = 1, plot = TRUE,
xlab = "Time (s)", ylab = "Energy",
type = "l", bty = "l", ...)

Arguments

Details

The discrete version of the Teager-Kaiser operator is computed according to:

y_{n} = x_{n}^{2/m} - (x_{n-M} \times x_{n+M})^{1/m}

,
with m the exponent parameter and M the lag parameter which both are usually equal to 1 for a conventional operator.
The Teaser-Kaiser operator can be used to track amplitude modulations (AM) and/or frequency modulations (FM).
See examples.

Value

This function returns a two-column matrix, the first column is time and the second column includes the successive energy values.
m/2 NA values are added at the start and end of the vector.

Author(s)

Jerome Sueur

References

Kvedalen, E. (2003). Signal processing using the Teager Energy Operator and other nonlinear operators. University of Oslo, Department of Informatics, PhD Thesis, x + 100 p.

See Also

env, ifreq.

Examples

op <- par(mfrow=c(2,1))

## sinusoid AM 
s1 <- synth(f=8000, d=0.1, cf=200, am=c(100,10), output="Wave") 
oscillo(s1)
TKEO(s1)
## linear AM decrease
s2 <- synth(f=8000, d=0.1, cf=200, shape="decr", output="Wave") 
oscillo(s2)
TKEO(s2)
## sinusoid FM
s3 <- synth(f=8000, d=0.1, cf=200, fm=c(150,50,0,0,0), output="Wave") 
oscillo(s3)
TKEO(s3)
## linear FM increase
s4 <- synth(f=8000, d=0.1, cf=200, fm=c(0,0,600,0,0), output="Wave") 
oscillo(s4)
TKEO(s4)
## AM and FM
s5 <- synth(f=8000, d=0.1, cf=200, am=c(100,10), fm=c(150,50,0,0,0), output="Wave")
oscillo(s5)
TKEO(s5)
par(op)

Statistics on time and frequency STFT contours

Description

This function returns statistics based on STFT time and frequency contours.

Usage

acoustat(wave, f, channel = 1, wl = 512, ovlp = 0, wn = "hanning",
tlim = NULL, flim = NULL,
aggregate = sum, fraction = 90,
plot = TRUE, type = "l", ...)

Arguments

Details

The principle of acoustat is as follows:

1. Compute the short-term Fourier transform (STFT) with usual parameters (wl for window length, ovlp for overlap of successive windows, and wn for the name of window shape).

2. This results in a time * frequency matrix.

3. Compute an aggregation function (specified with the argument aggregate set by default to sum) accross rows and columns of time * frequency matrix.)

4. This results in two components: (i) the time contour, and (ii) the frequency contour.

5. Each contour is considered as a probability mass function (PMF) and transformed into a cumulated distribution function (CDF).

6. Measures are extracted from each CDF: median (M), initial percentile (P1) value, terminal percentile (P2) value, interpercentile range (IPR). P1, P2 and IPR are defined using a fraction parameter (fraction) that sets the percent of the contour amplitude to be captured by the initial and terminal percentile values. A fraction of 50% would result in the familiar quartiles and interquartile range. An energy fraction of 80% would return the 10th and 90th percentile values, and the width of the range in between.

Value

The function returns a list with 10 items:

Note

acoustat was originally developped in Matlab language by Kurt Fristrup and XXXX Watkins (1992) .
The R function was kindly checked by Kurt Fristrup.

Author(s)

Jerome Sueur

References

Fristrup, K. M. and Watkins, W. A. 1992. Characterizing acoustic features of marine animal sounds. Woods Hole Oceanographic Institution Technical Report WHOI-92-04.

See Also

meanspec, specprop

Examples

data(tico)
note <- cutw(tico, from=0.5, to=0.9, output="Wave")
## default setting
acoustat(note)
## change the percentile fraction
acoustat(note, fraction=50)
## change the STFT parameters
acoustat(note, wl=1024, ovlp=80)
## change the function to compute the aggregate contours
## standard deviation instead of sum   
acoustat(note, aggregate=sd)
## direct time and frequency selection     
acoustat(tico, tlim=c(0.5,0.9), flim=c(3,6))
## some useless graphical changes
acoustat(note, type="o", col="blue") 

Add or insert a silence section

Description

Add or insert a silence section to a time wave.

Usage

addsilw(wave, f, channel = 1, at = "end", choose = FALSE, d = NULL,
plot = FALSE, output = "matrix", ...)

Arguments

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

oscillo, cutw,deletew, fadew,pastew, mutew,revw, zapsilw

Amplitude filter

Description

This function deletes all signal which amplitude is below a selected threshold.

Usage

afilter(wave, f, channel = 1, threshold = 5, plot = TRUE,
listen = FALSE, output = "matrix", ...)

Arguments

Details

The threshold value is in % relative to the maximal value of wave. Signal inferior to this value is clipped.

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Note

This function is used as an argument (threshold) in the following functions: autoc, csh, dfreq, timer and zc.

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

oscillo

Examples

data(orni)
op<-par(mfrow=c(2,1))
afilter(orni,f=22050)
title(main = "threshold level = 5")
afilter(orni,f=22050,threshold=0.5,colwave="blue")
title(main = "threshold level = 0.5")
par(op)

Water tank minimum resonant and cutoff frequencies

Description

This function computes the resonant and cutoff frequencies when recording in a given aquarium according to the criteria explained in Akamatsu et al. (2002)

Usage

akamatsu(Lx, Ly, Lz, mode = c(1,1,1),
         c = 148000,  plot = FALSE, xlab = "Frequency (kHz)",
         ylab = "Attenuation distance (cm)", ...)

Arguments

Details

From Akamatsu et al. (2002):

1. Resonant frequency

The calculated resonant frequencies of a rectangular glass tank with the dimension of Lx , Ly , and Lz (in centimeters) can be described by the following equation:

f^{rectangular}_{lmn} = \frac{c}{2} \times \sqrt{\left(\frac{l}{L_{x}}\right)^2 + \left(\frac{m}{L_{y}}\right)^2 + \left(\frac{n}{L_{z}}\right)^2}

where c is the sound velocity (cm/s) and each l, m, n reprents an integer, and the combination of these paramameters designates the 'mode number'. The mode (1, 1, 1) represents the resonance wave of minimum frequency. The mode (2, 1, 1) represents one of the higher order of resonant component and has additional node of the soundpressure level at the middle of the X axis, i.e., Lx/2.


2. Cutoff frequency

The cutoff frequency can be calculated as follows:

f^{rectangular}_{cutoff} = \frac{c}{2} \times \sqrt{ \left(\frac{1}{L_{y}}\right)^2 + \left(\frac{1}{L_{z}}\right)^2}

3. Attenuation distance

The theoretical attenuation distance D can be expressed in function of the cutoff frequency and the projected frequency following:

D^{rectangular}(f) = 2 \times log_{10} \times \frac{c}{4 \pi f^{rectangular}_{cutoff}} \times \frac{1}{\sqrt{1-\left(\frac{f}{f^{rectangular}_{cutoff}}\right)^2}}

Value

A list of two items:

Author(s)

Camille Desjonqueres

References

Akamatsu T, Okumura T, Novarini N, Yan HY (2002) Emprical refinements applicable to the recording of fish sounds in small tanks. Journal of the Acoustical Society of America, 112, 3073-3082.

Examples

akamatsu(60, 30, 40)

Amplitude modulation analysis of a time wave

Description

This function computes the Fourier analysis of a time wave envelope. This allows to detect periodicity, in particular those generated by amplitude modulations.

Usage

ama(wave, f, channel = 1, envt = "hil", wl = 512, plot = TRUE, type = "l", ...)

Arguments

Details

This function is based on env and meanspec.
The envelope of wave is first computed and the spectrum of this envelope is then processed. All env and meanspec arguments can be set up. Be sure to set up wl large enough if you want to detect low amplitude modulation periodicity.

Value

If plot is FALSE, ama returns a numeric vector corresponding to the computed spectrum. If peaks is not NULL, ama returns a list with two elements:

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

env, fma, meanspec

Examples

data(orni)
# detection of the main amplitude modulation in a cicada song:
# one with a 0.258 kHz frequency (due to pulses in the echemes)
# one with a 2.369 kHz frequency (fundamental frequency)
ama(orni,f=22050,wl=1024)
# these amplitude modulations can be identify with a cursor:
ama(orni,f=22050,wl=1024,identify=TRUE)

Generate sound intensity attenuation data

Description

This function generates dB data following theoretical spherical attenuation of sound.

Usage

attenuation(lref, dref = 1, dstop, n, plot = TRUE,
xlab = "Distance (m)", ylab = "dB", type = "l", ...)

Arguments

Value

If plot is FALSE return a numeric vector with the data generated.

Note

Sound attenuation in a free, unbounded medium behaves in accordance with the inverse square law. attenuation generates data following this rule from a reference point where sound intensity level (SIL) or sound pressure level (SPL) is known. Such theoretical data can be compared with experimental data collected in a real environment.

Author(s)

Jerome Sueur

References

Hartmann, W. M. 1998 Signals, sound and sensation. New York: Springer.

See Also

convSPL, moredB

Examples

# theoretical attenuation up to 150 m of a 100 dB/1m sound source
attenuation(lref=100, dref=1, dstop=150, n=200)

Reading and interpreting Audiomoth file name

Description

This function reads and decomposes the files names generated by an Audiomoth device, audio digal recorders produced by the society Open Acoustic Devices.

Usage

audiomoth(x, tz = "")

Arguments

Details

The digital recorder Audiomoth produced by Open Acoustic Devices (https://www.openacousticdevices.info/) generates .wav files which names contains information about the time of recording. The information is encoded in hexadecimal (e.g. "5E9089F0"). The function audiomoth decodes this information so that time of recording can be retrieved in numeric or time format.

Value

The function returns a data.frame with the following columns:

Note

For the time zone see the 607 time zone names stored in OlsonNames.
The file names of Audiomoth may change with time. There is no guarantee that the function will be updated on time.

Author(s)

Jerome Sueur

References

See Open Acoustic Devices website for details regarding the Audiomoth: https://www.openacousticdevices.info/.

See Also

audiomoth.rename, as.POSIXct, OlsonNames, songmeter

Examples

## HEXADECIMAL EXAMPLES (OLD FORMAT)
## recording done on Friday 10 April 2020 16:54:44 UTC
## computer time zone (local time, Europe, Paris for the test)
audiomoth("5E90A4D4.WAV")
## UTC
audiomoth("5E90A4D4.WAV", tz="UTC")
## GMT (= UTC as UTC and GMT are synonyms)    
audiomoth("5E90A4D4.WAV", tz="GMT")
## UTC -2
audiomoth("5E90A4D4.WAV", tz="Etc/GMT-2")
## in Asia, Japan
audiomoth("5E90A4D4.WAV", tz="Japan")
## in South-America, Cayenne
audiomoth("5E90A4D4.WAV", tz="America/Cayenne")  
## several files
filenames <- c("5E914ED0.WAV", "5E915128.WAV",
"5E915380.WAV", "5E9155D8.WAV", "5E915830.WAV",
"5E915A88.WAV", "5E915CE0.WAV", "5E915F38.WAV",
"5E916190.WAV", "5E9163E8.WAV")
audiomoth(filenames)
## YYYYMMDD_HHMMSS.WAV FORMAT (ACTUAL FORMAT)
## single file
audiomoth("20230715_150000.wav")
## several files
filenames <- c("20230715_150000.wav", "20230715_151500.wav",
"20230715_153000.wav", "20230715_154500.wav")
audiomoth(filenames)

Rename audiomoth files in a readable format

Description

This function renames or copies files created with an Audiomoth device in a readable format including the data and time of recording.

Usage

audiomoth.rename(dir, overwrite = FALSE, tz = "", prefix = "")

Arguments

Details

The format of the new file names follows the format of the SongMeter SM2/SM4 deveices: PREFIX_YYYYMMDD_HHMMSS.wav.

Value

1 logical vector indicating which operation succeeded for each of the files attempted.

Note

Be careful if you overwrite the files.

Author(s)

Jerome Sueur

See Also

audiomoth, songmeter

Short-term autocorrelation of a time wave

Description

This function returns the fundamental frequency of a harmonic time wave. This is achieved by computing a correlation of the signal with itself after a time delay.

Usage

autoc(wave, f, channel = 1, wl = 512, fmin, fmax, threshold = NULL, plot = TRUE,
xlab = "Time (s)", ylab = "Frequency (kHz)", ylim = c(0, f/2000), pb =
FALSE, ...)

Arguments

Details

'fmin' and 'fmax' can help by reducing computing time but can also produce less accurate results.

Value

When plot is FALSE, autoc returns a two-column matrix, the first column corresponding to time in seconds (x-axis) and the second column corresponding to to fundamental frequency in kHz (y-axis).
NA corresponds to pause sections in wave (see threshold).

Author(s)

Jerome Sueur sueur@mnhn.fr and Thierry Aubin thierry.aubin@u-psud.fr

References

Hopp, S. L., Owren, M. J. and Evans, C. S. (Eds) 1998. Animal acoustic communication. Springer, Berlin, Heidelberg.

See Also

ceps, acf

Examples

data(sheep)
# fundamental frequency of a sheep
res <- autoc(sheep, f=8000, threshold=5, fmin=100, fmax=700, plot=FALSE)
spectro(sheep, f=8000, ovlp=75, scale=FALSE)
points(res, pch=20)
legend(0.5, 3.6, "Fundamental frequency", pch=20, bty=0, cex=0.7)

Beep sound

Description

Generate a simple beep to be used as an alert, for instance at the end of a loop of when ending up a long script.

Usage

beep(d = 0.5, f = 8000, cf = 1000)

Arguments

Value

Nothing returned, a pure tone sound is played back. The default duration is 0.5 s and the default frequency is 1000 Hz

Note

The function uses listen of seewave which calls play of tuneR. You might need to set up your sound player with setWavPlayer of tuneR.

Author(s)

Jerome Sueur

Examples

## Not run: 
# default settings
beep()
# change the duration and the frequency
beep(d=1, cf=880)

## End(Not run)

Butterworth frequency filter

Description

This function is a Butterworth frequency filter that filters out a selected frequency section of of a time wave (low-pass, high-pass, low-stop, high-stop, bandpass or bandstop frequency filter).

Usage

bwfilter(wave, f, channel = 1, n = 1, from = NULL, to = NULL,
bandpass = TRUE, listen = FALSE, output = "matrix")

Arguments

Details

The order of the filter determines the value of the roll-off value, that is the dB decrease per octave of the transfer function. A filter of order n will have a transfer function with a roll-off value of - n*6 dB.

Value

A new wave is returned. The class of the returned object is set with the argument output.

Note

This function mainly uses the functions filter() and filtfilt() from the package signal

Author(s)

Jerome Sueur, functions filter() and filtfilt() from the package signal.

References

Stoddard, P. K. (1998). Application of filters in bioacoustics. In: Hopp, S. L., Owren, M. J. and Evans, C. S. (Eds), Animal acoustic communication. Springer, Berlin, Heidelberg,pp. 105-127.

See Also

ffilter, bwfilter, preemphasis, lfs, afilter

Examples

require(signal)
f <- 8000
a <- noisew(f=f, d=1)
## low-pass
# 1st order filter
res <- bwfilter(a, f=f, n=1, to=1500)
# 8th order filter
res <- bwfilter(a, f=f, n=8, to=1500)
## high-pass
res <- bwfilter(a, f=f, from=2500)
## band-pass
res <- bwfilter(a, f=f, from=1000, to=2000)
## band-stop
res <- bwfilter(a, f=f, from=1000, to=2000,bandpass=FALSE)

Continuous coherence function between two time waves

Description

This function returns a two-dimension coherence representation between two time waves. The function corresponds to a sliding coherence function along the two signals. This produces a 2-D density plot. An amplitude contour plot can be overlaid.

Usage

ccoh(wave1, wave2, f, channel = c(1,1), wl = 512, ovlp = 0, plot = TRUE,
grid = TRUE, scale = TRUE, cont = FALSE,
collevels = seq(0, 1, 0.01), palette = reverse.heat.colors,
contlevels = seq(0, 1, 0.01), colcont = "black",
colbg="white", colgrid = "black",
colaxis = "black", collab="black",
xlab = "Time (s)", ylab = "Frequency (kHz)",
scalelab = "Coherence",
main = NULL,
scalefontlab = 1, scalecexlab =0.75, axisX = TRUE, axisY = TRUE,
flim = NULL, flimd = NULL,
...)

Arguments

Details

Coherence is a frequency domain function computed to show the degree of a relationship between two signals. The value of the coherence function ranges between zero and one, where a value of zero indicates there is no causal relationship between the signals. A value of one indicates the existence of linear frequency response between the two signals. This can be used, for instance, to compare the input and output signals of a system.
Any colour palette can be used. In particular, it is possible to use other palettes coming with seewave: temp.colors, reverse.gray.colors.1, reverse.gray.colors.2, spectro.colors, reverse.terrain.colors, reverse.topo.colors, reverse.cm.colors corresponding to the reverse of terrain.colors, topo.colors, cm.colors.
Use locator to identify points.

Value

This function returns a list of three items:

Note

This function is based on spec.pgram, contour and filled.contour. See spectro for graphical changes.

Author(s)

Jerome Sueur sueur@mnhn.fr but this function is mainly based on spec.pgram by Martyn Plummer, Adrian Trapletti and B.D. Ripley

See Also

coh, spectro, spec.pgram.

Examples

wave1<-synth(d=1,f=4000,cf=500)
wave2<-synth(d=1,f=4000,cf=800)
ccoh(wave1,wave2,f=4000)

Cepstrum or real cepstrum

Description

This function returns the cepstrum of a time wave allowing fundamental frequency detection.

Usage

ceps(wave, f, channel = 1, phase = FALSE, wl = 512, at = NULL, from = NULL, to = NULL,
tidentify = FALSE, fidentify = FALSE, col = "black", cex = 1, plot = TRUE,
qlab = "Quefrency (bottom: s, up: Hz)", alab = "Amplitude",
qlim = NULL, alim = NULL, type = "l", ...)

Arguments

Details

The cepstrum of a time wave is the inverse Fourier transform of the logarithm of the Fourier transform. The cepstrum of a wave s is then calculated as follows:

C(s) = Re[FFT^{-1}(\log{(|FFT(s)|)]}

The independent variable of a cepstral graph is called the quefrency. The quefrency is a measure of time, though not in the sense of a signal in the time domain. A correspondence with the frequency domain is obtained by simply computing the reverse of the temporal x coordinate. For instance if a peak appears at 0.005 s, this reveals a frequency peak at 200 Hz (=1/0.005). This explain the two scales plotted when plot is TRUE.
If at, from or to are FALSE then ceps computes the cepstrum of the whole signal.
When using tidentify or tidentify, press ‘stop’ tools bar button to return values in the console.

Value

When plot is FALSE, ceps returns the cesptral profile as a two-column matrix, the first column corresponding to quefrency (x-axis) and the second corresponding to amplitude (y-axis).

Warning

The argument peaks is no more available (version > 1.5.6). See the function fpeaks for peak(s) detection.

Note

Cepstral analysis is mainly used in speech processing. This analysis allows to extract the fundamental frequency, see the examples.
This function is based on fft.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Oppenheim, A.V. and Schafer, R.W. 2004. From frequency to quefrency: a history of the cepstrum. Signal Processing Magazine IEEE, 21: 95-106.

See Also

cepstro, fund, autoc

Examples

data(sheep)
par(mfrow=c(2,1))
# phase not taken into account
ceps(sheep,f=8000,at=0.4,wl=1024)
# phase taken into account
ceps(sheep,f=8000,at=0.4,wl=1024, phase=TRUE)

2D-cepstrogram of a time wave

Description

This function returns a two-dimension cepstrographic representation of a time wave. The function corresponds to a short-term cepstral transform. An amplitude contour plot can be overlaid.

Usage

cepstro(wave, f, channel = 1, wl = 512, ovlp = 0, plot = TRUE, grid = TRUE,
scale = TRUE, cont = FALSE, collevels = seq(0, 1, 0.01),
palette = reverse.heat.colors, contlevels = seq(0, 1, 0.01),
colcont = "black", colbg="white", colgrid = "black",
colaxis = "black", collab = "black",
xlab = "Time (s)", ylab = "Quefrency (ms)",
scalelab = "Amplitude", main = NULL, scalefontlab = 1, scalecexlab = 0.75,
axisX = TRUE, axisY = TRUE, tlim = NULL, qlim = NULL, ...)

Arguments

Details

It is unfortunately not possible to turn the y-axis to a frequency scale.
See spectro for the use of the graphical arguments.

Value

This function returns a list of three items:

Note

This function is based on ceps.

Author(s)

Jerome Sueur sueur@mnhn.fr.

References

Oppenheim, A.V. and Schafer, R.W. 2004. From frequency to quefrency: a history of the cepstrum. Signal Processing Magazine IEEE, 21: 95-106.

See Also

ceps, fund, autoc

Examples

data(sheep)
sheepc <- cutw(sheep, f=8000, from = 0.19, to = 2.3)
cepstro(sheepc,f=8000)

Coherence between two time waves

Description

This function returns the frequency coherence between two time waves.

Usage

coh(wave1, wave2, f, channel=c(1,1), plot =TRUE, xlab = "Frequency (kHz)",
ylab = "Coherence", xlim = c(0,f/2000), type = "l",...)

Arguments

Details

Coherence is a frequency domain function computed to show the degree of a relationship between two signals. The value of the coherence function ranges between zero and one, where a value of zero indicates there is no causal relationship between the signals. A value of one indicates the existence of linear frequency response between the two signals. This can be used, for instance, to compare the input and output signals of a system.

Value

When plot is FALSE, this coh returns a two-column matrix, the first column being the frequency axis in kHz (x-axis) and the second column being the coherence (y-axis).

Note

This function is based on spec.pgram.

Author(s)

Jerome Sueur sueur@mnhn.fr but this function is based on spec.pgram by Martyn Plummer, Adrian Trapletti and B.D. Ripley.

See Also

ccoh, spectro, spec.pgram.

Examples

wave1<-synth(d=1,f=4000,cf=500)
wave2<-synth(d=1,f=4000,cf=800)
coh(wave1,wave2,f=4000)

Comb filter

Description

This function processes a feedforward comb filter and plots a spectrogram of the filtered wave asso- ciated with the frequency response of the filter.

Usage

combfilter(wave, f, channel = 1, alpha, K, units = c("samples", "seconds"),
plot = FALSE, output = "matrix", ...)

Arguments

Details

A comb filter consists in adding a delayed version of a signal to itself resulting in constructive and destructive interference. The feedforward version of a comb filter can be written following:

y(n) = x(n) + \alpha \times x(n - K)

where alpha is the scaling factor and K the delay length. The frequency response of the filter is obtained with:

H(f) = \sqrt{(1+\alpha^2)+2 \times \cos(\omega K)}

The frequency response is periodic. The depth of the cycles is controlled with alpha and the number of cycles with K.

Value

A new wave is returned. The class of the returned object is set with the argument output.

Note

Setting K to high values may generate unwanted results.
The feedback form of the combfilter is not implemented yet.

Author(s)

Jerome Sueur

See Also

combfilter, fir, squarefilter, drawfilter, ffilter, bwfilter

Examples

## Not run: 
f <- 44100
## chirp
s1 <- synth(f=f, cf=1, d=2, fm=c(0,0,f/2,0,0), out="Wave")
combfilter(s1, alpha=1, K=50, plot=TRUE)
## harmonic sound
s2 <- synth(f=f, d=2, cf=600, harmonics=rep(1, 35), output="Wave")
combfilter(s2, alpha=1, K=10, plot=TRUE)
## noise, units in seconds
s3 <- noisew(d=2, f=44100, out="Wave")
combfilter(s3, alpha=0.5, K=1e-4, units="seconds", plot=TRUE)

## End(Not run)

Convert sound pressure level in other units

Description

This function converts sound pressure level (in dB) in sound power (Watt), intensity (Watt/m2) and pressure (Pa). By default, these conversions are applied to air-borne sound.

Usage

convSPL(x, d = 1, Iref = 10^-12, pref = 2*10^-5)

Arguments

Value

convSPL returns a list containing three components:

Note

Iref and pref correspond to a 1 kHz sound in air.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Hartmann, W. M. 1998 Signals, sound and sensation. New York: Springer.

See Also

moredB, dBweight, attenuation

Examples

# conversion of two SPL measurements taken at 0.5 m from the source
convSPL(c(80,85),d=0.5) 

Cross-correlation between two time wave envelopes

Description

This function tests the similarity between two time wave envelopes by returning their maximal correlation and the time shift related to it.

Usage

corenv(wave1, wave2, f, channel=c(1,1), envt="hil", msmooth = NULL, ksmooth = NULL,
ssmooth = NULL, plot = TRUE, plotval = TRUE,
method = "spearman", col = "black", colval = "red",
cexval = 1, fontval = 1, xlab = "Time (s)",
ylab = "Coefficient of correlation (r)", type = "l", pb = FALSE, ...)

Arguments

Details

Successive correlations between the envelopes of wave1 and wave2 are computed when regularly sliding forward and backward wave2 along wave1.
The maximal correlation is obtained at a particular shift (time offset). This shift may be positive or negative.
The higher smooth is set up, the faster will be the computation but less precise the results will be.
The corresponding p value, obtained with cor.test, is plotted. Inverting wave1 and wave2 may give slight different results.

Value

If plot is FALSE, corenv returns a list containing four components:

Author(s)

Jerome Sueur

See Also

env,corspec,covspectro, cor,cor.test.

Examples

## Not run: 
data(orni)
# cross-correlation between two echemes of a cicada song
wave1<-cutw(orni,f=22050,from=0.3,to=0.4,plot=FALSE)
wave2<-cutw(orni,f=22050,from=0.58,to=0.68,plot=FALSE)
corenv(wave1,wave2,f=22050)

## End(Not run)

Cross-correlation between two frequency spectra

Description

This function tests the similarity between two frequency spectra by returning their maximal correlation and the frequency shift related to it.

Usage

corspec(spec1, spec2, f = NULL, mel = FALSE, plot = TRUE, plotval = TRUE,
method = "spearman", col = "black", colval = "red",
cexval = 1, fontval = 1, xlab = NULL,
ylab = "Coefficient of correlation (r)", type="l",...)

Arguments

Details

It is important not to have data in dB.
Successive correlations between spec1 and spec2 are computed when regularly shifting spec2 towards lower or higher frequencies.
The maximal correlation is obtained at a particular shift (frequency offset). This shift may be positive or negative.
The corresponding p value, obtained with cor.test, is plotted.
Inverting spec1 and spec2 may give slight different results, see examples.

Value

If plot is FALSE, corspec returns a list containing four components:

Author(s)

Jerome Sueur, improved by Laurent Lellouch

References

Hopp, S. L., Owren, M. J. and Evans, C. S. (Eds) 1998. Animal acoustic communication. Springer, Berlin, Heidelberg.

See Also

spec, meanspec, corspec, covspectro, cor, cor.test.

Examples

## Not run: data(tico)
## compare the two first notes spectra
a<-spec(tico,f=22050,wl=512,at=0.2,plot=FALSE)
c<-spec(tico,f=22050,wl=512,at=1.1,plot=FALSE)
op<-par(mfrow=c(2,1), mar=c(4.5,4,3,1))
spec(tico,f=22050,at=0.2,col="blue")
par(new=TRUE)
spec(tico,f=22050,at=1.1,col="green")
legend(x=8,y=0.5,c("Note A", "Note C"),lty=1,col=c("blue","green"),bty="o")
par(mar=c(5,4,2,1))
corspec(a,c, ylim=c(-0.25,0.8),xaxs="i",yaxs="i",las=1)
par(op)
## different correlation methods give different results...
op<-par(mfrow=c(3,1))
corspec(a,c,xaxs="i",las=1, ylim=c(-0.25,0.8))
title("spearmann correlation (by default)")
corspec(a,c,xaxs="i",las=1,ylim=c(0,1),method="pearson")
title("pearson correlation")
corspec(a,c,xaxs="i",las=1,ylim=c(-0.23,0.5),method="kendall")
title("kendall correlation")
par(op)
## inverting x and y does not give exactly similar results
op<-par(mfrow=c(2,1),mar=c(2,4,3,1))
corspec(a,c)
corspec(c,a)
par(op)
## mel scale
require(tuneR)
data(orni)
orni.mel <- melfcc(orni, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
orni.mel.mean <- apply(orni.mel$aspectrum, MARGIN=2, FUN=mean)
tico.mel <- melfcc(tico, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
tico.mel.mean <- apply(tico.mel$aspectrum, MARGIN=2, FUN=mean)
corspec(orni.mel.mean, tico.mel.mean, f=22050, mel=TRUE, plot=TRUE)

## End(Not run)

Covariance between two spectrograms

Description

This function tests the similarity between two spectrograms by returning their maximal covariance and the time shift related to it.

Usage

covspectro(wave1, wave2, f, channel = c(1,1), wl = 512, wn = "hanning", n,
plot = TRUE, plotval = TRUE,
method = "spearman", col = "black", colval = "red", cexval = 1,
fontval = 1, xlab = "Time (s)",
ylab = "Normalised covariance (cov)", type = "l", pb = FALSE, ...)

Arguments

Details

Successive covariances between the spectrogram of wave1 and the spectrogram of wave2 are computed when regularly sliding forward and backward wave2 along wave1.
The maximal covariance is obtained at a particular shift (time offset). This shift may be positive or negative.
n sets in how many steps wave2 will be slided along wave1. Time process can be then decreased by setting low n value.
Inverting wave1 and wave2 may give slight different results.

Value

If plot is FALSE, covspectro returns a list containing three components:

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Hopp, S. L., Owren, M. J. and Evans, C. S. (Eds) 1998. Animal acoustic communication. Springer, Berlin, Heidelberg.

See Also

corspec, corenv, spectro, cor,

Examples

# covariance between two notes of a birdsong
## Not run: 
data(tico)
note1<-cutw(tico, f=22050, from=0.5, to=0.9)
note2<-cutw(tico, f=22050, from=0.9, to=1.3)
covspectro(note1,note2,f=22050,n=37)

## End(Not run)

Crest factor and visualization

Description

This function returns the crest factor and localizes the different crest(s).

Usage

crest(wave, f, channel = 1, plot = FALSE, col = 2, cex = 3, symbol = "*", ...)

Arguments

Details

The crest factor of a time series s is calculated according to:

C = \frac{max(s)}{rms(s)}

with rms the root-mean-square (see rms).

Value

The function returns a list of three items

Note

There might be several crests (maxima) along the time wave but there is a single crest factor.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Hartmann, W. M. 1998 Signals, sound and sensation. New York: Springer.

See Also

oscillo, rms

Examples

data(tico)
crest(tico, f=22050)
# see the crest location and change the default graphical parameters
crest(tico, f=22050, plot=TRUE, sym="-")

Continuous spectral entropy

Description

This function computes the continuous spectral entropy (H) of a time wave.

Usage

csh(wave, f, channel = 1, wl = 512, wn = "hanning", ovlp = 0,
fftw = FALSE, threshold = NULL,
plot = TRUE, xlab = "Times (s)", ylab = "Spectral Entropy",
ylim = c(0, 1.1), type = "l", ...)

Arguments

Details

See sh for computing method.

Value

When plot is FALSE, csh returns a two-column matrix, the first column being time in seconds (x-axis) and the second column being the spectral entropy (y-axis) computed along time.
NA corresponds to pause sections in wave (see threshold).

Note

The spectral entropy of a noisy signal will tend towards 1 whereas the spectral entropy of a pure tone signal will tend towards 0.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Toh, A. M., Togneri, R. & Nordholm, S. 2005 Spectral entropy as speech features for speech recognition. Proceedings of PEECS, pp. 60-65.

See Also

sh, th

Examples

data(orni)
csh(orni,f=22050,wl=512,ovlp=50)
# using the threshold argument can lead to some edge effets
# here sh=1 at the end of echemes
csh(orni,f=22050,wl=512,ovlp=50,threshold=5)

Cut a frequency spectrum

Description

This function can be used to select (cut) a specific part of a frequency spectrum.

Usage

cutspec(spec, f = NULL, flim, mel = FALSE, norm = FALSE, PMF = FALSE)

Arguments

Value

A new spectrum is returned. The class of the returned object is the one of the input object (spec)

Note

The sampling frequency f is not necessary if spec has been obtained with either spec or meanspec.
This function can be used before calling analysis function like sh or sfm. See examples.

Author(s)

Jerome Sueur, improved by Laurent Lellouch

See Also

spec, meanspec

Examples

data(orni)
a <- meanspec(orni,f=22050,plot=FALSE)
b <- cutspec(a,flim=c(4,8))
## quick check with a plot
plot(b,type="l")
## effects on spectral entropy
sfm(a)
sfm(b)
## mel scale
require(tuneR)
mel <- melfcc(orni, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
melspec.mean <- apply(mel$aspectrum, MARGIN=2, FUN=mean)
c <- cutspec(melspec.mean, f=22050, flim=c(4000,8000), mel=TRUE)

Cut a section of a time wave

Description

This function selects and cuts a section of data describing a time wave. Original and cut sections can be plotted as oscillograms for comparison.

Usage

cutw(wave, f, channel=1, from = NULL, to = NULL, choose = FALSE,
plot = FALSE, marks = TRUE, output="matrix", ...)

Arguments

Details

If plot is TRUE returns a two-frame plot with both original and cut sections.

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Author(s)

Jerome Sueur

See Also

oscillo, addsilw,deletew, fadew,mutew,pastew,revw, zapsilw

Examples

# a 0.4 s section in a bird song
data(tico)
a<-cutw(tico,f=22050,from=0.5,to=0.9)
oscillo(a,22050)
# a direct way to see what has been cut
cutw(tico,f=22050,from=0.5,to=0.9,plot=TRUE)

dB colour scale for a spectrogram display

Description

This function displays a vertical or horizontal dB colour scale to be used with spectro plots.

Usage

dBscale(collevels, palette = spectro.colors, side = 4,
textlab = "Amplitude\n(dB)", cexlab = 0.75,
fontlab = 1, collab = "black", colaxis = "black",...)

Arguments

Note

This function, based on filled.contour by Ross Ihaka, is not supposed to be used by itself but as a legend of spectro.
Any colour palette can be used. In particular, it is possible to use other palettes coming with seewave: rev.gray.colors.1, rev.gray.colors.2, rev.heat.colors, rev.terrain.colors, rev.topo.colors, rev.cm.colors corresponding to the reverse of heat.colors, terrain.colors, topo.colors, cm.colors.

Author(s)

Jerome Sueur sueur@mnhn.fr and Caroline Simonis csimonis@mnhn.fr.

See Also

spectro.

Examples

data(pellucens)
# place the scale on the left and not on the right as spectro() does
def.par <- par(no.readonly = TRUE)
layout(matrix(c(1, 2), nc = 2), widths = c(1, 5))
par(mar=c(5,3,4,2))
dBscale(collevels=seq(-30,0,1),side=2)
par(mar=c(5,4,4,2))
spectro(pellucens, f=22050,wl=512,scale=FALSE)
par(def.par)
# place the scale on the top and not on the right as spectro() does
def.par <- par(no.readonly = TRUE)
layout(matrix(c(0,1,2,2), nc = 2, byrow=TRUE),widths=c(1,2),heights=(c(1,5.5)))
par(mar=c(0.5,3,4,2))
dBscale(collevels=seq(-30,0,1), textlab = "",side=3)
mtext("Amplitude (dB)",side=2,line = 1,at=0.6,cex=0.75)
par(mar=c(5,4,0.5,2))
spectro(pellucens, f=22050,wl=512,scale=FALSE)
par(def.par)   

dB weightings

Description

This function returns the four most common dB weightings.

Usage

dBweight(f, dBref = NULL)

Arguments

Details

By default, the function returns four weightings. When dBref is not NULL then the function returns the conversion from a dB reference level to four dB weighting levels.

Value

dBweight returns a list of five items corresponding to five dB weightings.

Note

The transfer equations used here come from Wipipedia but they were originally coming from the appendix of an international standard on the design performance of sound level meters IEC 651:1979 (Neil Glenister, pers. com.).

Author(s)

Jerome Sueur sueur@mnhn.fr, Zev Ross, and Andrey Anikin

References

https://en.wikipedia.org/wiki/A-weighting, https://en.wikipedia.org/wiki/ITU-R_468_noise_weighting

See Also

convSPL, moredB

Examples

# weight for a 50 Hz frequency
dBweight(f=50)
# A weight for the 1/3 Octave centre frequencies.
dBweight(f=c(20,25,31.5,40,50,63,80,100,125,160,200,250,
315,400,500,630,800,1000,1500,
1600,2000,2500,3150,4000,5000,
6300,8000,10000,12500,16000,20000))$A
# correction for a 50 Hz sound emitted at 100 dB
dBweight(f=50, dB=100)
# weighting curves plot
f <- seq(10,20000,by=10)
par(las=1)
plot(f, dBweight(f)$A, type="n", log="x",
xlim=c(10,10^5),ylim=c(-80,20),xlab="",ylab="",xaxt="n",yaxt="n")
abline(v=c(seq(10,100,by=10),seq(100,1000,by=100),
seq(1000,10000,by=1000),seq(10000,100000,by=10000),
c(100,1000,10000,100000)),col="lightgrey",lty=2)
abline(v=c(100,1000,10000,100000),col="grey")
abline(h=seq(-80, 20, 20),col="grey")
par(new=TRUE)
plot(f, dBweight(f)$A, type="l", log="x",
xlab="Frequency (Hz)", ylab="dB",lwd=2, col="blue", xlim=c(10,10^5),ylim=c(-80,20))
title(main="Acoustic weighting curves (10 Hz - 20 kHz)")
lines(x=f, y=dBweight(f)$B, col="green",lwd=2)
lines(x=f, y=dBweight(f)$C, col="red",lwd=2)
lines(x=f, y=dBweight(f)$D, col="black",lwd=2)
legend("bottomright",legend=c("dB(A)","dB(B)","dB(C)","dB(D)"),
lwd=2,col=c("blue","green","red","black"),bty="o",bg="white")

Delete a section of a time wave

Description

This function selects and delete a section of data describing a time wave. Original section and section after deletion can be plotted as oscillograms for comparison.

Usage

deletew(wave, f, channel = 1, from = NULL, to = NULL, choose = FALSE, plot = FALSE,
marks = TRUE, output = "matrix", ...)

Arguments

Details

If plot is TRUE returns a two-frame plot with both original and resulting sections.

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

oscillo, addsilw,cutw, fadew, mutew, pastew, revw, zapsilw

Examples

# deletion a 0.4 s section in a bird song
data(tico)
a<-deletew(tico,f=22050,from=0.5,to=0.9)
oscillo(a,22050)
# a direct way to see what has been cut
deletew(tico,f=22050,from=0.5,to=0.9,plot=TRUE)

Dominant frequency of a time wave

Description

This function gives the dominant frequency (i. e. the frequency of highest amplitude) of a time wave.

Usage

dfreq(wave, f, channel = 1, wl = 512, wn = "hanning", ovlp = 0, fftw=  FALSE, at =
NULL, tlim = NULL, threshold = NULL, bandpass = NULL, clip = NULL,
plot = TRUE, xlab = "Times (s)", ylab = "Frequency (kHz)",
ylim = c(0, f/2000), ...)

Arguments

Value

When plot is FALSE, dfreq returns a two-column matrix, the first column corresponding to time in seconds (x-axis) and the second column corresponding to to dominant frequency in kHz (y-axis).
NA corresponds to pause sections in wave (see threshold).

Note

This function is based on fft.

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

spec, meanspec,spectro.

Examples

data(tico)
f <- 22050
# default
dfreq(tico,f)
# using the amplitude threshold and changing the graphical output
dfreq(tico, f, ovlp=50,threshold=5, type="l", col=2)
# using 'at' argument for specific positions along the time axis
dfreq(tico, f, at=c(0.25, 0.75, 1.2, 1.6))
dfreq(tico, f, at=seq(0.5, 1.4, by=0.005), threshold=5)
# a specific number of measures on a single note
dfreq(tico, f, at=seq(0.5, 0.9, len=100), threshold=5, xlim=c(0.5,0.9))
# overlap on spectrogram
# and use of 'clip' argument to better track the dominant frequency
# in noisy conditions
op <- par()
ticon <- tico@left/max(tico@left) + noisew(d=length(tico@left)/f, f)
spectro(ticon, f)
res <- dfreq(ticon, f, clip=0.3, plot=FALSE)
points(res, col=2, pch =13)
par(op)

Difference between two cumulative frequency spectra

Description

This function compares two distributions (e.g. two frequency spectra) by computing the difference between two cumulative frequency spectra

Usage

diffcumspec(spec1, spec2, f = NULL, mel = FALSE, 
plot = FALSE, type = "l", lty = c(1, 2), col = c(2, 4, 8),
flab = NULL, alab = "Cumulated amplitude",
flim = NULL, alim = NULL,
title = TRUE, legend = TRUE, ...)

Arguments

Details

Both spectra are transformed into cumulative distribution functions (CDF).
Spectral difference is then computed according to:

D_{cf}(x, y) = \frac{\sum_{i=1}^{n}|X_{i} - Y_{i}|}{n}, with with X and Y the spectrum CDFs, and D \in [0,1].

Value

A numeric vector of length 1 returning the difference between the two spectra. No unit.

Note

This metric is sensitive not only to the spectral overlap between but also to the mean frequential distance between the different frequency peaks.

Author(s)

Laurent Lellouch, Jerome Sueur

References

Lellouch L, Pavoine S, Jiguet F, Glotin H, Sueur J (2014) Monitoring temporal change of bird communities with dissimilarity acoustic indices. Methods in Ecology and Evolution, in press.

See Also

kl.dist, ks.dist, simspec, diffspec, logspec.dist, itakura.dist

Examples

## Hz scale
data(tico)
data(orni)
orni.hz <- meanspec(orni, plot=FALSE)
tico.hz <- meanspec(tico, plot=FALSE)
diffcumspec(orni.hz, tico.hz, plot=TRUE)
## mel scale
require(tuneR)
orni.mel <- melfcc(orni, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
orni.mel.mean <- apply(orni.mel$aspectrum, MARGIN=2, FUN=mean)
tico.mel <- melfcc(tico, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
tico.mel.mean <- apply(tico.mel$aspectrum, MARGIN=2, FUN=mean)
diffcumspec(orni.mel.mean, tico.mel.mean, f=22050, mel=TRUE, plot=TRUE)

Difference between two amplitude envelopes

Description

This function estimates the surface difference between two amplitude envelopes.

Usage

diffenv(wave1, wave2, f, channel = c(1,1), envt = "hil", msmooth = NULL, ksmooth = NULL,
plot = FALSE, lty1 = 1, lty2 = 2, col1 = 2, col2 = 4, cold = 8,
xlab = "Time (s)", ylab = "Amplitude", ylim = NULL, legend = TRUE, ...)

Arguments

Details

D is a Manhattan distance (l1 norm).
Envelopes of both waves are first transformed as probability mass functions (PMF).
Envelope difference is then computed according to:

D = \frac{\sum{|env1-env2|}}{2}, with D \in [0,1].

Value

The difference is returned. This value is without unit. When plot is TRUE, both envelopes and their difference surface are plotted on the same graph.

Note

This method can be used as a relative distance estimation between different envelopes.

Author(s)

Jerome Sueur sueur@mnhn.fr.

References

Sueur, J., Pavoine, S., Hamerlynck, O. & Duvail, S. (2008) - Rapid acoustic survey for biodiversity appraisal. PLoS ONE, 3(12): e4065.

See Also

env, corenv, diffspec, diffwave

Examples

data(tico) ; tico <- tico@left
data(orni) ; orni <- orni@left
# selection in tico of two waves with similar duration
tico2<-tico[1:length(orni)]
diffenv(tico2,orni,f=22050,plot=TRUE)
# smoothing the envelope gives a better graph but slightly changes the result
diffenv(tico2,orni,f=22050,msmooth=c(20,0),plot=TRUE)

Difference between two frequency spectra

Description

This function estimates the surface difference between two frequency spectra.

Usage

diffspec(spec1, spec2, f = NULL, mel = FALSE,
plot = FALSE, type="l", 
lty=c(1, 2), col =c(2, 4, 8),
flab = NULL, alab = "Amplitude",
flim = NULL, alim = NULL, title = TRUE, legend = TRUE, ...)

Arguments

Details

D is a Manhattan distance (l1 norm).
Both spectra are first transformed as probability mass functions (PMF).
Spectral difference is then computed according to:

D = \frac{\sum{|spec1-spec2|}}{2}, with D \in [0,1].

, with 0 < D < 1.

Value

The difference is returned. This value is without unit. When plot is TRUE, both spectra and their difference surface are plotted on the same graph.

Note

This method can be used as a relative distance estimation between different spectra.
The dB value obtained can be very different from the one visually estimated when looking at the graph (plot=TRUE).

Author(s)

Jerome Sueur, Sandrine Pavoine and Laurent Lellouch

References

Sueur, J., Pavoine, S., Hamerlynck, O. and Duvail, S. (2008). Rapid acoustic survey for biodiversity appraisal. PLoS One, 3(12): e4065.

See Also

spec, meanspec, corspec, simspec, diffcumspec, diffenv, kl.dist, ks.dist, logspec.dist, itakura.dist

Examples

a <- noisew(f=8000,d=1)
b <- synth(f=8000,d=1,cf=2000)
c <- synth(f=8000,d=1,cf=1000)
d <- noisew(f=8000,d=1)
speca <- spec(a,f=8000,wl=512,at=0.5,plot=FALSE)
specb <- spec(b,f=8000,wl=512,at=0.5,plot=FALSE)
specc <- spec(c,f=8000,wl=512,at=0.5,plot=FALSE)
specd <- spec(d,f=8000,wl=512,at=0.5,plot=FALSE)
diffspec(speca,speca,f=8000)
#[1] 0 => similar spectra of course !
diffspec(speca,specb)
diffspec(speca,specc,plot=TRUE)
diffspec(specb,specc,plot=TRUE)
diffspec(speca,specd,plot=TRUE)
## mel scale
require(tuneR)
data(orni)
data(tico)
orni.mel <- melfcc(orni, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
orni.mel.mean <- apply(orni.mel$aspectrum, MARGIN=2, FUN=mean)
tico.mel <- melfcc(tico, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
tico.mel.mean <- apply(tico.mel$aspectrum, MARGIN=2, FUN=mean)
diffspec(orni.mel.mean, tico.mel.mean, f=22050, mel=TRUE, plot=TRUE)

Difference between two time waves

Description

This function estimates the difference between two waves by computing the product between envelope surface difference and frequency surface difference.

Usage

diffwave(wave1, wave2, f, channel = c(1,1), wl = 512, envt = "hil",
msmooth = NULL, ksmooth = NULL)

Arguments

Details

D is a Manhattan distance (l1 norm).
This function computes the product between the values obtained with diffspec and diffenv functions.
This then gives a global (time and frequency) estimation of dissimilarity.
The frequency mean spectrum and the amplitude envelope needed for computing respectively diffspec and diffenv are automatically generated. They can be controlled through wl, msmooth and ksmooth arguments respectively.
See examples below and examples in diffspec and diffenv for implications on the results.

Value

A single value varying between 0 and 1 is returned. The value has no unit.

Note

This method can be used as a relative distance estimation between different waves.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Sueur, J., Pavoine, S., Hamerlynck, O. & Duvail, S. (2008) - Rapid acoustic survey for biodiversity appraisal. PLoS ONE, 3(12): e4065.

See Also

diffspec, diffenv

Examples

data(tico) ; tico <- tico@left
data(orni) ; orni <- orni@left
# selection in tico to have two waves of similar duration (length)
tico <- tico[1:length(orni)]
diffwave(tico,orni,f=22050)
# changing the frequency parameter (wl)
diffwave(tico,orni,f=22050,wl=1024)
# changing the temporal parameter (msmooth)
diffwave(tico,orni,f=22050,msmooth=c(20,0))

Time series discretisation

Description

This function transforms a numeric (time) series into a sequence of symbols

Usage

discrets(x, symb = 5, collapse = TRUE, plateau=1)

Arguments

Details

The function partitions the numeric (time) series into a sequence of finite number of symbols. These symbols result of the comparaison of each series value with its temporal neighbours.
They are two discretisations available:
when symb is set to 3, each value will be replaced by either:
- I if the series is Increasing,
- D if the series is Decreasing,
- F if the series remains Flat,
when symb is set to 5, each value will be replaced by either:
- I if the series is Increasing,
- D if the series is Decreasing,
- F if the series remains Flat,
- P if the series shows a Peak,
- T if the series shows a Trough.

The argument plateau can be used to control the way a plateau is encoded. A plateau is an elevated flat region that can be either considered a 'flat peak' encoded as PF...FP (plateau = 1) or as an increase, a flat region and a decrease encoded as IF...FD (plateau = 1. The default value (plateau = 1) refers to Cazelles et al. (2004).

Value

A character string of length 1 if collapse is TRUE. Otherwise, a character string of length n-2 if symbol=5 (the first and last values cannot be replaced with a symbol) or n-1 if symbol=3 (the first value cannot be replaced with a symbol.)

Author(s)

Jerome Sueur, improved by Laurent Lellouch

References

Cazelles, B. 2004 Symbolic dynamics for identifying similarity between rhythms of ecological time series. Ecology Letters, 7: 755-763.

See Also

symba

Examples

# a random variable
discrets(rnorm(30))
discrets(rnorm(30),symb=3)
# a frequency spectrum
data(tico)
spec1<-spec(tico,f=22050,at=0.2,plot=FALSE)
discrets(spec1[,2])

Draw the amplitude envelope of a time wave

Description

This function lets the user modifying the amplitude envelope of a time wave by drawing it with the graphics device

Usage

drawenv(wave, f, channel = 1, n = 20, plot = FALSE, listen = FALSE, output = "matrix")

Arguments

Details

The function first plots an oscillogram view of wave.
The user has then to choose points on the positive side of the y-axis (amplitude). The junction of these points will draw a new amplitude envelope.
The order of points along the x-axis (time) is not important but points cannot be cancelled. When this process is finished the new time wave is returned in the console or as an oscillogram in a second graphics device if plot is TRUE.
The function uses locator.

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

setenv, env, synth

Examples

## Not run: 
a<-synth(d=1,f=22050,cf=1000)
# drawenv(a,f=22050,plot=TRUE)
# choose points on the oscillogram view to draw a new enveloppe
# stop (ESC on Windows; right mouse button on Linux)
# check the result on the second graphics device opened thanks to plot=TRUE

## End(Not run)

Draw the amplitude profile of a frequency filter

Description

This function lets the user drawing the amplitude profile of a frequency filter.

Usage

drawfilter(f, n = 256, continuous = TRUE, discrete = TRUE)

Arguments

Details

If the same frequency of a discrete filter is selected twice then the sum of the amplitudes of the two selections is used. If both arguments continuous and discrete are set to TRUE and if frequencies selected overlap between the two filters then only the frequencies of the discrete filter are considered.

Value

The function returns a two-column matrix, the first column is the frequency in kHz and the second column is the amplitude of the filter.

Note

This function can be used to prepare bandpass or bandstop custom filters to be used with fir and ffilter. See examples.

Author(s)

Laurent Lellouch

See Also

fir, squarefilter, combfilter, ffilter, drawenv

Examples

## Not run: 
f <- 8000
a <- noisew(f=f, d=1)
## bandpass continuous and discrete
cont.disc <- drawfilter(f=f/2)
a.cont.disc <- fir(a, f=f, custom=cont.disc)
spectro(a.cont.disc, f=f)
## bandpass continuous only
cont <- drawfilter(f=f/2, discrete=FALSE)
a.cont <- fir(a, f=f, custom=cont)
spectro(a.cont, f=f)
## bandstop continuous only
cont.stop <- drawfilter(f=f/2, discrete=FALSE)
a.cont.stop <- fir(a, f=f, custom=cont.stop, bandpass=FALSE)
spectro(a.cont.stop, f=f)
## bandpass discrete only
disc <- drawfilter(f=f/2, continuous=FALSE)
a.disc <- fir(a, f=f, custom=disc, bandpass=FALSE)
spectro(a.disc, f=f)

## End(Not run)

Duration of a time wave

Description

Returns the duration (in second) of a time wave

Usage

duration(wave, f, channel=1)

Arguments

Value

A numeric vector of length 1 returning the duration in second.

Author(s)

Jerome Sueur

Examples

data(tico)
duration(tico)

Dynamic oscillogram

Description

This graphical function displays a time wave as an windowed oscillogram.

Usage

dynoscillo(wave, f, channel = 1, wd = NULL, wl = NULL, wnb = NULL, title = TRUE, ...)

Arguments

Details

The arguments wd, wl and wn have to be used isolated, not in conjunction. They basically do the same, ie they set the duration of the zooming window that is slided along the signal. For instance, for a 5 seconds sound with a sampling rate (f) at 44.1 kHz, wl = 4096 is equivalent to wd = 4096 / 44100 = 0.093 s and equivalent to wnb = 5*4096 / 44100 = 53.

Note

This function requires the package rpanel.

Author(s)

Jerome Sueur

See Also

oscillo, oscilloST, dynspec.

Examples

## Not run: 
require(rpanel)
data(tico)
dynoscillo(tico, wn=4)
## End(Not run)

Dynamic sliding spectrum

Description

This function plots dynamically a sliding spectrum along a time wave. This basically corresponds to a short-term Fourier transform.

Usage

dynspec(wave, f, channel = 1, wl = 512, wn = "hanning", zp = 0,
ovlp = 0, fftw = FALSE, norm = FALSE, dB = NULL,  dBref = NULL, plot = TRUE,
title = TRUE, osc = FALSE, 
tlab = "Time (s)", flab = "Frequency (kHz)",
alab = "Amplitude", alim = NULL, flim = c(0, f/2000),
type = "l", from = NULL, to = NULL, envt = NULL,
msmooth = NULL, ksmooth = NULL, colspec = "black",
coltitle = "black", colbg = "white", colline = "black",
colaxis = "black", collab = "black", cexlab = 1,
fontlab = 1, colwave = "black",
coly0 = "lightgrey", colcursor = "red", bty = "l")

Arguments

Details

Use the slider panel to move along the time wave.
Use the argument norm if you wish to have each spectrum normalised, i.e. with values between 0 and 1 or maximised to 0 dB when dB is TRUE.
The function requires the package rpanel that is based on the package tcltk.

Value

This function returns a list of three items:

Note

This function is very similar to a spectrogram. See the Details of spectro for some information regarding the short term Fourier transform.

Author(s)

Jerome Sueur and Caroline Simonis

See Also

spectro, spectro3D, wf, spec, dynspectro, fft, oscillo, env.

Examples

## Not run: 
data(sheep)
require(rpanel)
dynspec(sheep,f=8000,wl=1024,ovlp=50,osc=TRUE)

## End(Not run)

Dynamic sliding spectrogramn

Description

This function plots dynamically a sliding spectrogram along a time wave.

Usage

dynspectro(wave, f, channel = 1, slidframe = 10,
wl = 512, wn = "hanning", zp = 0, ovlp = 75,
fftw = FALSE, dB = TRUE, plot = TRUE,
title = TRUE, osc = FALSE,
tlab = "Time (s)", flab = "Frequency (kHz)", alab = "Amplitude",
from = NULL, to = NULL,
collevels = NULL, palette = spectro.colors,
envt = NULL, msmooth = NULL, ksmooth = NULL,
coltitle = "black", colbg = "white", colline = "black",
colaxis = "black", collab = "black", cexlab = 1,
fontlab = 1, colwave = "black",
coly0 = "lightgrey", colcursor = "red", bty = "l")

Arguments

Details

Use the slider panel to move along the time wave.
The function requires the package rpanel that is based on the package tcltk.
The function is mainly written for inspecting long sounds.
The function is based on image for fast display when spectro is based on filled.contour. Displaying the amplitude envelope with the argument envt can slow down significantly the display.

Value

This function returns a list of three items:

Note

This function is very similar to a spectrogram. See the Details of spectro for some information regarding the short term Fourier transform.

Author(s)

David Pinaud and Jerome Sueur

See Also

spectro, spectro3D, wf, spec, dynspec, fft, oscillo, env.

Examples

## Not run: 
data(sheep)
require(rpanel)
dynspectro(sheep, ovlp=95, osc=TRUE)

## End(Not run)

Echo generator

Description

This function generate echoes of a time wave.

Usage

echo(wave, f, channel = 1, amp, delay, plot = FALSE,
listen = FALSE, output = "matrix", ...)

Arguments

Details

amp and delay should strictly have the same length corresponding to the number of desired echoes.

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Note

This function is based on a convolution (convolve) between the input wave and a pulse echo filter.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Stoddard, P. K. (1998). Application of filters in bioacoustics. In: Hopp, S. L., Owren, M. J. and Evans, C. S. (Eds), Animal acoustic communication. Springer, Berlin, Heidelberg,pp. 105-127.

See Also

synth

Examples

# generation of the input wave
a <- synth(f=11025,d=1,cf=2000,shape="tria",am=c(50,10),fm=c(1000,10,1000,0,0))
# generation of three echoes
# with respectively a relative amplitude of 0.8, 0.4, and 0.2
# and with a delay of 1s, 2s, and 3s  from the beginning of the input wave
aecho <- echo(a,f=11025,amp=c(0.8,0.4,0.2),delay=c(1,2,3))
# another echo with time delays overlapping with the input wave
aecho <- echo(a,f=11025,amp=c(0.4,0.2,0.4),delay=c(0.6,0.8,1.5))

Amplitude envelope of a time wave

Description

This function returns the absolute or Hilbert amplitude envelope of a time wave.

Usage

env(wave, f, channel = 1, envt = "hil", 
msmooth = NULL, ksmooth = NULL, ssmooth = NULL,
asmooth = NULL,
fftw = FALSE, norm = FALSE,
plot = TRUE, k = 1, j = 1, ...)

Arguments

Details

When envt is set as "abs", the amplitude envelope returned is the absolute value of wave.
When envt is set as "hil", the amplitude envelope returned is the modulus (Mod) of the analytical signal of wave obtained through the Hilbert transform (hilbert).

Value

Data are returned as one-column matrix when plot is FALSE.

Note

Be aware that smoothing with either msmooth or ksmooth changes the original number of points describing wave.

Author(s)

Jerome Sueur. Implementation of 'fftw' argument by Jean Marchal and Francois Fabianek. Implementation of 'asmooth' by Thibaut Marin-Cudraz.

See Also

oscillo,hilbert

Examples

data(tico)
# Hilbert amplitude envelope
env(tico)
# absolute amplitude envelope
env(tico, envt="abs")
# smoothing with a 10 points and 50% overlaping mean sliding window
env(tico, msmooth=c(10,50))
# smoothing kernel
env(tico, ksmooth=kernel("daniell",10))
# sum smooth
env(tico, ssmooth=50)
# autocorrelation smooth
env(tico, asmooth=50)
# overplot of oscillographic and envelope representations
oscillo(tico)
par(new=TRUE)
env(tico, colwave=2)

Export sound data

Description

Export sound data as a text file that can be read by a sound player like 'Goldwave'

Usage

export(wave, f = NULL, channel = 1, filename = NULL, header=TRUE, ...)

Arguments

Details

Creates a new text file with a header describing the main features of the sound (wave). For instance, for a 2 s sound with a sampling frequency of 8000 Hz, the header will be: [ASCII 8000Hz, Channels: 1, Samples: 160000, Flags: 0]. This type of file can be read by sound players like Goldwave (http://www.goldwave.com/).

Author(s)

Jerome Sueur sueur@mnhn.fr

Examples

a<-synth(f=8000,d=2,cf=2000,plot=FALSE)
export(a,f=8000)
unlink("a.txt")

Fade in and fade out of a time wave

Description

This function applies a “fade in” and/or a “fade out” to a time wave following a linear, exponential or cosinus-like shape.

Usage

fadew(wave, f, channel = 1, din = 0, dout = 0, shape = "linear", plot = FALSE,
listen = FALSE, output = "matrix", ...)

Arguments

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

oscillo, addsilw, cutw, deletew,mutew, pastew, revw, zapsilw

Examples

a<-noisew(d=5,f=4000)
op<-par(mfrow=c(3,1))
fadew(a,f=4000,din=1,dout=2,plot=TRUE,title="Linear",cexlab=0.8)
fadew(a,f=4000,din=1,dout=2,shape="exp",plot=TRUE,title="Exponential shape",
    colwave="blue",coltitle="blue",cexlab=0.8)
fadew(a,f=4000,din=1,dout=2,shape="cos",plot=TRUE,title="Cosinus-like shape",
    colwave="red",coltitle="red",cexlab=0.8)
par(op)

Frequency bands plot

Description

This graphical function returns a frequency spectrum as a bar plot.

Usage

fbands(spec, f = NULL, bands = 10, width = FALSE, mel = FALSE, plot = TRUE,
xlab = NULL, ylab = "Relative amplitude", ...)

Arguments

Details

The function proceeds as follows

• divides the spectrum in bands. The limits of the bands are set with the argument bands. There are two options:

  • you set a number of bands with equal size by giving a single value to bands. For instance, setting bands to a value of 10 will slice the spectrum in 10 equal parts and return 10 local peaks.

  • you set the limits of the bands. This is achieve by giving a numeric vector to bands. The limits can follow a regular or irregular series. For instance attributing the vector c(0,2,4,8) will generate the following bands [0,2[, [2,4[, [4,8] kHz. Be aware that the last value should not exceed half the sampling frequency used to obtain the spectrum spec.

• uses the function barplot.

Value

A two-column matrix, the first column corresponding to the frequency values (x-axis, mean of the bars limits) and the second column corresponding to height values (y-axis) of the bars.

Note

The value below bars is the mean between the corresponding frequency limits.

Author(s)

Jerome Sueur, improved by Laurent Lellouch

See Also

meanspec, spec, barplot.

Examples

data(sheep)
spec <- meanspec(sheep, f=8000, plot=FALSE)
# default plot
fbands(spec)
# setting a specific number of bands
fbands(spec, bands=6)
#setting specific regular bands limits
fbands(spec, bands=seq(0,4,by=0.25))
# some plot tuning
op <- par(las=1)
fbands(spec, bands=seq(0,4,by=0.1),
       horiz=TRUE, col=heat.colors(41),
       xlab="", ylab="",
       cex.axis=0.75, cex.names = 0.75,
       axes=FALSE)
par(op)
# showing or not the width of the bands
oct <- octaves(440,3)/1000
op <- par(mfrow=c(2,1))
fbands(spec, bands=oct, col="blue")
fbands(spec, bands=oct, width = TRUE, col="red")
par(op)
# kind of horizontal zoom
op <- par(mfrow=c(2,1))
fbands(spec, bands=seq(0,4,by=0.2), col=c(rep(1,10),
   rep("orange",5),rep(1,5)), main="all frequency range")
fbands(spec, bands=seq(2,3,by=0.2),
   col="orange", main="a subset or zoom in")
par(op)
# kind of dynamic frequency bands
specs <- dynspec(sheep, f=8000, plot= FALSE)$amp
out <- apply(specs, f=8000, MARGIN=2,
     FUN = fbands, bands = seq(0,4,by=0.2),
     col = 1, ylim=c(0,max(specs)))
# mel scale
require(tuneR)
mel <- melfcc(sheep, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
melspec.mean <- apply(mel$aspectrum, MARGIN=2, FUN=mean)
melspec.mean <- melspec.mean/max(melspec.mean) # [0,1] scaling 
fbands(melspec.mean, f=8000, bands=8)

Doppler effect

Description

This function computes the altered frequency of a moving source due to the Doppler effect.

Usage

fdoppler(f, c = 340, vs, vo = 0, movs = "toward", movo = "toward")

Arguments

Details

The altered frequency f' is computed according to:

f{'} = f\times{\frac{c \pm v_{o}}{c \pm v_{s}}}

with f = original frequency produced by the source (in Hz or kHz),
vs = speed of the source,
vo = speed of the observer.

Value

The altered frequency is returned in a vector.

Note

You can use wasp to have exact values of c. See examples.

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

wasp

Examples

# a 400 Hz source moving toward or away from the observer at 85 m/s
fdoppler(f=400,vs=85)
# [1] 533.3333
fdoppler(f=400,vs=85,movs="away")
# [1] 320
# use wasp() if you wish to have exact sound speed at a specific temperature
fdoppler(f=wasp(f=400,t=25)$c, vs=85)
# [1] 461.8667
# Doppler effect at different source speeds
f<-seq(1,10,by=1); lf<-length(f)
v<-seq(10,300,by=20); lv<-length(v)
res<-matrix(numeric(lf*lv),ncol=lv)
for(i in 1:lv) res[,i]<-fdoppler(f=f,vs=v[i])
op<-par(bg="lightgrey")
matplot(x=f,y=res,type="l",lty=1,las=1,col= spectro.colors(lv),
xlab="Source frequency (kHz)", ylab="Altered frequency (kHz)")
legend("topleft",legend=paste(as.character(v),"m/s"),
lty=1,col= spectro.colors(lv))
title(main="Doppler effect at different source speeds")
par(op)

Frequency filter

Description

This function filters out a selected frequency section of a time wave (low-pass, high-pass, low-stop, high-stop, bandpass or bandstop frequency filter).

Usage

ffilter(wave, f, channel = 1, from = NULL, to = NULL, bandpass = TRUE,
custom = NULL, wl = 1024, ovlp = 75, wn = "hanning", fftw = FALSE,
rescale=FALSE, listen=FALSE, output="matrix")

Arguments

Details

A short-term Fourier transform is first applied to the signal (see spectro), then the frequency filter is applied and the new signal is eventually generated using the reverse of the Fourier Transform (istft).
There is therefore neither temporal modifications nor amplitude modifications.

Value

If plot is FALSE, a new wave is returned. The class of the returned object is set with the argument output.

Author(s)

Jerome Sueur

See Also

afilter,lfs,fir, preemphasis, combfilter, bwfilter

Examples

a<-noisew(f=8000,d=1)
# low-pass
b<-ffilter(a,f=8000,to=1500)
spectro(b,f=8000,wl=512)
# high-pass
c<-ffilter(a,f=8000,from=2500)
spectro(c,f=8000,wl=512)
# band-pass
d<-ffilter(a,f=8000,from=1000,to=2000)
spectro(d,f=8000,wl=512)
# band-stop
e<-ffilter(a,f=8000,from=1500,to=2500,bandpass=FALSE)
spectro(e,f=8000,wl=512)
# custom
myfilter1<-rep(c(rep(0,64),rep(1,64)),4)
g<-ffilter(a,f=8000,custom=myfilter1)
spectro(g,f=8000)

Near field and far field limits

Description

This function helps in knowing whether you are working in the near or far field.

Usage

field(f, d)

Arguments

Details

Areas very close to the sound source are in the near-field where the contribution of particle velocity to sound energy is greater thant that of sound pressure and where these components are not in phase. Sound propagation properties are also different near or far from the source. It is therefore important to know where the microphone was from the source.
To know this, the product k*d is computed according to:

k\times{d} = \frac{f}{c}\times{d}

with d = distance from the source (m), f = frequency (Hz) and c = sound celerity (m/s).
If k*d is greatly inferior 1 then the microphone is in the near field.
The decision help returned by the function follows the rule:
far field:

k\times{d} > 1

between near and far field limits:

0.1 \leq k\times{d} \leq 1

near field:

k\times{d} < 0.1

.

Value

A list of two values is returned:

Note

This function works for air-borne sound only.

Author(s)

Jerome Sueur sueur@mnhn.fr

Examples

# 1 kHz near field at 1 cm from the source
field(f=1000,d=0.01)
# playing with distance from source and sound frequency
op<-par(bg="lightgrey")
D<-seq(0.01,0.5,by=0.01); nD<-length(D)
F<-seq(100,1000,by=25); nF<-length(F)
a<-matrix(numeric(nD*nF),nrow=nD)
for(i in 1:nF) a[,i]<-field(f=F[i],d=D)$kd
matplot(x=D,y=a,type="l",lty=1,col= spectro.colors(nF),
  xlab="Distance from the source (m)", ylab="k*d")
title("Variation of the product k*d with distance and frequency")
text(x=c(0.4,0.15),y=c(0.02,1), c("Near Field","Far Field"),font=2)
legend(x=0.05,y=1.4,c("100 Hz","1000 Hz"),lty=1,
  col=c(spectro.colors(nF)[1],spectro.colors(nF)[nF]),bg="grey")
abline(h=0.1)
par(op)

Finite Impulse Response filter

Description

This function is a FIR filter that filters out a selected frequency section of a time wave (low-pass, high-pass, low-stop, high-stop, bandpass or bandstop frequency filter).

Usage

fir(wave, f, channel = 1, from = NULL, to = NULL, bandpass = TRUE, custom = NULL,
wl = 512, wn = "hanning", rescale=FALSE, listen = FALSE, output = "matrix")

Arguments

Details

This function is based on the reverse of the Fourier Transform (fft) and on a convolution (convolve) between the wave to be filtered and the impulse filter.

Value

A new wave is returned. The class of the returned object is set with the argument output.

Author(s)

Jerome Sueur

References

Stoddard, P. K. (1998). Application of filters in bioacoustics. In: Hopp, S. L., Owren, M. J. and Evans, C. S. (Eds), Animal acoustic communication. Springer, Berlin, Heidelberg,pp. 105-127.

See Also

ffilter, bwfilter, preemphasis, lfs, afilter

Examples

a<-noisew(f=8000,d=1)
# low-pass
b<-fir(a,f=8000,to=1500)
spectro(b,f=8000)
# high-pass
c<-fir(a,f=8000,from=2500)
spectro(c,f=8000)
# band-pass
d<-fir(a,f=8000,from=1000,to=2000)
spectro(d,f=8000)
# band-stop
e<-fir(a,f=8000,from=1500,to=2500,bandpass=FALSE)
spectro(e,f=8000)
# custom filter manually generated
myfilter1<-rep(c(rep(0,32),rep(1,32)),4)
g<-fir(a,f=8000,custom=myfilter1)
spectro(g,f=8000)
# custom filter generated using spec()
data(tico)
myfilter2<-spec(tico,f=22050,at=0.7,wl=512,plot=FALSE)
b<-noisew(d=1,f=22050)
h<-fir(b,f=22050,custom=myfilter2)
spectro(h,f=22050)

Frequency modulation analysis

Description

This function computes the Fourier analysis of the instantaneous frequency of a time wave. This allows to detect periodicity in frequency modulation.

Usage

fma(wave, f, channel = 1, threshold = NULL, plot = TRUE, ...)

Arguments

Details

This function is based on ifreq and spec.
The instantaneous frequency of wave is first computed and the spectrum of this frequency modulation is then processed. All env and spec arguments can be set up.

Value

If plot is FALSE, fma returns a numeric vector corresponding to the computed spectrum. If peaks is not NULL, fma returns a list with two elements:

Author(s)

Jerome Sueur sueur@mnhn.fr

See Also

ifreq, hilbert, spec, ama

Examples

# a sound with a 1 kHz sinusoid FM
a<-synth(d=1, f=8000, cf=1500, fm=c(1000,1000,0,0,0), output="Wave")
fma(a)

Frequency peak detection

Description

This function searches for peaks of a frequency spectrum.

Usage

fpeaks(spec, f = NULL,
nmax = NULL, amp = NULL, freq = NULL, threshold = NULL,
mel =FALSE,
plot = TRUE, title = TRUE,
xlab = NULL, ylab = "Amplitude",
labels = TRUE, digits = 2,
legend = TRUE, collab = "red", ...)

Arguments

Details

Here are some details regarding the different selection parameters:

• nmax: this parameter is to be used if you wish to get a specific number of peaks. The peaks selected are those with the highest slopes. It then does not work in conjunction with the other parameters.

• freq: this parameter allows to remove from the selection successive peaks with a small frequency difference. Imagine you have two successive peaks at 1200 Hz and 1210 Hz and at 0.5 and 0.25 in amplitude. If you set freq to 50 Hz, then only the first peak will be kept.

• amp: this parameter allows to remove from the selection peaks with low slopes. You can make the selection on both slopes or on a single one. Imagine you have an asymetric peak with a 0.01 left slope and a 0.02 right slope. The peak will be discarded for the following settings: both values higher than 0.02 (e.g. amp = c(0.03,0.04)), the first value higher than 0.01 (e.g. amp = c(0.02,0.001)), the second value higher than 0.02 (e.g. amp = c(0.001,0.03)). If you do not want apply the selection on one of the slope use 0. For instance, a selection on the left slope only will be achieved with: amp = c(0.02,0).

• threshold: this parameter can be used to do a rough selection on the spectrum. Peaks with an amplitude value (not a slope) lower than this threshold will be automatically discarded. This can be useful when you want to remove peaks of a low-amplitude background noise.

Value

A two-column matrix, the first column corresponding to the frequency values (x-axis) and the second column corresponding to the amplitude values (y-axis) of the peaks.

Note

You can also use fpeaks with other kind of spectrum, for instance a cepstral spectrum. See examples.

Author(s)

Jerome Sueur and Amandine Gasc

See Also

localpeaks, meanspec, spec

Examples

data(tico)
spec <- meanspec(tico, f=22050, plot=FALSE)
specdB <- meanspec(tico, f=22050, dB="max0", plot=FALSE)
# all peaks
fpeaks(spec)
# 10 highest peaks
fpeaks(spec, nmax=10)
# highest peak (ie dominant frequency)
fpeaks(spec, nmax=1)
# peaks that are separated by more than 500 Hz
fpeaks(spec, freq=500)
# peaks with a left slope higher than 0.1
fpeaks(spec, amp=c(0.1,0))
# peaks with a right slope higher than 0.1
fpeaks(spec, amp=c(0,0.1))
# peaks with left and right slopes higher than 0.1
fpeaks(spec, amp=c(0.1,0.1))
# peaks above a 0.5 threshold
fpeaks(spec, threshold=0.5)
# peaks of a dB spectrum with peaks showing slopes higher than 3 dB
fpeaks(specdB, amp=c(3,3))
# comparing different parameter settings
meanspec(tico, f=22050)
col <- c("#ff000090","#0000ff75","#00ff00")
cex <- c(2,1.25,1.5)
pch <- c(19,17,4)
title(main="Peak detection \n (spectrum with values between 0 and 1)")
res1 <- fpeaks(spec, plot = FALSE)
res2 <- fpeaks(spec, amp=c(0.02,0.02), plot =FALSE)
res3 <- fpeaks(spec, amp=c(0.02,0.02), freq=200, plot = FALSE)
points(res1, pch=pch[1], col=col[1], cex=cex[1])
points(res2, pch=pch[2], col=col[2], cex=cex[2])
points(res3, pch=pch[3], col=col[3], cex=cex[3])
legend("topright", legend=c("all peaks","amp", "amp & freq"), pch=pch,
pt.cex=cex, col=col, bty="n")
# example with a cepstral spectrum
data(sheep)
res <- ceps(sheep,f=8000,at=0.4,wl=1024,plot=FALSE)
fpeaks(res, nmax=4, xlab="Quefrency (s)")
# melscale
require(tuneR)
mel <- melfcc(sheep, nbands = 256, dcttype = "t3", fbtype = "htkmel", spec_out=TRUE)
melspec.mean <- apply(mel$aspectrum, MARGIN=2, FUN=mean)
melspec.mean <- melspec.mean/max(melspec.mean) # [0,1] scaling 
fpeaks(melspec.mean, nmax=4, f=8000, mel=TRUE)
fpeaks(melspec.mean, freq=4, f=8000, mel=TRUE) # freq in Hz!
fpeaks(melspec.mean, threshold=0.3, f=8000, mel=TRUE)
fpeaks(melspec.mean, amp=c(0.1,0.1), f=8000, mel=TRUE)

Fourier transform windows

Description

Generates different Fourier Transform windows.

Usage

ftwindow(wl, wn = "hamming",
        correction = c("none", "amplitude", "energy"))

Arguments

Value

A vector of length wl.

Note

Try the example to see windows shape.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Harris, F.J., 1978. On the use of windows for harmonic analysis with the discrete Fourier Transform. Proceedings of the IEEE, 66(1): 51-83.

See Also

covspectro, dfreq, meanspec, spec, spectro, spectro3D

Examples

a<-ftwindow(512)
b<-ftwindow(512,wn="bartlett")
c<-ftwindow(512,wn="blackman")
d<-ftwindow(512,wn="flattop")
e<-ftwindow(512,wn="hanning")
f<-ftwindow(512,wn="rectangle")
all<-cbind(a,b,c,d,e,f)
matplot(all,type="l",col=1:6,lty=1:6)
legend(legend=c("hamming","bartlett","blackman","flattop","hanning","rectangle"),
x=380,y=0.95,col=1:6,lty=1:6,cex=0.75)

Fundamental frequency track

Description

This function estimates the fundamental frequency through a short-term cepstral transform.

Usage

fund(wave, f, channel = 1, wl = 512, ovlp = 0, fmax = f/2, threshold = NULL,
at = NULL, from = NULL, to = NULL, 
plot = TRUE, xlab = "Time (s)", ylab = "Frequency (kHz)",
ylim = c(0, f/2000), pb = FALSE, ...)

Arguments

.

Value

When plot is FALSE, fund returns a two-column matrix, the first column corresponding to time in seconds (x-axis) and the second column corresponding to to fundamental frequency in kHz (y-axis).
NA corresponds to pause sections in wave (see threshold).
No plot is produced when using at.

Note

This function is based on ceps.

Author(s)

Jerome Sueur sueur@mnhn.fr.

References

Oppenheim, A.V. and Schafer, R.W. 2004. From frequency to quefrency: a history of the cepstrum. Signal Processing Magazine IEEE, 21: 95-106.

See Also

cepstro, ceps, autoc

Examples

data(sheep)
# estimate the fundamental frequency at a single position
fund(sheep, f=8000, fmax=300, at=1, plot=FALSE)
# track the fundamental frequency along time
fund(sheep,f=8000,fmax=300,type="l")
# with 50% overlap between successive sliding windows, time zoom and 
# amplitude filter (threshold)
fund(sheep,f=8000,fmax=300,type="b",ovlp=50,threshold=5,ylim=c(0,1),cex=0.5)
# overlaid on a spectrogram
spectro(sheep,f=8000,ovlp=75,zp=16,scale=FALSE,palette=reverse.gray.colors.2)
par(new=TRUE)
fund(sheep,f=8000,fmax=300,type="p",pch=24,ann=FALSE,
  xaxs="i",yaxs="i",col="black",bg="red",threshold=6)

Gammatone filter

Description

Generate gammatone filter in the time domain (impulse response).

Usage

gammatone(f, d, cfreq, n = 4, a = 1, p = 0, output = "matrix")

Arguments

Details

The gammatone function in the time domain (impulse response) is obtained with:

g(t) = a \times t^{n-1} \times e^{-2\pi \beta t} \times \cos(2 \pi cf t + \phi)

with a the amplitude, t time, n the filter order, cf the center frequency, \phi the initial phase.
The parameter \beta is the equivalent rectangular bandwidth (ERB) bandwidth which varies according to the center frequency cf following:

\beta = 24.7 \times (4.37 \times \frac{cf}{1000} + 1)

Value

A wave is returned. The class of the returned object is set with the argument output.

Note

Use the FFT based function, as spec or meanspec, to get the filter in the frequency domain. See examples.

Author(s)

Jerome Sueur

References

Holdsworth J, Nimmo-Smith I, Patterson R, Rice P (1988) Implementing a gammatone filter bank. Annex C of the SVOS Final Report: Part A: The Auditory Filterbank, 1, 1-5.

See Also

melfilterbank

Examples

## gammatone filter in the time domain (impulse response)
f <- 44100
d <- 0.05
res <- gammatone(f=f, d=d, cfreq=440, n=4)
## time display
oscillo(res, f=f)
## frequency display
spec(res, f=f)
## generate and plot a bank of 32 filters from 500 to 10000 Hz
n <- 32
cfreq <- round(seq(500, 10000, length.out=n))
res <- matrix(NA, nrow=f*d/2, ncol=n)
for(i in 1:n){
    res[,i] <- spec(gammatone(f=f, d=d, cfreq=cfreq[i]), f=f, dB="max0", plot=FALSE)[,2]
}
x <- seq(0,f/2,length.out=nrow(res))/1000
plot(x=x, y=res[,1],
     xlim=c(0,14), ylim=c(-60,0),
     type="l", col=2, las=1, 
     xlab="Frequency (kHz)", ylab="Relative amplitude (dB)")
for(i in 2:n) lines(x, res[,i], col=2)
## use the frequency domain to filter a white noise input
## here around the center frequency 2000 Hz
res <- gammatone(f=f, d=d, cfreq=2000, n=4)
gspec <- spec(res, f=f, plot=FALSE)[,2]
nw <- noisew(f=44100, d=1)
nwfilt <- fir(nw, f=44100, wl=length(gspec)*2, custom=gspec) 
spectro(nwfilt, f=f)

2D-spectrogram of a time wave using ggplot2

Description

This function returns a ggplot object to draw a spectrogram with the package ggplot2. This is an alternative to spectro.

Usage

ggspectro(wave, f, tlab = "Time (s)",
flab = "Frequency (kHz)", alab = "Amplitude\n(dB)\n", ...)

Arguments

Details

This function return the fist layer (data and aesthetic mapping) of a ggplot2 plot.
See the example section to understand how to build a spectrogram and consult ggplot2 help to get what you exactly need.
There is no way to plot the oscillogram as spectro does.

Value

A ggpot layer.

Note

This function requires ggplot2 package.

Author(s)

Jerome Sueur

References

Wickham H (2009) – ggplot2: elegant graphics for data analysis. UseR! Springer.

See Also

spectro, spectro3D, dynspec

Examples

## Not run: 
require(ggplot2)
## first layer
v <- ggspectro(tico, ovlp=50)
summary(v)
## using geom_tile ##
v + geom_tile(aes(fill = amplitude)) + stat_contour()
## coordinates flip (interest?)
v + geom_tile(aes(fill = amplitude)) + stat_contour() + coord_flip()
## using stat_contour ##
# default (not nice at all)
v + stat_contour(geom="polygon", aes(fill=..level..))
# set up to 30 color levels with the argument bins
(vv <- v + stat_contour(geom="polygon", aes(fill=..level..), bins=30))
# change the limits of amplitude and NA values as transparent
vv + scale_fill_continuous(name="Amplitude\n(dB)\n", limits=c(-30,0), na.value="transparent")
# Black-and-white theme
(vv + scale_fill_continuous(name="Amplitude\n(dB)\n", limits=c(-30,0),
  na.value="transparent", low="white", high="black") + theme_bw())
# Other colour scale (close to spectro() default output)
v + stat_contour(geom="polygon", aes(fill=..level..), bins=30)
  + scale_fill_gradientn(name="Amplitude\n(dB)\n", limits=c(-30,0),
  na.value="transparent", colours = spectro.colors(30))

## End(Not run)

Hilbert transform and analytic signal

Description

This function returns the analytic signal of a time wave through Hilbert transform.

Usage

hilbert(wave, f, channel = 1, fftw = FALSE)

Arguments

Details

The analytic signal is useful to get the amplitude envelope (see argument henv of oscillo and the instantaneous phase or frequency (see ifreq) of a time wave.

Value

hilbert returns the analytic signal as a complex matrix. The imaginary part of this matrix is the Hilbert transform.

Note

To get the Hilbert component only, use Im(Hilbert(wave)).

Author(s)

Jonathan Lees jonathan.lees@unc.edu. Implementation of 'fftw' argument by Jean Marchal and Francois Fabianek.

References

Mbu Nyamsi, R. G., Aubin, T. & Bremond, J. C. 1994 On the extraction of some time dependent parameters of an acoustic signal by means of the analytic signal concept. Its application to animal sound study. Bioacoustics, 5: 187-203.

See Also

ifreq

Examples

a<-synth(f=8000, d=1, cf=1000)
aa<-hilbert(a, f=8000)

Instantaneous frequency

Description

This function returns the instantaneous frequency (and/or phase) of a time wave through the computation of the analytic signal (Hilbert transform).

Usage

ifreq(wave, f, channel = 1, phase = FALSE, threshold = NULL,
plot = TRUE, xlab = "Time (s)", ylab = NULL,
ylim = NULL, type = "l", ...)

Arguments

Details

The instantaneous phase is the argument of the analytic signal obtained throught the Hilbert transform.
The instantaneous phase is then unwrapped and derived against time to get the instantaneous frequency.
There may be some edge effects at both start and end of the time wave.

Value

If plot is FALSE, ifreq returns a list of two components:

Note

This function is based on the analytic signal obtained with the Hilbert transform (see hilbert).
The function requires the package signal.
The matrix describing the instantaneous phase has one more row than the one describing the instantaneous frequency.

Author(s)

Jerome Sueur sueur@mnhn.fr

References

Mbu Nyamsi, R. G., Aubin, T. & Bremond, J. C. 1994 On the extraction of some time dependent parameters of an acoustic signal by means of the analytic signal concept. Its application to animal sound study. Bioacoustics, 5: 187-203.

See Also

hilbert, zc

Examples

# generate a sound with sine and linear frequency modulations
a<-synth(d=1, f=8000, cf=1500, fm=c(200,10,1000,0,0))
# plot on a single graphical device the instantaneous frequency and phase
op<-par(mfrow=c(2,1))
ifreq(a,f=8000,main="Instantaneous frequency")
ifreq(a,f=8000,phase=TRUE,main="Instantaneous phase")
par(op)

Inverse of the short-term Fourier transform

Description

This function returns a wave object from a complex STFT matrix by computing the inverse of the short-term Fourier transform (STFT)

Usage

istft(stft, f, wl, ovlp=75, wn="hanning", output = "matrix")

Arguments

Details

The function is based on the inverse of the FFT (see fft) and on the overlap add (OLA) method.
The overlap percentage must satisfy the Perfect Reconstruction OLA-constraint. For the most windows, this constraint is:

ovlp = 100 \times (1 - \frac{1}{4 \times n}),

with n being a positive integer.
A default value is set to 75%. We suggest not to change it.

Value

A new wave is returned. The class of the returned object is set with the argument output.

Note

The stft input data must be complex.
This function is used by ffilter, lfs to respectively filter in frequency and shift in frequency a sound.
The function can be used to reconstruct or modify a sound. See examples.

Author(s)

Original Matlab code by Hristo Zhivomirov (Technical University of Varna, Bulgaria), translated and adapted to R by Jerome Sueur

See Also

spectro, ffilter, lfs

Examples

## Not run: 
# STFT and iSTFT parameters
wl <- 1024
ovlp <- 75
# reconstruction of the tico sound from the stft complex data matrix
data(tico)
data <- spectro(tico, wl=wl, ovlp=ovlp, plot=FALSE, norm=FALSE, dB=NULL, complex=TRUE)$amp
res <- istft(data, ovlp=ovlp, wn="hanning", wl=wl, f=22050, out="Wave")
spectro(res)
# a strange frequency filter
n <- noisew(d=1, f=44100)
data <- spectro(n, f=44100, wl=wl, ovlp=ovlp, plot=FALSE, norm=FALSE, dB=NULL, complex=TRUE)$amp
data[64:192, 6:24] <- 0 
nfilt <- istft(data, f=8000, wl=wl, ovlp=ovlp, output="Wave")
spectro(nfilt, wl=wl, ovlp=ovlp)

## End(Not run)

Itakuro-Saito distance

Description

Compare two distributions (e.g. two frequency spectra) by computing the Itakuro-Saito distance

Usage

itakura.dist(spec1, spec2, scale=FALSE)

Arguments

Details

The Itakura-Saito (I-S) distance is a non-symmetric measure of the difference between two probability distributions. It is here adapted for frequency spectra. The distance is asymmetric, ie computing the I-S distance between spec1 and spec2 is not the same as computing it between spec2 and spec1. A symmetry can be obtained by calculating the mean between the two directions.
The distance is obtained following:

D_{I-S}(spec1 \Vert spec2) = \sum{\frac{spec1}{spec2} - log(\frac{spec1}{spec2}) - 1}

Value

The function returns a list of three items:

If scale = TRUE the distance is divided by the length of spec1 (or spec2).

Note

The function works for both Hz and (htk-)mel scales.

Author(s)

Jerome Sueur, improved by Laurent Lellouch

See Also

kl.dist, ks.dist, logspec.dist, simspec, diffspec

Examples

# Comparison of two spectra
data(tico)
tico1 <- spec(tico, at=0.65, plot=FALSE)
tico2 <- spec(tico, at=1.1, plot=FALSE)
itakura.dist(tico1, tico2)
itakura.dist(tico1, tico2, scale=TRUE) 

Kullback-Leibler distance

Description

Compare two distributions (e.g. two frequency spectra) by computing the Kullback-Leibler distance

Usage

kl.dist(spec1, spec2, base = 2)

Arguments

Details

The Kullback-Leibler distance or relative entropy is a non-symmetric measure of the difference between two probability distributions. It is here adapted for frequency spectra. The distance is asymmetric, ie computing the K-L distance between spec1 and spec2 is not the same as computing it between spec2 and spec1. A symmetry can be obtained by calculating the mean between the two directions.
The distance is obtained following:

D_{K-L}(spec1 \Vert spec2) = \sum{spec1 \times log(\frac{spec1}{spec2})}

Value

The function returns a list of three items:

Note

The base of the logarithm can be changed using the argument base. When sets to base 2, the information is measured in units of bits. When sets to base e, the information is measured in nats.
The function works for both Hz and (htk-)mel scales.

Author(s)

Jerome Sueur, improved by Laurent Lellouch

References

Kullback, S., Leibler, R.A. (1951). On information and sufficiency. Annals of Mathematical Statistics, 22: 79-86

See Also

ks.dist, logspec.dist, simspec, diffspec

Examples

# Comparison of two spectra
data(tico)
tico1 <- spec(tico, at=0.65, plot=FALSE)
tico2 <- spec(tico, at=1.1, plot=FALSE)
kl.dist(tico1, tico2)    # log2 (binary logarithm)
kl.dist(tico1, tico2, base=exp(1))  # ln (natural logarithm)

Kolmogorov-Smirnov distance

Description

This function compares two distributions (e.g. two frequency spectra) by computing the Kolmogorov-Smirnov distance

Usage

ks.dist(spec1, spec2, f = NULL, mel = FALSE,
plot = FALSE, type = "l",
lty = c(1, 2), col = c(2, 4), 
flab = NULL, alab = "Cumulated amplitude",
flim = NULL, alim = NULL,
title = TRUE, legend = TRUE, ...)

Arguments

Details

The Kolmogorov distance is the maximal distance between the cumulated spectra. The function returns this distance and the corresponding frequency. This is an adaptation of the statistic computed by the non-parametric Kolmogorov-Smirnov test (see ks.test).

Value

The function returns a list of two items

Note

There is no p-value associated to the K-S distance.
If no frequency is provided, only the distance D.