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

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


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

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


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

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

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

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

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


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

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

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

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

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

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

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


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

I would need additional information about your attempt (e.g., example code) to help with your first question, but here's my input for the division-by-zero question. First, do you really need to divide by $q$? In the problems where these transforms come up, you can often get away with working with $rf(r)$ and $q\hat{f}(q)$ in their own right. Decide if you ...


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.


2

This is due to an implicit time shift, which corresponds to a phase shift, perceived here as a sign inversion. The signal you are Fourier transforming is symmetric around t=0 and this is why you should expect a positive power spectrum. However, the discrete Fourier transform (DFT) of a time-series $x_n$, which is what the FFT alorithm implements, is ...


2

Since your histogram is very sparse, your density field is approximately a sum of weighted delta functions: $$\rho(\vec r) = \sum_{j=1}^{N_{\mathrm{bins}}} w_j\, \delta(\vec r - \vec r_j)$$ where $r_j$ is the center of each bin, and $w_j$ is the weight. Its Fourier transform is then simply $$\begin{align} \hat\rho(\vec k) &= \int \exp(-i\vec k\cdot\vec ...


2

You need to pay attention that unless properly padded the Multiplication in the Frequency Domain (DFT) applies Circular Convolution while you're after Linear Convolution. For practical examples and more information have a look on my answers: Kernel Convolution in Frequency Domain - Cyclic Padding. Question About Linear and Circular Convolution - 1D and 2D.


1

I work with these types of transforms all the time. My tool of choice is not the full FFT, but rather a discrete Fourier sine transform (DST). Not all Fourier analysis packages have one, but it looks like NAG's C06RAF will do the trick. The example code supplied by NAG should be a good starting point, but here's the outline I follow: Make sure my $q$ and $r$...


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