ndft: Non-uniform Discrete Fourier Transform

Description

Compute the one-dimensional Non-uniform Discrete Fourier Transform (NDFT). This is needed if the data is sampled in non-uniform intervals.

Usage

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

Arguments

vector of real or complex function values

vector of points in direct space, typically time or position coordinates. If inverse=FALSE, x must have the same length as f.

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.

Value

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

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.

Examples

Run this code

# NOT RUN {
# 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)

# }