1
$\begingroup$

Assume I have a continuous function $f(x)$. Its spectral representation is $F(X)=\mathbb{F}(f)$. I need to calculate another function $f'(x)$ which has a spectral representation of $F'(X)=g(F(x))$, so basically $f'(x)={\mathbb{F}}^{-1}(g(\mathbb{F}(f)))$, where $\mathbb{F}$ is a fourier transform, $\mathbb{F}^{-1}$ is inverse fourier transform. $g()$ is a relatively well-behaved smooth function, but possibly nonlinear.

Now, in continuous case it's simple : just calculate $F(X)$ by doing Fourier transform of $f$, calculate $g(F(X))$, calculate $f'$ using inverse Fourier transform. However, for practical purposes f needs to be discretized into, say, $N$ points. And the spectrum is then $N$ harmonics. And after applying $g()$ to each harmonic, the result can produce a non-integer harmonic, basically. For example, assume I have 10 frequencies with 10Hz between 2 frequencies. And my $g(X)=1.5X$, in this case a frequency of 30Hz would map onto 45Hz that does not "align" with any other frequency in the spectrum.

What I tried to do is to do some sort of interpolation: for each frequency $X$, I calculate the frequency that transforms into it using $g^{-1}(X)$ and then interpolate that in the frequency domain. However, the exact way of interpolating a frequency has proven to be crucial because even linear interpolation vs sinc interpolation produce distinctly different results.

What I noticed is that sinc interpolation in spectral domain produces correct result, but it applies a window to the initial signal with the shape of the discretization area. As in, it acts as if the initial signal is not periodic (which is actually desirable in my case), but this window is somehow shifted and since I don't know why it originates to start with, I don't know how to control it.

$\endgroup$
1
$\begingroup$

I have actually came up with my own solution for this problem. I'm not sure it's the best way to do it, but I have deduced it from first principles, so I thought maybe it's worth sharing or getting a feedback on it. Basically the main idea is as follows:

$f'$ is a sum of its harmonics $F'(X)$, so $f'(x)=\frac{1}{M}\sum_{m}X'_m e^{i2\pi mx}$. Each harmonic $m$ is a sine wave that has an amplitude $X'_m$ that can be represented as $$X'_m=\int_0^1 e^{-i2\pi mx}f'(x)dx\ (1)$$, but $f'(x)$ in turn can also be reconstructed from $X_n$ via: $$f'(x)=\frac{1}{N}\sum_n X_n e^{i2\pi g(x)n}\ (2)$$ where function $g()$ is the transformation of the basis function $e^{i2\pi xk}$ from the non-primed basis to the primed one. In my case, this transformation is analytical and I can excactly transform each harmonic with $g()$.

Now I can substitute $(2)$ into $(1)$: $$X'_m=\frac{1}{N}\sum_n X_n\int_0^1 e^{-i2\pi mx} e^{i2\pi g(x)n}dx\ (3)$$

So far these formulas are exact and the only error introduced is from discretization of both $f()$ and $f'()$ on a finite set of points. However, the summation in (3) has a computational cost of $O(N)$ per harmonic, which is prohibitively expensive. But the cool part that I found is that only a few values of $n$ will have a significant contribution.

So in order to accelerate this step, I followed this rationale: if a frequency number $n$ on non-primed basis is transformed exactly into some other integer frequency number $m$ of the primed basis, then their convolution is exactly 1 and all other convolutions will be exactly zero. But often an integer harmonic $m$ will transform to a non-integer $n=g(m)$, for example if a sin wave from bin 1 shifts into a bin 4.2f, you can expect it to have the most contribution from bin 4, bin 5 and maybe some other bins in their vicinity. So in practice it suffices to limit the summation from the entire $\sum_N$ to summation over a small vicinity of $g(n)$, for example by using floor or rounding function: $[floor(g(n))-\Delta..floor(g(n))+\Delta]$ where $\Delta$ is the radius of summation. This range becomes the entire range $[0..N-1]$ if $\Delta=N/2$, but in my experiments $Delta$ of 2 to 10 has proven to be perfectly enough.

This method also totally works for $N \neq M$, so it can be used for transformations involving re-sampling. This formula also works neatly for both finite and infinite bases as long as the integral is calculated correctly for finite basis functions defined for example in a box: just make sure that box of the non-primed basis is transformed for each harmonic according to its $g()$. This allows to naturally implement transformations such as shifting and scaling without periodic repetition of the non-primed signal.

Overall cost of the algorithm is $O(FFT_N)+O(IFFT_M)+O(M\Delta)$ so as long as $\Delta < \log(M)$, this additional convolution step won't increase the asymptotic cost of the overall algorithm.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.