Last data update: 2014.03.03

 R: Computes Short Time Fourier Transforms
stftR Documentation

Computes Short Time Fourier Transforms

Description

Processes a dataset, creating an object contained processed time-frequency analyses. These can then be plotted.

Usage

stft(X, start=0, end=1, length=NULL,  time.format = c("auto"), 
    type = c("mv", "svm", "sum"), mv.indices, date.col, 
    reassign = TRUE,plot.it = FALSE,...)

Arguments

X

The dataset to be processed.

start, end, length, time.format

A specification for the segment to process, as in get.intervals.

type

The type of STFT to compute.

mv.indices

For type = "mv" or type = "sum", the indices to process and the order to process them in.

date.col

logical. Whether the first column should be ignored and treated as a timestamp. If unset, is automatically chosen.

reassign

logical. If TRUE, compute the time-reassigned STFT. For type %in% c("mv", "sum"), this is done with the first coordinate in mv.indices.

plot.it

logical. Whether to plot the STFT immediately when processing is complete, using the default plot.stft options.

...

Additional optional arguments to control the STFT computation. These are:

win: Window size in seconds for STFT computation. Increased window size mean better frequency resolution, but poorer time resolution. Defaults to 10 seconds.

inc: Increment between successive time steps for processing. Defaults to win/2.

coef: Number of fourier frequencies to compute. Small values will remove the higher frequencies from the processed object. Defaults to the maximum, win/2.

wtype: String giving the name of a window function, providing coefficients for filtering before processing. "hanning.window" is the default, with "uniform.window" also available.

freq: Sampling frequency of data set. If not given, is taken from X itself, or assumed to be 1 if unavailable.

center: If TRUE (Default), center the data in each window before processing is done. Useful for avoiding excessively large DC offset coefficients in results.

calc.null: If TRUE (Defaults to FALSE), compute a 'null' STFT by resampling the data completely, then doing a STFT.

pvalues: If TRUE (Defaults to FALSE) Compute bootstrapped pvalues for each position by resampling within each window and applying a wilcox test.

time: Allows the user to set an overriding timestamp vector to be used for processing.

quiet: If TRUE, suppress output.

Details

This function accepts input in a variety of forms and computes short time fourier transforms to extract frequency structure from the data.

X may be an array, a vector, or an AccData object. If date.col is TRUE, the first column of an array X would be used to determine timestamps. Otherwise indices would be used. If date.col is not set, the function will attempt to determine whether the first column is timestamp-like. The timestamp column is removed from X (and so not included in consideration of mv.indices, for instance).

With vectors, the basic method is to compute the STFT by creating windows of size win seconds every inc seconds, and computing the fourier transform. With multi-dimensional data and AccData, processing is done on the dimensions that are in mv.indices, or the first three non-date columns if that is unavailable. Three methods are possible:

1. type = "mv": The one dimensional method is first applied to each of the chosen column indices. These are then collated by choosing, for each time-frequency combination, the maximum such value across each of the indices.

2. type = "svm": The SVM is computed first for each time step by computing the square rooted sum of squares. This is then dealt with using the one dimensional method.

3. type = "sum": As in "mv", the 1d method is applied. The square of the modulus of the result is then summed and square rooted.

If reassign is set, the time reassigned stft is also computed for the first element of mv.indices or the svm as appropriate, by using finite differencing. This gives potentially better resolution results for data with a clear signal component.

Value

A "stft" class object - a list with the following components:

call

The function call.

type

Type of STFT computed.

values

Mod of FFT computed, with each row corresponding to a specific time increment.

increment, windowsize, center, sampling.frequency

Various control parameters used in the computation.

null.logmean, null.logsd

Log of the square rooted mean and standard deviation of the Mod FFT squared for the randomised data, if calc.null = TRUE.

p.values

Wilcoxian pvalues, if pvalues = TRUE.

principals

Principal frequencies.

frequency

Frequencies at which FFT is computed.

time

Timestamps for FFT windows.

LGD

Local group delay matrix for reassigned STFT.

CIF

Channelized instantaneous frequency matrix for reassigned STFT.

Acknowledgements

The initial implementation of this function is based on that in package e1071. The reassigned STFT implementation is due to Nelson (2001), via Fulop (2006).

References

Fulop, S.A. & Fitz, K. (2006). Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications J Acoustical Society of America 119(1), 360–371.

