ndft: Non-uniform Discrete Fourier Transform

In obreschkow/cooltools: Practical Tools for Scientific Computations and Visualizations

Non-uniform Discrete Fourier Transform

Description

Compute the one-dimensional Non-uniform Discrete Fourier Transform (NDFT). This is needed if the data are sampled at irregularly spaced times.

Usage

ndft(
  f,
  x = seq(0, length(f) - 1)/length(f),
  nu = seq(0, length(f) - 1),
  inverse = FALSE,
  weighing = TRUE,
  simplify = TRUE
)

Arguments

f	vector of real or complex function values
x	vector of points in direct space, typically time or position coordinates. If inverse=FALSE, x must have the same length as f.
nu	vector of frequencies in units of [1/units of x]. If inverse=TRUE, nu must have the same length as f.
inverse	logical flag; if TRUE, the inverse Fourier transform is performed.
weighing	logical flag; if TRUE, irregularly spaced evaluations of f will be weighted proportionally to their bin width in x (if inverse=FALSE) or nu (if inverse=FALSE).
simplify	logical flag; if TRUE, the complex output array will be simplified to a real array, if it is real within the floating point accuracy.

Details

The one-dimensional NDFT of a vector f=(f_1,...,f_N) is defined as

F_j=\sum_i w_i f_i exp(-2\pi i*x_i*\nu_j)

where w_i are optional weights, proportional to the interval around x_i, only used if weighing=TRUE. Likewise, the inverse NDFT is defined as

f_i=\sum_j w_j F_j exp(+2\pi i*x_i*nu_j)

where w_j are optional weights, proportional to the interval around \nu_j. In this implementation NDFTs are computed using a brute force algorithm, scaling as O(N*N), which is considerably worse than the O(N)*log(N) scaling of FFT algorithms. It is therefore important to pick the required frequencies wisely to minimise computing times.

Value

Returns a vector of the same length as x (if inverse=FALSE) or nu (if inverse=TRUE).

Author(s)

Danail Obreschkow

See Also

fft

Examples

# Define an example signal
nu1 = 1 # [Hz] first frequency
nu2 = 8 # [Hz] second frequency in the signal
s = function(t) sin(2*pi*nu1*t)+0.7*cos(2*pi*nu2*t+5)

# Discretize signal
N = 50 # number of samples
t.uniform = seq(0,N-1)/N
t.nonuniform = t.uniform^1.3
s.uniform = s(t.uniform)
s.nonuniform = s(t.nonuniform)

# Plot signal
oldpar = par(mfrow = c(1, 2))
curve(s,0,1,500,xaxs='i',main='Time signal',xlab='Time t',ylab='s(t)',col='grey')
points(t.uniform,s.uniform,pch=16,cex=0.8)
points(t.nonuniform,s.nonuniform,pch=4,col='blue')
legend('topright',c('Continuous signal','Uniform sample','Non-uniform sample'),
       lwd=c(1,NA,NA),pch=c(NA,16,4),col=c('grey','black','blue'),pt.cex=c(1,0.8,1))

# Uniform and non-uniform DFT
nu = seq(0,N-1) # discrete frequencies
spectrum.uniform = stats::fft(s.uniform)
spectrum.nonuniform = ndft(s.nonuniform,t.nonuniform,nu)
spectrum.wrong = stats::fft(s.nonuniform)

# Evaluate power
power.uniform = Mod(spectrum.uniform)^2
power.nonuniform = Mod(spectrum.nonuniform)^2
power.wrong = Mod(spectrum.wrong)^2

# Plot DFT and NDFT up to Nyquist frequency
plot(nu,power.uniform,pch=16,cex=0.8,xlim=c(0,N/2),xaxs='i',
               main='Power spectrum',xlab=expression('Frequency'~nu~'[Hz]'),ylab='Power')
points(nu,power.nonuniform,pch=4,col='blue')
points(nu,power.wrong,pch=1,col='red')
abline(v=c(nu1,nu2),col='grey',lty=2)
legend('topright',c('DFT of uniform sample','NDFT of non-uniform sample',
'DFT of non-uniform sample (wrong)','Input frequencies'),
       lwd=c(NA,NA,NA,1),lty=c(NA,NA,NA,2),pch=c(16,4,1,NA),
       col=c('black','blue','red','grey'),pt.cex=c(0.8,1,1,NA))
par(oldpar)

obreschkow/cooltools documentation built on Nov. 16, 2024, 2:46 a.m.