41

There is a wide variety of techniques for non-uniform FFT, and the most efficient ones are all meant for exactly your case: quasi-uniform samples. The basic idea is to smear the unevenly sampled sources onto a slightly finer ("oversampled") uniform grid though local convolutions against Gaussians. A standard FFT can then be run on the oversampled uniform ...


15

For integrations with many variables, the Monte Carlo method usually is a decent fit. Its error decreases as $O (\sqrt{N})$ where N is the number of equidistributed points selected. Of course this is not good for low dimension (1D and 2D) spaces where high order methods exist. Most of these deterministic methods, however, take a large number of points in ...


11

Have you thought of using Barycentric Interpolation? A detailed description on how to do it efficiently for Chebyshev nodes is given in Section 5 of this paper. This is actually an exact evaluation of the Chebyshev interpolant. If you're evaluating a polynomial of degree $n$ at $m$ nodes, the cost is in $\mathcal O(nm)$. Update Another alternative, if you ...


11

Your matrix $A$ isn't a circulant matrix- it's just Toeplitz. Furthermore, your $a$ vector doesn't have the "-1" in it anywhere, so you clearly don't have sufficient information. The standard remedy for this is to embed $A$ in a larger matrix which actually is a circulant matrix by padding the vectors with 0's and correctly handling the "-1" (which ...


10

Let me first answer all the questions: What is the theoretical convergence rate for an FFT Poison solver? The theoretical convergence is exponential as long as the solution is sufficiently smooth. How fast should this energy converge? The Hartree energy $E_H$ should converge exponentially for a sufficiently smooth solution. If the solution is less ...


10