Nelson. D.J. (2001). Cross-spectral methods for processing speech J Acoustical Society of America 110(1), 2575-2592.

See Also

plot.stft, get.intervals

Examples

#Some artificial data
time = 1:5000
#sum of two sine curves at 0.3 Hz and 0.05 Hz
f1 = 0.3; f2 = 0.05
sin1 = sin(time * f1 * 2*pi)
sin2 = sin(time * f2 * 2*pi)
#add a bit of noise
signal = sin1 + sin2 + 1*rnorm(5000)
#non-reassigned
stft(signal, plot = TRUE, reassign = FALSE, win = 100)
#reassigned
stft(signal, plot = TRUE, reassign = TRUE, win = 100)

#add a third component: varying frequency.
stft(signal + sin( cumsum(seq(f2, f1, length = 5000))*2*pi), plot = TRUE, reassign = TRUE, win = 100)

#Real data
binfile  = system.file("binfile/TESTfile.bin", package = "GENEAread")[1]

#Read in the entire file, calibrated
procfile<-read.bin(binfile)
#Default is mv
stft(procfile, plot.it = TRUE)
#Try sum?
stft(procfile, plot.it = TRUE, type = "sum", reassign = FALSE)

#Just look at the last 50% of the data
stft(procfile, start = 0.5, plot.it = TRUE)

#not reassigned, svm
stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE)
#a narrower 5 second window means better time resolution
stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE, win = 5)
#choose increments so as not to overlap
stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE, win = 5, inc = 5)
#uniform windows
stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE, wtype = "uniform.window")

#Svm, reassigned, quietly
obj = stft(procfile, type = "svm", quiet = TRUE)
plot(obj, cex = 3, showmax = FALSE, mode = "pval")

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(GENEAread)
Loading required package: bitops
GENEAread 1.1.1 loaded

> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/GENEAread/stft.Rd_%03d_medium.png", width=480, height=480)
> ### Name: stft
> ### Title: Computes Short Time Fourier Transforms
> ### Aliases: stft
> ### Keywords: ts
> 
> ### ** Examples
> 
> #Some artificial data
> time = 1:5000
> #sum of two sine curves at 0.3 Hz and 0.05 Hz
> f1 = 0.3; f2 = 0.05
> sin1 = sin(time * f1 * 2*pi)
> sin2 = sin(time * f2 * 2*pi)
> #add a bit of noise
> signal = sin1 + sin2 + 1*rnorm(5000)
> #non-reassigned
> stft(signal, plot = TRUE, reassign = FALSE, win = 100)
================================================================================
STFT object:
70-01-01 00:00:50.000 (Thu)  to  70-01-01 01:22:30.000 (Thu) 
99 increments of 50 s 
Window size:  100 ( 100 s ) -> f resolution:  0.01 Hz

