5
$\begingroup$

I have a power spectrum $P(\omega) = 10/\omega^2$ and I want to simulate a random time-series realization of this process $f(t)$ and a time-shifted copy of that realization $f(t- \tau)$. To do this, I want to the Fourier Shift theorem $$\mathcal{F}[f(t - \tau)](\omega) = e^{-2\pi j \omega \tau}\mathcal{F}[f(t)](\omega)$$

To get $f(t)$ I first discretize the frequency domain, get the amplitude of the frequency components by taking $A(\omega) = \sqrt{2P(\omega)\Delta \omega}$ and then multiply them by a random complex phase uniformly spaced on $[0, 2\pi]$. Finally I perform the Fourier Transform and take the real part to get a random realization of $f(t)$. This overall looks something like $$ f(m\Delta t) = Re\left\{\sum_{n=1}^N A(\omega_n)e^{j n m 2\pi/N} e^{j\theta}\right\} $$ where $\theta = \mathcal{U}[0,2\pi]$.

To get the shifted copy, I use the shift theorem and get

$$ f(m\Delta t - \tau) = Re\left\{\sum_{n=1}^N A(\omega_n)e^{j n m 2\pi/N} e^{j\theta}e^{-2\pi \omega_n \tau}\right\} $$ where $\theta = \mathcal{U}[0,2\pi]$.

P[w_] := 10^1 w^-2
tot = 2^10;
wmax = .1;
dw = wmax/tot;
dt = N[(2 \[Pi])/(2 dw tot)]
An = Table[2 Sqrt[P[n dw] dw], {n, 1, tot}];
theta = RandomVariate[UniformDistribution[{0, 2 \[Pi]}], tot];
Bn = Exp[I theta];
Cn = Table[Exp[-2 \[Pi] I n dw dt 20], {n, 1, tot}];

f1 = Re[Fourier[An*Bn, FourierParameters -> {1, -1}]];
f2 = Re[Fourier[An*Bn*Cn, FourierParameters -> {1, -1}]];

ListPlot[{f1, RotateRight[f2, 20]}, Joined -> True]

Below is the output of this calculation, and not only do the time-series not line up after I've shifted them, but there are also differences in the amplitude. I would like help in understanding why there are differences in amplitude, why the shifts don't line up, and how to correctly apply the Fourier Shift Theorem.

Two time-series, one shifted relative to the other using the Fourier Shift Theorem. The time-series do not overlap as expected and there are amplitude changes. Y-coordinate is amplitude, X-coordinate is time.

$\endgroup$

1 Answer 1

2
$\begingroup$

This is all surprisingly tricky, so I show what I would do, rather than attempting to fix your version.

A number of general points:

  • You have to use consistent times and frequencies. Here I have unit spacing in time and consequently a total frequency range of 1. Note that the frequencies appear in the spectrum in a shifted order.
  • You (may) need to worry about aliasing. Until you know what you are doing, ensure that your spectra drop to zero rapidly at the frequency discontinuity (where the minimum and maximum frequency meet)

Define the parameters of the problem, and a power spectrum that drops off rapidly.

npoints = 256;
bandwidth = 1/24;
spectrum = Exp[-π (f/bandwidth/2)^2];

Define a function that moves the zero frequency to the edge of the data set (or back again, assuming npoints is even).

fftshift = RotateLeft[#, Length[#]/2] &;

It is important that we define the frequency of our spectral values correctly

framp = fftshift[Range[-npoints/2, npoints/2 - 1]/npoints];

Use complex Gaussian noise

noise = {1, I} . RandomReal[NormalDistribution[0, 1], {2, npoints}];

and shape its spectrum

spec1 = noise*Sqrt[spectrum /. f -> framp];

ListLinePlot[Transpose[{fftshift[framp], fftshift[Abs[spec1]^2]}], 
 PlotRange -> All]

Power Spectrum

Define a signal with this power spectrum

signal1 = InverseFourier[spec1, FourierParameters -> {1, -1}];

ListLinePlot[Transpose[ReIm[signal1]]]

Noise with required spectrum

(You could check its autocorrelation function, but I've omitted that step).

Define a delay

dt = 7.5;

and apply this delay in the frequency domain

phaseterm = Exp[I*2*Pi*framp*dt];
spec2 = spec1*phaseterm;

Create a delayed version of the signal

signal2 = InverseFourier[spec2, FourierParameters -> {1, -1}];

and see that (the real part) is indeed delayed

ListLinePlot[Re /@ {signal1, signal2}]

With and without delay

Zoom in to verify that the delay has the requested value

ListLinePlot[Re /@ {signal1, signal2}, PlotRange -> {{150, 180}, All}]

Zoomed

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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