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 ...


8

You may find the expansion of a plane wave in spherical waves to be helpful here: $$ e^{i\mathbf{k}\cdot\mathbf{x}} = 4\pi\sum_{l=0}^\infty\sum_{m=-l}^l i^lj_l(kr)Y_{lm}(\theta,\phi)Y_{lm}^*(\vartheta,\varphi) $$ where $\theta$, $\phi$ are the angular variables for $\mathbf{x}$ and $\vartheta$, $\varphi$ for $\mathbf{k}$; the radial functions $j_l$ are the ...


8

FFT returns a complex array that has the same dimensions as the input array. The output array is ordered as follows: Element 0 contains the zero frequency component, F0. The array element F1 contains the smallest, nonzero positive frequency, which is equal to 1/(Ni Ti), where Ni is the number of elements and Ti is the sampling interval. F2 corresponds to ...


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-...


7

Expanding upon my comment above, if you only need a few digits of accuracy you can probably use the method of stationary phase. We can follow the procedure on Wikipedia. We can write the transform as follows: $$ f(t) = \int_{-\infty}^{\infty} g(\omega)e^{i\phi(\omega,t)}d\omega, $$ where $\phi(\omega,t)=C\omega^2-\omega t$. The principal contribution of the ...


7

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) ...


6

There are three issues that are likely to cause such problems in pseudospectral methods: Gibbs oscillations Aliasing Time step too large In any case you likely develop oscillations in the solution until some point ends up with a negative density, resulting in a NaN when computing the pressure or sound speed or some other term. The solution to 3 is obvious, ...


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

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-...


5

Rather than having us debug your homework for you, I recommend that you read this blog post by Jake Vanderplas that contains a wonderful explanation of the FFT algorithm along with a highly readable Python implementation.


5

Let's suppose you have a function $f$ in $L^{2}$, and for simplicity, let's suppose it's a single-variable function defined over the whole real line. If its Fourier transform $\mathcal{F}(f)$ is $\hat{f}$, then $\mathcal{F}(f')(\xi) = 2\pi i \xi\hat{f}(\xi)$. If you're calculating the $n$th derivative, iterating the process above yields $\mathcal{F}(f^{(n)})...


5

That is not very difficult actually. Let $\Lambda\subset\Bbb R^n$ be the lattice on which you sampled your points, and let $\Phi:\Bbb R^n\to\Bbb R^n$ be a linear transformation under which $\Bbb Z^n\to\Lambda$ isomorphically (i.e. the columns of the matrix of $\Phi$ generate $\Lambda$). Let $f:\Bbb R^n\to\Bbb C$ be the function you are interested in, and of ...


5

The interpolation indeed affects the Fourier transform. @Steve already gives the correct answer in general, but I want to give you an example that helps the intuition more. Think for example that you have samples a sine function on a set of equally spaced points. If you did a (discrete) FFT on these points alone, you would of course recover a Fourier ...


4