------ 
{ stft(X = signal, reassign = FALSE, plot.it = TRUE, win = 100) }
> #reassigned
> stft(signal, plot = TRUE, reassign = TRUE, win = 100)
================================================================================
STFT object:
70-01-01 00:00:50.000 (Thu)  to  70-01-01 01:22:30.000 (Thu) 
99 increments of 50 s 
Window size:  100 ( 100 s ) -> f resolution:  0.01 Hz
[Reassign]
------ 
{ stft(X = signal, reassign = TRUE, plot.it = TRUE, win = 100) }
> 
> #add a third component: varying frequency.
> stft(signal + sin( cumsum(seq(f2, f1, length = 5000))*2*pi), plot = TRUE, reassign = TRUE, win = 100)
================================================================================
STFT object:
70-01-01 00:00:50.000 (Thu)  to  70-01-01 01:22:30.000 (Thu) 
99 increments of 50 s 
Window size:  100 ( 100 s ) -> f resolution:  0.01 Hz
[Reassign]
------ 
{ stft(X = signal + sin(cumsum(seq(f2, f1, length = 5000)) * 2 *      pi), reassign = TRUE, plot.it = TRUE, win = 100) }
> 
> #Real data
> binfile  = system.file("binfile/TESTfile.bin", package = "GENEAread")[1]
> 
> #Read in the entire file, calibrated
> procfile<-read.bin(binfile)
Loading required package: mmap
Number of pages in binary file: 104 
Calculated page references... 
Processing...
================================================================================
Processing took: 0.067 secs .
Loaded 31200 records (Approx  2 MB of RAM)
12-05-23 16:47:50.000 (Wed)  to  12-05-23 16:53:01.990 (Wed) 
Warning message:
In FUN(newX[, i], ...) : NAs introduced by coercion
> #Default is mv
> stft(procfile, plot.it = TRUE)
Extracting time interval:  16:47:50 16:53:01 
================================================================================
================================================================================
================================================================================
STFT object:
12-05-23 16:47:54.990 (Wed)  to  12-05-23 16:52:54.990 (Wed) 
61 increments of 5 s 
Window size:  1000 ( 10 s ) -> f resolution:  0.1 Hz
[MV- 123 ][Reassign]
------ 
{ stft(X = procfile, plot.it = TRUE) }
> #Try sum?
> stft(procfile, plot.it = TRUE, type = "sum", reassign = FALSE)
Extracting time interval:  16:47:50 16:53:01 
================================================================================
================================================================================
================================================================================
STFT object:
12-05-23 16:47:54.990 (Wed)  to  12-05-23 16:52:54.990 (Wed) 
61 increments of 5 s 
Window size:  1000 ( 10 s ) -> f resolution:  0.1 Hz
[SUM- 123 ]
------ 
{ stft(X = procfile, type = "sum", reassign = FALSE, plot.it = TRUE) }
> 
> #Just look at the last 50% of the data
> stft(procfile, start = 0.5, plot.it = TRUE)
Extracting time interval:  16:50:25 16:53:01 
================================================================================
================================================================================
================================================================================
STFT object:
12-05-23 16:50:30.980 (Wed)  to  12-05-23 16:52:55.980 (Wed) 
30 increments of 5 s 
Window size:  1000 ( 10 s ) -> f resolution:  0.1 Hz
[MV- 123 ][Reassign]
------ 
{ stft(X = procfile, start = 0.5, plot.it = TRUE) }
> 
> #not reassigned, svm
> stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE)
Extracting time interval:  16:47:50 16:53:01 
================================================================================
STFT object:
12-05-23 16:47:54.990 (Wed)  to  12-05-23 16:52:54.990 (Wed) 
61 increments of 5 s 
Window size:  1000 ( 10 s ) -> f resolution:  0.1 Hz
[SVM]
------ 
{ stft(X = procfile, type = "svm", reassign = FALSE, plot.it = TRUE) }
> #a narrower 5 second window means better time resolution
> stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE, win = 5)
Extracting time interval:  16:47:50 16:53:01 
================================================================================
STFT object:
12-05-23 16:47:52.490 (Wed)  to  12-05-23 16:52:57.490 (Wed) 
123 increments of 2.5 s 
Window size:  500 ( 5 s ) -> f resolution:  0.2 Hz
[SVM]
------ 
{ stft(X = procfile, type = "svm", reassign = FALSE, plot.it = TRUE,      win = 5) }
> #choose increments so as not to overlap
> stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE, win = 5, inc = 5)
Extracting time interval:  16:47:50 16:53:01 
================================================================================
STFT object:
12-05-23 16:47:52.490 (Wed)  to  12-05-23 16:52:57.490 (Wed) 
62 increments of 5 s 
Window size:  500 ( 5 s ) -> f resolution:  0.2 Hz
[SVM]
------ 
{ stft(X = procfile, type = "svm", reassign = FALSE, plot.it = TRUE,      win = 5, inc = 5) }
> #uniform windows
> stft(procfile, type = "svm", reassign = FALSE, plot.it = TRUE, wtype = "uniform.window")
Extracting time interval:  16:47:50 16:53:01 
================================================================================
STFT object:
12-05-23 16:47:54.990 (Wed)  to  12-05-23 16:52:54.990 (Wed) 
61 increments of 5 s 
Window size:  1000 ( 10 s ) -> f resolution:  0.1 Hz
[SVM]
------ 
{ stft(X = procfile, type = "svm", reassign = FALSE, plot.it = TRUE,      wtype = "uniform.window") }
> 
> #Svm, reassigned, quietly
> obj = stft(procfile, type = "svm", quiet = TRUE)
Extracting time interval:  16:47:50 16:53:01 
> plot(obj, cex = 3, showmax = FALSE, mode = "pval")
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>


Created & Maintained by Osamu Ogasawara (osamu.ogasawara@gmail.com) and IMSBIO