In signal processing, aliasing is avoided by sending a signal through a low pass filter before sampling. Jack Poulson already explained one technique for non-uniform FFT using truncated Gaussians as low pass filters. One inconvenient feature of truncated Gaussians is that even after you have decided on the grid spacing for the FFT (=the sampling rate in ...


9

17. Lots and lots of work has gone into good fft implementations and it's unlikely you'll be able to reliably outperform a good fft library. For example, fftw "automatically adapts itself to your machine, your cache, the size of your memory, the number of registers, and all the other factors that normally make it impossible to optimize a program for more ...


9

This is a linear integral equation (because the right hand side is linear in $h$) and presumably the problem is ill-posed (for one, because of the multiplication by $\sin(\theta)$ but also because the "integral kernel" is smooth). Methods to solve this (approximately) are general Galerkin methods or collocation methods. For Galerkin methods you take a ...


8

This is more anecdotal evidence than demonstrated proof, but it appears that existing implementations for FFT's, such as FFTW, have a limit to their scaling ability. When we started using LAMMPS's ${\bf k}$-space solvers in very large systems ($O(10^7)$ atoms), we found that the scaling continued, so long as we were able to keep the number of processors ...


8

Sparse grid quadrature is an alternative approach to integrate in higher dimensions. Quadrature relies on evaluating a weighted sum of function values at specific "optimal" points. Traditional quadrature uses a tensor product grid construction in higher dimensions, which means that you would have to evaluate the function at an exponentially growing ...


8

I see several issues: The DFT computed with fft puts the zero mode at the beginning of the array, and if you want to compute the derivative, it is necessary to apply fftshift/ifftshift to the array N to make sure the derivative is correct. It is easy to see for yourself what the correct expression is by working it out with pen and paper, and see also the ...


7

Michael Pippig at the University of Chemnitz (Germany) has implemented an MPI-parallelized FFT that uses FFTW in the background. This might help you: http://www-user.tu-chemnitz.de/~mpip/software.php?lang=en It is using the algorithm proposed by Plimpton from Sandia National Labs as suggested by Eldila's comment.


7

The frequency content of the interpolated signal is significantly influenced by the interpolation basis. If you have a band-limited function that you have adequately sampled (i.e. satisfying Nyquist criterion), interpolating with any function that is not band-limited to the same frequency will indeed introduce high-frequency noise. Unfortunately, exact band-...


6

The fast multipole method (FMM) is $O(n)$ and has much lower communication requirements, so it provides a highly scalable discrete Fourier transform. Edelman, McCorquodale, and Toledo (1999) The Future Fast Fourier Transform? analyze this approach and conclude that FMM would be preferable to conventional FFT at large scale. Note that the FMM is only ...


6

If you are interested in software I can recommend the NFFT library (in C with an interface to MATLAB) which can be found here. Note that there is also a PFFT library for parallel FFT computation and even a PNFFT library for parallel non-equispaces FFTs by the same developers.


6

In a word, yes. A least squares problem arises in the context of an over determined system. That means that you have more equations than unknowns and would be the case if you had more data points than Fourier coefficients that you're trying to find. This would give you the "best fit" Fourier series of a given order. For example, you could find the best fit ...


6

This turned out much longer than I planned so I hope it is useful and can be a resource for others with the same questions. I believe you mean to ask primarily about the discrete Fourier transform (in which case your definition (2) for $\hat{f}_\alpha$ is not right, for the DFT it should be $\hat{f}_\alpha = \sum_{j=0}^{2M-1} x_j \exp(- i\alpha\frac{\pi j}{...


6

You are solving Poisson's equation: $$\nabla^2 \Phi = 4\pi G \rho.$$ Notice that if $\Phi$ is a solution, then so is $\Phi+C$ for any constant $C$. Furthermore, $C$ will have no effect on your simulation since the forces depend only on derivatives of $\Phi$. The divide-by-zero issue is just a manifestation of this degeneracy. In the Fourier-...


6

From your link, consider the definition of the round off error and the statement "Since the exact solution must satisfy the discretized equation exactly, the error must also satisfy the discretized equation." This is actually only true if the PDE is homogeneous, that is, if we can write it in the form $\mathcal{L}(u;u_t,u_x,u_{xx},\ldots)=0$, with all terms ...


5

Yes, it is true, the CVX package (disclosure: I am the author) does not support the use of the fft() command in constraint and objective expressions. However, it's relatively simple to get the equivalent result, in both the 1D and 2D cases. First, 1D. Given a vector $x\in\mathbb{R}^n$, the DFT of $x$ is equal to $Wx$, where $W$ is the so-called DFT matrix. ...


5

In the context of PDEs, it's important to recognize that the value of $n$ for a required 1D FFT will typically grow like the $d$th root of the total number of grid points, where the dimensionality $d$ is most often 3. Searching for "parallel FFT" or "pseudospectral scalability" on Google Scholar yields a wealth of information that I'm unqualified to assess. ...


5

There are several papers by Markus Püschel on his web site here that discuss Cooley-Tukey-like (so I'm guessing "fast") algorithms for lattice transforms, such as DFTs on triangular and hexagonal 2-D lattices. In the triangular case, he calls the DFT the discrete triangle transform (DTT). Markus has a code called SPIRAL that automatically generates code for ...


5

This is pure haar scaling function approximation of $f(t)$ $$f(t)=\sum_{k=-\infty}^{\infty}a_k\phi(2^jt-k)$$ where $$\phi(t) = \begin{cases}1 \quad; 0 \leq t < 1,\\0\quad;\mbox{otherwise.}\end{cases}$$ is haar scaling function, $j$ scaling aparameter, $k$ translation parameter and the coefficients $$a_k=\int f(t)\phi(2^jt-k)dt$$ So on our case if $f(t)...


5

The easiest way to generate a divergence-free 2D function is to use a streamfunction, $\psi$. Where $$ u=\frac{\partial\psi}{\partial y}\quad \rm{and}\quad v=-\frac{\partial \psi}{\partial x}$$ Once you generate a random $\psi$, you can then use your Fourier representation (or finite differences) to compute the two derivatives. I haven't tried this kind of ...


5

You can write a curve as a parametric equation ($x(t)$, $y(t)$) for a range say $0 < t < 1$. Then you can have a Fourier decomposition of both $x(t)$ and $y(t)$. A rectangle can be written as such a parametric equation. However, you want to decompose a given function as a sum of functions from a specific class. There is a whole field devoted to this. ...


5

1) You have misinterpreted the order the fft2 output elements are stored. This is a very common mistake, since it is very easy to get confused. Thankfully, there is an easy way to avoid such bugs. Arrange the Fourier space vectors the "natural way" % Fourier space vectors kx = 2*pi/Lx*[-Nx/2:Nx/2-1]'; ky = 2*pi/Ly*[-Ny/2:Ny/2-1]'; Then, make sure the fft2 ...


5

The point is to guarantee that the Fourier coefficients exist, not to make really fine distinctions about what might happen if the guarantee is violated. Furthermore, you would be hard pressed to find a reasonable function for which this is an issue, if it exists at all. The class of integrable functions is already so wide as to include everything you might ...


5

Running your code, it seems like your pulse looks kinda like this: (sorry for not adding units to the plots, I used the same as you, i.e. t is in fs and w in rad/fs) So, the FWHM is not correct (should be 4 fs) and the angular carrier/central frequency is messed up (should be ca. 2.3 / ca. 0.37 per fs). I think there is two things not quite right here: 1) ...


4

As an alternative you could look into using the Goertzel Algorithm to directly compute the frequency components you are interested in.


Only top voted, non community-wiki answers of a minimum length are eligible