A brute-force approach is to pad the signal and the kernel with enough zeros that there would be no aliasing when doing circular convolution: AIFFourier = Fourier[PadRight[AIF, 1024]]; kernelCurve = FillArray[0.005, 1, mat]; kernelCurveFourier = Fourier[PadRight[kernelCurve, 1024]]; resSource = InverseFourier[AIFFourier*kernelCurveFourier]; ListPlot[...


4

This is part of the (complex-valued) Fourier transform for which there is (provably) no more efficient way than the Fast Fourier Transform (FFT) if you want to compute the integrals for at least a significant fraction of all frequencies $n$ between zero and the largest frequency you care about. In other words, while your intuition may tell you that ...


4

For large values of $k$, the Hankel transform depends primarily on $f(r)$ at small values of $r$. Your smallest values of $r$ aren't small enough to get good coverage. For example at $k=1 \times 10^{5}$, and your smallest value of $r=4 \times 10^{-5}$, $kr=4$, so you've basically skipped over the first cycle of the Bessel function. Since you know that ...


4

It looks to me that the FFT convolution algorithm is doing what is expected here. Remember that you work with discretized signals, and therefore, discretized convolution. There are a few ways to implement the discrete convolution numerically. The slight differences basically come from the way to deal with the finiteness of each signal. The discrete ...


4

Let $f(\omega)$ be your power spectrum. Then maybe something like $$ \frac{\|f\|_{L^\infty}\|f\|_{L^0}}{\|f\|_{L^1}} = \frac{\mathrm{max}_{\omega\in\Omega} f(\omega)\cdot|\omega_{max}-\omega_{min}|}{\int_\Omega|f(\omega)|d\omega}. $$ I know that $L^0$ isn't really good notation but I think it is useful for presenting this. This quantity is minimized by ...


4

Maxim Umansky’s answer describes the storage convention of the FFT frequency components in detail, but doesn’t necessarily explain why the original code didn’t work. There are three main problems in the code: x = linspace(0,2*pi,N): By constructing your spatial domain like this, your x values will range from $0$ to $2\pi$, inclusive! This is a problem ...


3

Seems like I'm a bit late, but I was just answering a similar question. The first thing here, is that you probably want to pay attention to the definition of the w-axis. The FFT-algorithm will consider your signal to start at w=0, so the way you defined frequencies, your carrier frequency will be wrong. But back to your question: the websites you used to ...


3

You might want to take a look at the Eigen C++ matrix class library. http://eigen.tuxfamily.org/index.php?title=Main_Page There is an FFT class in the unsupported section of the library but my impression is that it is relatively mature. Here is a snippet of code using this FFT class with the BOOST multi-precision classes. #include <vector> #...


3

It's a question of function spaces. If your function is in $L^2$, then you know that you can approximate the function using a Fourier series because the Fourier basis (i.e., the sines and cosines) are dense in $L^2$. In other words, you know that if you take more and more terms of the series, then you will approximate the function better and better. On the ...


3

The accepted answer above is very misleading and @nat-chouf has it correct. In the spirit of the original question, running a pseudo-spectral mode of isotropic turbulence, then zero-ing out $k>\frac{2}{3}k_\textrm{max}$ at each time-step will absolutely remove energy from the model, assuming that energy is moved to these higher wavenumbers.


3

I recommend you first read answer here to find out why $\tilde{E}(\mathbf{k}) \neq \frac{1}{2} \tilde{\mathbf{u}}(\mathbf{k}) \cdot \tilde{\mathbf{u}}(\mathbf{k})$. In order to find velocity profile from Fourier transform ($\mathscr{F}$) of kinetic energy ($E(\mathbf{r}) = \frac{1}{2} \mathbf{u}(\mathbf{r}) \cdot \mathbf{u}(\mathbf{r})$), you could obtain it ...


2

Isn't this the same as the real part of the Fourier transform of f(x) $$ \begin{align} I(n) &= \int^1_0dxf(x)\cos(n\pi x) \\ &=\Re{{\int^1_0dx f(x) e^{-j n \pi x}}} \\ &= \Re\int^1_0dx f(x) (cos(n\pi x) - j\sin(n \pi x)) \\ &= \int^1_0dxf(x)\cos(n\pi x) \end{align} $$


2

If $f$ is very nice, then you can approximate $f$ with a sequence of piecewise polynomials and integrate over the resulting intervals exactly. This may be much, much cheaper than using an FFT. This is also true for approximations of $f$ by any functions where you can easily know or determine the exact antiderivative in each interval. If $f$ is a black box, ...


2

The problem can be partially solved following this tutorial. Given Fisher-Kolmogorov equation \begin{equation} u_t=u_{xx}+u(1-u) \end{equation} It can also be written as \begin{equation} u_t = u_{xx} + u\,v\\ v_t = v_{xx} - u\,v \end{equation} where $v = (1-u)$ Solving the equation \begin{equation} u_t = u_{xx} \\ \end{equation} from $t$ to $\...


2

I think you are looking for mpi4py-fft, which is a Python package (BSD-2 licensed) with its wrappers on the serial FFTW library. From pretty extensive mpi4py documentation: Parallel FFTs are computed through a combination of global redistributions and serial transforms. The library also depends on mpi4py (developed also by one of the co-authors of mpi4py-...


2

I have tried the following code in Julia for $f(r) = exp(-r)$ with 1D FFT, which seems to be working somehow. So, could you compare your code with it and see if there is some difference of normalization? (Please note that fft() below makes no normalization, i.e., $1/sqrt(N)$ is not present, while the FFT routine in NAG library may be multiplying by some ...


2

Parseval's theorem tells you that the dot product between two vectors equals the dot product of the Fourier transforms, possibly up to a constant. Consequently, you do not need to transform back the result of $a\ast b$ and $c \ast d$ and should be able to do the overall operation with just 4 FFTs.


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